Rev 329 | Rev 331 | 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> |
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; |
193 | f9daq | 43 | int fData; |
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); |
319 | 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 | |||
267 | f9daq | 70 | Blab2::Blab2():fNfirst(0), fNeve(0), fData(0), 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.)", |
||
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 | |||
114 | int Blab2::Fill(std::vector<int> id, BParticle *p){ |
||
115 | for (int i=0; i< id.size(); i++){ |
||
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 |
||
193 | f9daq | 164 | BParticle p = *p1 + *p2; // Combine two particles into a new particle |
165 | if (p.InMassRange(min, max)){ |
||
272 | 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 | |||
188 | fList[_p2]->Clear(); |
||
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 |
||
198 | BParticle p = *p1 + *p2 + *p3; // Combine two particles into a new particle |
||
199 | if (p.InMassRange(min, max)){ |
||
200 | Fill(hid, &p); |
||
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 ); |
||
230 | Fill(hid, p); |
||
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 | |||
272 | TFile * f = new TFile(TString::Format("../../data/hadron-%d.root",fData)); // Open a data file |
||
273 | if(f->IsZombie()) { send_message(0,TString::Format("File %d not found\n",fData), 0 ); return; } |
||
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 |