Subversion Repositories f9daq

Rev

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