]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/tzero15tracks.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Destructors fixed to check for already existing objects before deleting (F.Bellini)
[u/mrichter/AliRoot.git] / TOF / tzero15tracks.C
CommitLineData
bcb56cb9 1void 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
2672] << 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
2692] << 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
271icle[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
314ified14);
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