Subversion Repositories f9daq

Rev

Rev 333 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
193 f9daq 1
#include <TClonesArray.h>
2
#include <TH1F.h>
3
#include <TStopwatch.h>
4
#include <TDatime.h>
5
#include <TString.h>
6
#include <TFile.h>
7
#include <TTree.h>
8
#include <TBranch.h>
272 f9daq 9
#include <TBufferJSON.h>
273 f9daq 10
#include <TMath.h>
319 f9daq 11
#include <vector>
193 f9daq 12
#include "BParticle.h"
13
#include "BEvent.h"
14
 
15
class Hdr{
16
public:
17
  int id;
18
  int len;
19
  int progress;
20
};
21
 
319 f9daq 22
std::vector<int> histogram(int n, ...){
23
std::vector<int> result;
24
int val = 0;
25
   va_list ap;
26
   int i;
193 f9daq 27
 
319 f9daq 28
   va_start(ap, n);
29
   for(i = 0; i < n; i++) {
30
      result.push_back(  va_arg(ap, int) );
31
   }
32
   va_end(ap);
33
return result;
34
}
35
 
193 f9daq 36
class Blab2 {
37
public:
329 f9daq 38
const char *names[12]={"photon", "electron", "pion", "muon", "kaon", "proton", "J/Psi", "D", "D*", "B", "Phi","Lambda0"};
272 f9daq 39
 
193 f9daq 40
UInt_t fNeve;
266 f9daq 41
UInt_t fNfirst;
267 f9daq 42
UInt_t fPrint;
344 f9daq 43
TString fData;
193 f9daq 44
TH1F *fH[100];
266 f9daq 45
UInt_t fHtype[100];
193 f9daq 46
TClonesArray *fList[100];
47
 
48
Blab2();
49
~Blab2();
50
void Init();
51
void event();
52
void Process();
53
void h1d(const char *varname, const char *name, int nbins, double min, double max, int id );
54
int  selector(int pin, int charge, SIMPLEPID particlename,  int hid, int pout );
319 f9daq 55
int  selector(int pin, int charge, SIMPLEPID particlename,  std::vector<int> hid, int pout );
267 f9daq 56
int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, int hid, int id );
319 f9daq 57
int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
330 f9daq 58
int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, int hid, int id );
59
int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
193 f9daq 60
int  fix_mass(int id);
344 f9daq 61
int  Fill(std::vector<int> hid, BParticle *p);
193 f9daq 62
void plist(int i);
63
 
64
 
65
ClassDef ( Blab2, 1 )
66
};
67
 
68
ClassImp(Blab2)
69
 
344 f9daq 70
Blab2::Blab2():fNfirst(0), fNeve(0), fData(), fPrint(0) {
193 f9daq 71
 
72
 Process();
73
};
74
 
266 f9daq 75
 
193 f9daq 76
void Blab2::h1d(const char *varname, const char *name, int nbins, double min, double max, int id ){
266 f9daq 77
   TString svar(varname);
273 f9daq 78
   TString axis[]={"mass (GeV/c2)",
79
		"momentum (GeV/c)",
80
		"energy (GeV)","charge",
81
		"identity",
82
		"momentum (GeV/c)",
83
		"momentum (GeV/c)",
84
		"momentum (GeV/c)",
85
		"momentum (GeV/c)",
86
		"angle (deg.)",
344 f9daq 87
		"cos(theta)"};
266 f9daq 88
   fHtype[id] = 0;
89
   if (svar.Contains("GetMass"    )) fHtype[id]=0;
90
   if (svar.Contains("GetMomentum")) fHtype[id]=1;
91
   if (svar.Contains("GetEnergy"  )) fHtype[id]=2;
272 f9daq 92
   if (svar.Contains("GetCharge"  )) fHtype[id]=3;
93
   if (svar.Contains("GetPid"     )) fHtype[id]=4;
94
   if (svar.Contains("GetXMomentum")) fHtype[id]=5;
95
   if (svar.Contains("GetYMomentum")) fHtype[id]=6;
96
   if (svar.Contains("GetZMomentum")) fHtype[id]=7;
277 f9daq 97
   if (svar.Contains("GetTransverseMomentum")) fHtype[id]=8;
273 f9daq 98
   if (svar.Contains("GetTheta"))             fHtype[id]=9;
99
   if (svar.Contains("GetCosTheta"))          fHtype[id]=10;
266 f9daq 100
 
319 f9daq 101
   //fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
272 f9daq 102
   if (fHtype[id]==4) {
103
     fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), 11, 0, 11);
104
     for (int i=0;i<11;i++) fH[id]->GetXaxis()->SetBinLabel(i+1,names[i]);
105
   } else {
106
     fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
107
   }
108
 
109
 
193 f9daq 110
}
111
 
319 f9daq 112
 
113
 
344 f9daq 114
int Blab2::Fill(std::vector<int> id, BParticle *p){
319 f9daq 115
  for (int i=0; i< id.size(); i++){
344 f9daq 116
 int hid = id[i];
273 f9daq 117
  if (hid>=0 && fH[hid]) {
118
      double val;
266 f9daq 119
      switch (fHtype[hid]){
273 f9daq 120
      case 0 : val  = p->GetMass(); break;
121
      case 1 : val  = p->GetMomentum(); break;
122
      case 2 : val  = p->e(); break;
123
      case 3 : val  = p->charge(); break;
124
      case 4 : val  = p->pid(); break;
125
      case 5 : val  = p->px(); break;
126
      case 6 : val  = p->py(); break;
127
      case 7 : val  = p->pz(); break;
128
      case 8 : val  = p->GetTransverseMomentum(); break;
277 f9daq 129
      case 9 : val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; val = 180.0*TMath::ACos(val)/TMath::Pi(); break;
273 f9daq 130
      case 10: val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; break;
131
      default: val  = 0 ; break;
132
   }
266 f9daq 133
   fH[hid]->Fill(val);
134
}
319 f9daq 135
 
136
   }
266 f9daq 137
 
138
return 0;
139
}
319 f9daq 140
int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, int hid, int _p2 ){
141
std::vector<int> a;
142
return combiner(_p0,_p1,same,pid,min,max,a,_p2);
143
}
266 f9daq 144
 
330 f9daq 145
int Blab2::combiner3(int _p0, int _p1, int _p2, int same, SIMPLEPID pid, double min, double max, int hid, int _p3 ){
146
std::vector<int> a;
147
return combiner3(_p0,_p1,_p2, same,pid,min,max,a,_p3);
148
}
149
 
319 f9daq 150
int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, std::vector<int> hid, int _p2 ){
193 f9daq 151
   // Loop over all the particles in both lists.
152
 if (_p0 < 0 ) _p0 =0;
153
 if (_p1 < 0 ) _p1 =0;
154
 
155
 
156
 fList[_p2]->Clear();
157
 int nprt=0;
158
 
159
 for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
160
    // the second loop
267 f9daq 161
   // in the case the second parti
162
   for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {
163
      if (p1==p2) continue;     // do not use the same particle in the combinations
344 f9daq 164
      BParticle  p = *p1 + *p2; // Combine two particles into a new particle
193 f9daq 165
      if (p.InMassRange(min, max)){
344 f9daq 166
            Fill(hid, &p);
267 f9daq 167
 	    p.SetPid(pid); // set PID to particlename to fix the particle mass
272 f9daq 168
            p.SetEnergyFromPid();
169
	    TClonesArray& list = *fList[_p2];
193 f9daq 170
	    new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
272 f9daq 171
 
193 f9daq 172
      }
173
 
174
   }
175
 
176
 }
177
 return _p2;
178
}
179
 
180
 
330 f9daq 181
int Blab2::combiner3(int _p0, int _p1,int _p2, int same, SIMPLEPID pid, double min, double max, std::vector<int> hid, int _p3 ){
182
   // Loop over all the particles in both lists.
183
 if (_p0 < 0 ) _p0 =0;
184
 if (_p1 < 0 ) _p1 =0;
185
 if (_p2 < 0 ) _p2 =0;
186
 
187
 
331 f9daq 188
 fList[_p3]->Clear();
330 f9daq 189
 int nprt=0;
190
 
191
 for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
192
    // the second loop
193
   // in the case the second parti
194
   for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {
195
      if (p1==p2) continue;     // do not use the same particle in the combinations
196
      for(TIter next3 = (_p1!=_p2 && same==0) ?  TIter(fList[_p2]): TIter(next2) ; BParticle * p3 =(BParticle *) next3();) {
197
        if (p2==p3) continue;     // do not use the same particle in the combinations
344 f9daq 198
        BParticle  p = *p1 + *p2 + *p3; // Combine two particles into a new particle
330 f9daq 199
        if (p.InMassRange(min, max)){
344 f9daq 200
            Fill(hid, &p);
330 f9daq 201
 	    p.SetPid(pid); // set PID to particlename to fix the particle mass
202
            p.SetEnergyFromPid();
203
	    TClonesArray& list = *fList[_p3];
204
	    new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
205
        }
206
      }
207
 
208
   }
209
 
210
 }
211
 return _p3;
212
}
213
 
214
 
319 f9daq 215
int Blab2::selector(int pin, int charge, SIMPLEPID type ,  int  hid, int pout ){
216
std::vector<int> a;
217
return selector(pin,charge,type,a,pout);
218
}
219
int Blab2::selector(int pin, int charge, SIMPLEPID type ,  std::vector<int> hid, int pout ){
193 f9daq 220
 if (pin < 0 ) pin =0;
221
 
222
  fList[pout]->Clear();
223
  int nprt=0;
224
 
225
  for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();) {
272 f9daq 226
	if (p->charge()== charge || charge > 1){
227
          if ( p->pid()== type || type == ALL ) {
228
	    TClonesArray& list = *fList[pout];
229
	    new (list[nprt++]) BParticle ( *p );
344 f9daq 230
            Fill(hid, p);
272 f9daq 231
          }
193 f9daq 232
	}
233
  }
234
   return pout;
235
}
236
 
237
 
238
int Blab2::fix_mass(int pin){
239
   if (pin < 0 ) pin =0;
240
   for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();)  p->SetEnergyFromPid();
241
   return pin;
242
}
243
 
244
void Blab2::plist(int i){
245
  fList[i]= new TClonesArray( "BParticle", 500 );
246
}
247
Blab2::~Blab2(){};
248
 
249
 
250
void send_message(int id, TString msg, int progress){
251
static Hdr hdr;
252
 
253
   hdr.id = id;
254
   hdr.len= msg.Length();
255
   hdr.progress= progress;
256
   fwrite(&hdr,3*sizeof(int),1,stdout);
257
   fwrite(msg.Data(),hdr.len,1,stdout);
258
   fflush(stdout);
259
 
260
}
261
 
262
 
263
void Blab2::Process(){
264
 
267 f9daq 265
char sList[0xFFFF];
193 f9daq 266
for (int i=0;i<100;i++) fH[i]=NULL;
266 f9daq 267
for (int i=0;i<100;i++) fHtype[i]=0;
193 f9daq 268
for (int i=0;i<100;i++) fList[i]=NULL;
269
 
270
Init();
271
 
344 f9daq 272
TFile * f = new TFile(TString::Format("../../data/%s",fData.Data())); // Open a data file
273
if(f->IsZombie()) {  send_message(0,TString::Format("File %s not found\n",fData.Data()), 0 );  return; }
193 f9daq 274
TTree * t =(TTree *) f-> Get( "T"); // Obtain a pointer to a tree of "event" data in the file
275
BEvent * mevent = new BEvent(); // Create a  "BEvent" object where data from the file will be loaded
276
TBranch * branch = t-> GetBranch( "BEvent"); // Obtain a branch for "BEvent" in the tree
277
branch-> SetAddress(&mevent); // Register created "BEvent" object to the branch
266 f9daq 278
TH1F *fHnprt= new TH1F("h100", "Number of particles in the event;N particles;N events", 50, -0.5, 49.5);
279
 
272 f9daq 280
 
281
 
267 f9daq 282
send_message(0, TString::Format("<br>Number of Events in the file %lld<br>\n", t->GetEntries() ),0);
193 f9daq 283
TStopwatch timer;
284
timer.Start();
267 f9daq 285
int nev  = 0;
266 f9daq 286
int i    =TMath::Min(fNfirst, (UInt_t) t-> GetEntries());
287
int cNeve=TMath::Min(fNfirst+fNeve, (UInt_t) t-> GetEntries());
273 f9daq 288
int fPart = fPrint;
289
int totaltracks = 0;
193 f9daq 290
while (i<cNeve){
291
 t-> GetEntry(i); // Read the content of the event from the file
292
 fList[0]= mevent->GetParticleList();
266 f9daq 293
 
193 f9daq 294
 event();
295
 
267 f9daq 296
 int progress = (100*i)/cNeve;
297
 if (i%10000==0) send_message(2, TString::Format("Event %d\n",i), progress);
266 f9daq 298
 
299
 int nprt=0;
267 f9daq 300
 if (nev>100) fPrint = 0; // disable particle prints for huge numer of events
301
 if (fPrint) sprintf(sList,"Primary particle list for Event %d<br/><table class='plist' ><tr><th>N<th>px(GeV/c)<th>py(GeV/c)<th>pz(GeV/c)<th>p(GeV/c)<th>Energy(GeV)<th>Charge<th>ID<th></tr>", i);
266 f9daq 302
 for(TIter next(fList[0]); BParticle * p =(BParticle *) next();) {
303
   nprt++;
267 f9daq 304
   if (fPrint) sprintf(sList,"%s<tr><td>%d<td>%g<td>%g<td>%g<td>%g<td>%g<td>%1.0f<td>%s</tr>",sList,nprt, p->px(),p->py(),p->pz(),p->GetMomentum(),p->e(), p->charge(),names[p->pid()] );
266 f9daq 305
 }
273 f9daq 306
 if (fPart) fHnprt->Fill(nprt);
307
 totaltracks += nprt;
267 f9daq 308
 if (fPrint) {
309
   sprintf(sList,"%s</table>",sList);
310
   send_message(0, TString(sList),progress);
311
   nev++;
312
 }
193 f9daq 313
 mevent-> Clear();  // Clear the memory.
314
 for (int k=0;k<100;k++) if (fList[k]!=0) fList[k]->Clear();
315
 i++;
316
}
273 f9daq 317
double avgtracks=(i)?float(totaltracks)/i:0;
318
send_message(0, TString::Format("Number of events processed: %d<br/>\nNumber of particles: %d<br/>\nAverage number of particles per event %f<br/>\n", i, totaltracks, avgtracks ),100);
193 f9daq 319
 
273 f9daq 320
if (fPart) send_message(1,TBufferJSON::ConvertToJSON(fHnprt),100 );
268 f9daq 321
 
322
 
323
for (int i=0;i<100;i++) if (fH[i]!=0) send_message(1,TBufferJSON::ConvertToJSON(fH[i]),100 );
324
 
193 f9daq 325
TDatime d;
326
timer.Stop();
327
send_message(3, TString::Format("'%s', %d, %f, %f", d.AsSQLString(),i, timer.RealTime(), timer.CpuTime() ),100);
328
 
329
}
330
 
331