Subversion Repositories f9daq

Rev

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