Rev 273 | Rev 319 | 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; |
||
277 | f9daq | 80 | if (svar.Contains("GetTransverseMomentum")) 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; |
||
277 | f9daq | 109 | case 9 : val = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; val = 180.0*TMath::ACos(val)/TMath::Pi(); break; |
273 | f9daq | 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 |