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