]>
Commit | Line | Data |
---|---|---|
bcb56cb9 | 1 | void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10) |
2 | { | |
3 | ||
4 | // | |
5 | // Calculate the t0 distribution | |
6 | // | |
7 | // (primary particles reaching TOF) | |
8 | // | |
9 | ||
10 | ||
11 | // Dynamically link some shared libs | |
12 | if (gClassTable->GetID("AliRun") < 0) { | |
13 | gROOT->LoadMacro("loadlibs.C"); | |
14 | loadlibs(); | |
15 | } | |
16 | ||
17 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
18 | Int_t ngood; | |
19 | Int_t nmisidentified=0; | |
20 | Int_t nmisidentified0=0; | |
21 | Int_t nmisidentified1=0; | |
22 | Int_t nmisidentified2=0; | |
23 | Int_t nmisidentified3=0; | |
24 | Int_t nmisidentified4=0; | |
25 | Int_t nmisidentified5=0; | |
26 | Int_t nmisidentified6=0; | |
27 | Int_t nmisidentified7=0; | |
28 | Int_t nmisidentified8=0; | |
29 | Int_t nmisidentified9=0; | |
30 | Int_t nmisidentified10=0; | |
31 | Int_t nmisidentified11=0; | |
32 | Int_t nmisidentified12=0; | |
33 | Int_t nmisidentified13=0; | |
34 | Int_t nmisidentified14=0; | |
35 | Int_t nbytes = 0; | |
36 | Int_t ipartold = -1; | |
37 | Int_t dblehit = 0; | |
38 | Int_t j,hit,ipart; | |
39 | Int_t nhits; | |
40 | Int_t selected=0; | |
41 | Int_t istop=0; | |
42 | const Float_t timeresolution=globalTimeRes; // in [s] | |
43 | Float_t timeresolutioninns=timeresolution*(1.e+9); // convert in [ns] | |
44 | cout << "Global Time Resolution " << timeresolutioninns << " ns" << endl; // | |
45 | ||
46 | // output (histos!) filename | |
47 | char outFileName[100]; | |
48 | strcpy(outFileName,"ht015tr"); | |
49 | strcat(outFileName,datafile); | |
50 | ||
51 | gBenchmark->Start("TOFT0"); | |
52 | ||
53 | TH1F* htzerobest= new TH1F("htzerobest","T0 for best assignment",200,-1.,1.); | |
54 | TH1F* hchibest = new TH1F("hchibest","ChiSquare Min Distribution",80,0.,40.); | |
55 | TH1F* hchibestconflevel = new TH1F("hchibestconflevel","ChiSquare Min Confidence Level",10,0.,1.); | |
56 | ||
57 | Float_t t0best=999.; | |
58 | Int_t assparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}; | |
59 | Int_t truparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}; | |
60 | Float_t timeofflight[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
61 | Float_t momentum[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
62 | Float_t timezero[15]; | |
63 | Float_t weightedtimezero[15]; | |
64 | Float_t beta[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
65 | Float_t sqMomError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
66 | Float_t sqTrackError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
67 | Float_t massarray[3]={0.13957,0.493677,0.9382723}; | |
68 | Float_t chisquare=999.; | |
69 | Float_t tracktoflen[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
70 | TRandom* rnd = new TRandom(); | |
71 | TParticle *particle; | |
72 | AliTOFhitT0 *tofHit; | |
73 | ||
74 | TFile *file = TFile::Open(datafile,"old"); | |
75 | ||
76 | // Get AliRun object from file or create it if not on file | |
77 | gAlice = (AliRun*)file->Get("gAlice"); | |
78 | ||
79 | // Import the Kine and Hits Trees for the event evNumber in the file | |
80 | ||
81 | Int_t nallevents = (Int_t) gAlice->TreeE()->GetEntries(); | |
82 | ||
83 | for (Int_t nev=0; nev<nallevents;nev++) { | |
84 | ||
85 | Int_t nparticles = gAlice->GetEvent(nev); | |
86 | if (nparticles <= 0) return; | |
87 | // ngood=0; | |
88 | // Get pointers to Alice detectors and Hits containers | |
89 | ||
90 | AliDetector *TOF = gAlice->GetDetector("TOF"); | |
91 | TObjArray *Particles = gAlice->Particles(); | |
92 | ||
93 | if (TOF) TClonesArray *TOFhits = TOF->Hits(); | |
94 | ||
95 | TTree *TH = gAlice->TreeH(); | |
96 | Int_t ntracks = TH->GetEntries(); | |
97 | ||
98 | ||
99 | // Start loop on primary tracks in the hits containers | |
100 | Int_t lasttrack=-1; | |
101 | Int_t nset=0; | |
102 | for (Int_t track=lasttrack+1; track<ntracks;track++) { | |
103 | ||
104 | if(nset>=3) break; // check on the number of set analyzed | |
105 | ||
106 | gAlice->ResetHits(); | |
107 | nbytes += TH->GetEvent(track); | |
108 | ||
109 | // =======>Histogram TOF | |
110 | if (TOF) { | |
111 | nhits = TOFhits->GetEntriesFast(); | |
112 | ||
113 | for (hit=0;hit<nhits;hit++) { | |
114 | ||
115 | tofHit = (AliTOFhitT0*)TOFhits->UncheckedAt(hit); | |
116 | ipart = tofHit->GetTrack(); | |
117 | ||
118 | // check to discard the case when the same particle is selected more than one | |
119 | // time | |
120 | ||
121 | if (ipart != ipartold){ | |
122 | ||
123 | particle = (TParticle*)gAlice->Particle(ipart); | |
124 | ||
125 | Float_t idealtime=tofHit->GetTof(); | |
126 | // Float_t time=idealtime; | |
127 | Float_t time = rnd->Gaus(idealtime, timeresolution); | |
128 | Float_t toflen=tofHit->GetLen(); | |
129 | toflen=toflen/100.; // toflen given in m because GEANT gives cm! | |
130 | Int_t pdg = particle->GetPdgCode(); | |
131 | Int_t abspdg =TMath::Abs(pdg); | |
132 | Float_t idealmom = particle->P(); | |
133 | Float_t momres=idealmom*0.025; // 2.5% res token into account for all momenta | |
134 | //Float_t mom =idealmom+2.*(rnd->Rndm()-0.5)*momres; | |
135 | Float_t mom =rnd->Gaus(idealmom,momres); | |
136 | Bool_t isgoodpart=(abspdg==211 || abspdg==2212 || abspdg==321); | |
137 | ||
138 | time*=1.e+9; // tof given in nanoseconds | |
139 | if (particle->GetFirstMother() < 0 && isgoodpart && mom<=1.7 && mom>=1.25){ | |
140 | ||
141 | selected+=1; | |
142 | istop=selected; | |
143 | if(istop>15) break; | |
144 | Int_t index=selected-1; | |
145 | timeofflight[index]=time; | |
146 | tracktoflen[index]=toflen; | |
147 | momentum[index]=mom; | |
148 | switch (abspdg) { | |
149 | case 211: | |
150 | truparticle[index]=0; | |
151 | break ; | |
152 | case 321: | |
153 | truparticle[index]=1; | |
154 | break ; | |
155 | case 2212: | |
156 | truparticle[index]=2; | |
157 | break ; | |
158 | } | |
159 | ||
160 | } | |
161 | ipartold = ipart; | |
162 | ||
163 | if(istop==15){ // start analysis on current set | |
164 | nset+=1; | |
165 | lasttrack=track; | |
166 | istop=0; | |
167 | selected=0; | |
168 | cout << "starting t0 calculation for current set" << endl; | |
169 | for (Int_t i1=0; i1<3;i1++) { | |
170 | for (Int_t i2=0; i2<3;i2++) { | |
171 | for (Int_t i3=0; i3<3;i3++) { | |
172 | for (Int_t i4=0; i4<3;i4++) { | |
173 | for (Int_t i5=0; i5<3;i5++) { | |
174 | for (Int_t i6=0; i6<3;i6++) { | |
175 | for (Int_t i7=0; i7<3;i7++) { | |
176 | for (Int_t i8=0; i8<3;i8++) { | |
177 | for (Int_t i9=0; i9<3;i9++) { | |
178 | for (Int_t i10=0; i10<3;i10++) { | |
179 | for (Int_t i11=0; i11<3;i11++) { | |
180 | for (Int_t i12=0; i12<3;i12++) { | |
181 | for (Int_t i13=0; i13<3;i13++) { | |
182 | for (Int_t i14=0; i14<3;i14++) { | |
183 | for (Int_t i15=0; i15<3;i15++) { | |
184 | beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]); | |
185 | beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]); | |
186 | beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]); | |
187 | beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]); | |
188 | beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]); | |
189 | beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]); | |
190 | beta[6]=momentum[6]/sqrt(massarray[i7]*massarray[i7]+momentum[6]*momentum[6]); | |
191 | beta[7]=momentum[7]/sqrt(massarray[i8]*massarray[i8]+momentum[7]*momentum[7]); | |
192 | beta[8]=momentum[8]/sqrt(massarray[i9]*massarray[i9]+momentum[8]*momentum[8]); | |
193 | beta[9]=momentum[9]/sqrt(massarray[i10]*massarray[i10]+momentum[9]*momentum[9]); | |
194 | beta[10]=momentum[10]/sqrt(massarray[i11]*massarray[i11]+momentum[10]*momentum[10]); | |
195 | beta[11]=momentum[11]/sqrt(massarray[i12]*massarray[i12]+momentum[11]*momentum[11]); | |
196 | beta[12]=momentum[12]/sqrt(massarray[i13]*massarray[i13]+momentum[12]*momentum[12]); | |
197 | beta[13]=momentum[13]/sqrt(massarray[i14]*massarray[i14]+momentum[13]*momentum[13]); | |
198 | beta[14]=momentum[14]/sqrt(massarray[i15]*massarray[i15]+momentum[14]*momentum[14]); | |
199 | ||
200 | ||
201 | Float_t meantzero=0.; | |
202 | Float_t sumAllweights=0.; | |
203 | for (Int_t itz=0; itz<15;itz++) { | |
204 | sqMomError[itz]=((1.-beta[itz]*beta[itz])*0.025)*((1.-beta[itz]*beta[itz])*0.025)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz])); // this gives the square of the momentum error in nanoseconds | |
205 | sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); // total error for the current track | |
206 | sumAllweights+=1./sqTrackError[itz]; | |
207 | // redefining beta, it is useful in order to calculate t zero | |
208 | beta[itz]*=0.299792; | |
209 | timezero[itz]=(tracktoflen[itz]/beta[itz])-timeofflight[itz]; | |
210 | weightedtimezero[itz]=((tracktoflen[itz]/beta[itz])-timeofflight[itz])/sqTrackError[itz];// weighted time zero for current track | |
211 | meantzero+=weightedtimezero[itz]; | |
212 | } | |
213 | meantzero=meantzero/sumAllweights; // it is given in [ns] | |
214 | ||
215 | Float_t dummychisquare=0.; | |
216 | // calculate the chisquare for the current assignment | |
217 | for (Int_t itz=0; itz<15;itz++) { | |
218 | dummychisquare+=(timezero[itz]-meantzero)*(timezero[itz]-meantzero)/sqTrackError[itz]; | |
219 | } | |
220 | ||
221 | ||
222 | if(dummychisquare<=chisquare){ | |
223 | assparticle[0]=i1; | |
224 | assparticle[1]=i2; | |
225 | assparticle[2]=i3; | |
226 | assparticle[3]=i4; | |
227 | assparticle[4]=i5; | |
228 | assparticle[5]=i6; | |
229 | assparticle[6]=i7; | |
230 | assparticle[7]=i8; | |
231 | assparticle[8]=i9; | |
232 | assparticle[9]=i10; | |
233 | assparticle[10]=i11; | |
234 | assparticle[11]=i12; | |
235 | assparticle[12]=i13; | |
236 | assparticle[13]=i14; | |
237 | assparticle[14]=i15; | |
238 | chisquare=dummychisquare; | |
239 | t0best=meantzero; | |
240 | } | |
241 | } | |
242 | } | |
243 | } | |
244 | } | |
245 | } | |
246 | } | |
247 | } | |
248 | } | |
249 | } | |
250 | } | |
251 | ||
252 | } | |
253 | } | |
254 | } | |
255 | } | |
256 | } | |
257 | ||
258 | cout << "true" << truparticle[0] << truparticle[1] << truparticle[2] << truparticle[3] << truparticle[4] << truparticle[5] << truparticle[6] << truparticle[7] << truparticle[8] << truparticle[9] <<truparticle[10] << truparticle[11] << truparticle[12] << truparticle[13] << truparticle[14] <<endl; | |
259 | cout << "best" << assparticle[0] << assparticle[1] << assparticle[2] << assparticle[3] << assparticle[4] << assparticle[5] << assparticle[6] << assparticle[7] << assparticle[8] << assparticle[9] <<assparticle[10] << assparticle[11] << assparticle[12] << assparticle[13] << assparticle[14] << endl; | |
260 | if(truparticle[0]==assparticle[0] && truparticle[1]==assparticle[1] && truparticle[2]==assparticle[2] && truparticle[3]==assparticle[3] && truparticle[4]==assparticle[4]&& truparticle[5]==assparticle[5] && truparticle[6]==assparticle[6] && truparticle[7]==assparticle[7] && truparticle[8]==assparticle[8] && truparticle[9]==assparticle[9] &&truparticle[10]==assparticle[10] && truparticle[11]==assparticle[11] && truparticle[12]==assparticle[12] && truparticle[13]==assparticle[13] && truparticle[14]==assparticle[14]) ngood+=1; | |
261 | if(truparticle[0]!=assparticle[0]) nmisidentified0+=1; | |
262 | if(truparticle[1]!=assparticle[1]) nmisidentified1+=1; | |
263 | if(truparticle[2]!=assparticle[2]) nmisidentified2+=1; | |
264 | if(truparticle[3]!=assparticle[3]) nmisidentified3+=1; | |
265 | if(truparticle[4]!=assparticle[4]) nmisidentified4+=1; | |
266 | if(truparticle[5]!=assparticle[5]) nmisidentified5+=1; | |
267 | if(truparticle[6]!=assparticle[6]) nmisidentified6+=1; | |
268 | if(truparticle[7]!=assparticle[7]) nmisidentified7+=1; | |
269 | if(truparticle[8]!=assparticle[8]) nmisidentified8+=1; | |
270 | if(truparticle[9]!=assparticle[9]) nmisidentified9+=1; | |
271 | if(truparticle[10]!=assparticle[10]) nmisidentified10+=1; | |
272 | if(truparticle[11]!=assparticle[11]) nmisidentified11+=1; | |
273 | if(truparticle[12]!=assparticle[12]) nmisidentified12+=1; | |
274 | if(truparticle[13]!=assparticle[13]) nmisidentified13+=1; | |
275 | if(truparticle[14]!=assparticle[14]) nmisidentified14+=1; | |
276 | cout << "chisquare for current set" << chisquare << endl; | |
277 | //Float_t weight=1./chisquare; | |
278 | //htzerobest->Fill(t0best,weight); | |
279 | htzerobest->Fill(t0best); | |
280 | hchibest->Fill(chisquare); | |
281 | Double_t dblechisquare=(Double_t)chisquare; | |
282 | Float_t confLevel=(Float_t)TMath::Prob(dblechisquare,14); | |
283 | cout << "conf level " << confLevel << endl; | |
284 | hchibestconflevel->Fill(confLevel); | |
285 | chisquare=999.; | |
286 | t0best=999.; | |
287 | ||
288 | } // end for the current set. close if(istop==5) | |
289 | ||
290 | ||
291 | } // end condition on ipartold | |
292 | ||
293 | ||
294 | }// end loop on hits for the current track | |
295 | }// end TOF | |
296 | if(istop>=15) break; | |
297 | }// end loop on Tracks | |
298 | cout << "ngood for current event " << ngood << endl; | |
299 | ||
300 | ||
301 | }//end loop on events | |
302 | ||
303 | nmisidentified=(nmisidentified0+nmisidentified1+nmisidentified2+nmisidentified3+nmisidentified4+nmisidentified5+nmisidentified6+nmisidentified7+nmisidentified8+nmisidentified9+nmisidentified10+nmisidentified11+nmisidentified12+nmisidentified13+nmisidentified14); | |
304 | cout << "total misidentified " << nmisidentified << endl; | |
305 | TFile *houtfile = new TFile(outFileName,"recreate"); | |
306 | houtfile->cd(); | |
307 | ||
308 | htzerobest->Write(); | |
309 | hchibest->Write(); | |
310 | hchibestconflevel->Write(); | |
311 | houtfile->Close(); | |
312 | ||
313 | gBenchmark->Stop("TOFT0"); | |
314 | cout << "15 tracks:" << endl ; | |
315 | cout << " took " << gBenchmark->GetCpuTime("TOFT0") << " seconds in order to calculate T0 " | |
316 | << gBenchmark->GetCpuTime("TOFT0")/nallevents << " seconds per event " << endl ; | |
317 | cout << endl ; | |
318 | ||
319 | } | |
320 | ||
321 |