]>
Commit | Line | Data |
---|---|---|
06d06fbb | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | ||
17 | /////////////////////////////////////////////////////////////////////////////// | |
18 | // // | |
19 | // Time Projection Chamber // | |
20 | // Comparison macro for ESD // | |
21 | // responsible: | |
22 | // marian.ivanov@cern.ch // | |
23 | // | |
24 | // | |
25 | ||
26 | /* | |
27 | marian.ivanov@cern.ch | |
28 | Usage: | |
29 | ||
30 | ||
31 | ||
32 | gSystem->Load("libPWG1.so"); | |
33 | // | |
1a167b98 | 34 | AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0); |
06d06fbb | 35 | t2->Exec(); |
36 | ||
37 | ||
38 | TFile f("cmpESDTracks.root"); | |
06d06fbb | 39 | TTree * tree = (TTree*)f.Get("ESDcmpTracks"); |
1a167b98 | 40 | |
41 | AliTreeDraw comp; | |
06d06fbb | 42 | comp.SetTree(tree) |
43 | ||
44 | ||
45 | ||
46 | // | |
47 | //some cuts definition | |
48 | TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01") | |
49 | //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5") | |
50 | //TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9"); | |
51 | TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5"); | |
52 | TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5"); | |
53 | ||
54 | ||
55 | TCut crec("crec","fReconstructed==1"); | |
56 | TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25"); | |
57 | TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1"); | |
58 | ||
59 | TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1"); | |
60 | TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50"); | |
61 | TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)"); | |
62 | TCut cchi2("cchi2","fESDtrack.fITSchi2MIP[0]<7.&&fESDtrack.fITSchi2MIP[1]<5.&&fESDtrack.fITSchi2MIP[2]<7.&&fESDtrack.fITSchi2MIP[3]<7.5&&fESDtrack.fITSchi2MIP[4]<6.") | |
63 | ||
64 | ||
65 | // | |
66 | //example | |
c4d09608 | 67 | comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); |
68 | comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); | |
69 | comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); | |
70 | comp.T()->SetAlias("theta","MC.fTrackRef.Theta()"); | |
71 | comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); | |
72 | comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
73 | comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
06d06fbb | 74 | |
75 | ||
76 | TH1F his("his","his",100,0,20); | |
77 | TH1F hpools("hpools","hpools",100,-7,7); | |
78 | TH1F hfake("hfake","hfake",1000,0,150); | |
79 | TProfile profp0("profp0","profp0",20,-0.4,0.9) | |
80 | ||
81 | comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06) | |
82 | comp.fRes->Draw(); | |
83 | comp.fMean->Draw(); | |
84 | ||
85 | comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06) | |
86 | comp.fRes->Draw(); | |
87 | ||
88 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5) | |
89 | comp.fRes->Draw(); | |
90 | ||
91 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio<0.1",10,0.2,1.5) | |
92 | comp.fRes->Draw(); | |
93 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio>0.1",10,0.2,1.5) | |
94 | comp.fRes->Draw(); | |
95 | ||
c4d09608 | 96 | comp.T()->Draw("fESDtrack.fITSsignal/fESDtrack.fTPCsignal","fITSOn&&fTPCOn&&fESDtrack.fITSFakeRatio==0") |
06d06fbb | 97 | |
98 | TH1F his("his","his",100,0,20); | |
99 | TH1F hpools("hpools","hpools",100,-7,7); | |
100 | ||
101 | TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1); | |
102 | TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4); | |
103 | TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3); | |
104 | TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2); | |
105 | ||
c4d09608 | 106 | comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) |
107 | comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) | |
108 | comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) | |
109 | comp.T()->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) | |
06d06fbb | 110 | |
111 | ||
c4d09608 | 112 | comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") |
113 | comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") | |
114 | comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") | |
115 | comp.T()->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") | |
06d06fbb | 116 | |
c4d09608 | 117 | comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); |
118 | comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); | |
119 | comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); | |
120 | comp.T()->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1); | |
06d06fbb | 121 | |
122 | hdedx3->SetXTitle("P(GeV/c)"); | |
123 | hdedx3->SetYTitle("dEdx(unit)"); | |
124 | hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same"); | |
125 | ||
126 | comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8) | |
127 | ||
128 | TProfile prof("prof","prof",10,0.5,5); | |
129 | ||
130 | ||
131 | ||
132 | ||
133 | */ | |
134 | ||
135 | ||
136 | #include <stdio.h> | |
137 | #include <string.h> | |
138 | //ROOT includes | |
139 | #include "Rtypes.h" | |
140 | #include "TFile.h" | |
141 | #include "TTree.h" | |
06d06fbb | 142 | #include "TStopwatch.h" |
06d06fbb | 143 | #include "TVector3.h" |
d390cc7e | 144 | #include "TGeoManager.h" |
9f3282ed | 145 | //#include "Getline.h" |
1a167b98 | 146 | // |
06d06fbb | 147 | //ALIROOT includes |
1a167b98 | 148 | // |
06d06fbb | 149 | #include "AliRun.h" |
06d06fbb | 150 | #include "AliESDtrack.h" |
06d06fbb | 151 | #include "AliTPCParam.h" |
152 | #include "AliTPC.h" | |
06d06fbb | 153 | #include "AliTrackReference.h" |
06d06fbb | 154 | #include "AliTPCParamSR.h" |
155 | #include "AliTracker.h" | |
1a167b98 | 156 | #include "AliESDEvent.h" |
06d06fbb | 157 | #include "AliESD.h" |
158 | #include "AliESDfriend.h" | |
159 | #include "AliESDtrack.h" | |
160 | #include "AliTPCseed.h" | |
161 | #include "AliITStrackMI.h" | |
06d06fbb | 162 | #include "AliESDVertex.h" |
163 | #include "AliExternalTrackParam.h" | |
164 | #include "AliESDkink.h" | |
165 | #include "AliESDv0.h" | |
166 | #include "AliV0.h" | |
167 | // | |
168 | #include "AliTreeDraw.h" | |
76472f75 | 169 | #include "AliMCInfo.h" |
170 | #include "AliGenKinkInfo.h" | |
171 | #include "AliGenV0Info.h" | |
172 | ||
173 | ||
999d8278 | 174 | #include "AliESDRecInfo.h" |
175 | #include "AliESDRecV0Info.h" | |
176 | #include "AliESDRecKinkInfo.h" | |
06d06fbb | 177 | #include "AliRecInfoMaker.h" |
178 | ||
179 | ||
180 | ||
181 | ClassImp(AliRecInfoMaker) | |
182 | ||
183 | ||
184 | ||
185 | ||
186 | AliTPCParam * AliRecInfoMaker::GetTPCParam(){ | |
1a167b98 | 187 | // |
188 | // create TPC param | |
189 | // | |
06d06fbb | 190 | AliTPCParamSR * par = new AliTPCParamSR; |
191 | par->Update(); | |
192 | return par; | |
193 | } | |
194 | ||
195 | ||
196 | ||
52bbef74 | 197 | void AliRecInfoMaker::MakeAliases(TTree * tree) |
06d06fbb | 198 | { |
199 | // | |
200 | // aliases definition | |
201 | // | |
52bbef74 | 202 | tree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); |
203 | tree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); | |
204 | tree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); | |
205 | tree->SetAlias("theta","MC.fTrackRef.Theta()"); | |
206 | tree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); | |
207 | tree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
208 | tree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
06d06fbb | 209 | |
52bbef74 | 210 | tree->SetAlias("trddedx","(RC.fESDtrack.fTRDsignals[0]+RC.fESDtrack.fTRDsignals[1]+RC.fESDtrack.fTRDsignals[2]+RC.fESDtrack.fTRDsignals[3]+RC.fESDtrack.fTRDsignals[4]+RC.fESDtrack.fTRDsignals[5])/6."); |
06d06fbb | 211 | |
52bbef74 | 212 | tree->SetAlias("dtofmc2","fESDtrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)"); |
213 | tree->SetAlias("dtofrc2","(fESDtrack.fTrackTime[2]-fESDtrack.fTOFsignal)"); | |
214 | ||
215 | tree->SetAlias("psum","fESDtrack.fTOFr[4]+fESDtrack.fTOFr[3]+fESDtrack.fTOFr[2]+fESDtrack.fTOFr[1]+fESDtrack.fTOFr[0]"); | |
216 | tree->SetAlias("P0","fESDtrack.fTOFr[0]/psum"); | |
217 | tree->SetAlias("P1","fESDtrack.fTOFr[1]/psum"); | |
218 | tree->SetAlias("P2","fESDtrack.fTOFr[2]/psum"); | |
219 | tree->SetAlias("P3","fESDtrack.fTOFr[3]/psum"); | |
220 | tree->SetAlias("P4","fESDtrack.fTOFr[4]/psum"); | |
221 | tree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)"); | |
06d06fbb | 222 | } |
223 | ||
224 | ||
06d06fbb | 225 | //////////////////////////////////////////////////////////////////////// |
226 | AliRecInfoMaker::AliRecInfoMaker(const char* fnGenTracks, | |
cd875161 | 227 | const char* fnCmp, |
228 | const char* fnGalice, | |
229 | Int_t nEvents, Int_t firstEvent): | |
230 | ||
231 | fEventNr(0), //! current event number | |
232 | fNEvents(0), //! number of events to process | |
233 | fFirstEventNr(0), //! first event to process | |
234 | fFileCmp(0), //! output file with cmp tracks | |
235 | fTreeCmp(0), //! output tree with cmp tracks | |
236 | fTreeCmpKinks(0), //! output tree with cmp Kinks | |
237 | fTreeCmpV0(0), //! output tree with cmp V0 | |
238 | // | |
239 | fFileGenTracks(0), //! input files with generated tracks | |
240 | fTreeGenTracks(0), //! tree with generated tracks | |
241 | fTreeGenKinks(0), // tree with gen kinks | |
242 | fTreeGenV0(0), // tree with gen V0 | |
243 | // | |
244 | fLoader(0), //! pointer to the run loader | |
245 | // | |
246 | fIndexRecTracks(0), //! index of particle label in the TreeT_ESD | |
247 | fFakeRecTracks(0), //! number of fake tracks | |
248 | fMultiRecTracks(0), //! number of multiple reconstructions | |
249 | // | |
250 | fIndexRecKinks(0), //! index of particle label in treeesd | |
251 | fMultiRecKinks(0), //! number of multiple reconstructions | |
252 | fSignedKinks(0), //! indicator that kink was not fake | |
253 | // | |
254 | fIndexRecV0(0), //! index of particle label in treeesd | |
255 | fMultiRecV0(0), //! number of multiple reconstructions | |
256 | fSignedV0(0), //! indicator that kink was not fake | |
257 | // | |
258 | fRecArray(0), // container with rec infos | |
259 | fEvent(0), //!event | |
260 | fESDfriend(0), //!event friend | |
261 | // | |
262 | fParamTPC(0), //! AliTPCParam | |
263 | fNParticles(0), //! number of particles in the input tree genTracks | |
264 | fDebug(0), //! debug flag | |
265 | fNextTreeGenEntryToRead(0), //! last entry already read from genTracks tree | |
266 | fNextKinkToRead(0), //! last entry already read from genKinks tree | |
267 | fNextV0ToRead(0), //! last entry already read from genV0 tree | |
268 | // | |
269 | fMCInfo(0), //! MC information writen per particle | |
270 | fGenKinkInfo(0), //! MC information writen per Kink | |
271 | fGenV0Info(0), //! MC information writen per Kink | |
272 | fRecInfo(0), //! Rec. information writen per particle | |
273 | fFriend(0), //! friend track | |
274 | fRecKinkInfo(0), //! reconstructed kink info | |
275 | fRecV0Info(0) //! reconstructed kink info | |
06d06fbb | 276 | { |
1a167b98 | 277 | // AliRecInfoMaker - connencts the MC information with reconstructed information |
278 | // fnGenTracks - file with MC to be created before using AliGenInfoMaker | |
279 | // fnCmp - file name to be created | |
280 | // fnGalice - file with Loaders - usualy galice.root | |
281 | // | |
282 | // nEvent - number of event s to be analyzed | |
283 | // AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0); | |
284 | // | |
285 | ||
286 | ||
06d06fbb | 287 | Reset(); |
288 | // fFnGenTracks = fnGenTracks; | |
289 | // fFnCmp = fnCmp; | |
290 | sprintf(fFnGenTracks,"%s",fnGenTracks); | |
291 | sprintf(fFnCmp,"%s",fnCmp); | |
292 | ||
293 | fFirstEventNr = firstEvent; | |
294 | fEventNr = firstEvent; | |
295 | fNEvents = nEvents; | |
06d06fbb | 296 | // |
297 | fLoader = AliRunLoader::Open(fnGalice); | |
298 | if (gAlice){ | |
299 | //delete gAlice->GetRunLoader(); | |
300 | delete gAlice; | |
301 | gAlice = 0x0; | |
302 | } | |
303 | if (fLoader->LoadgAlice()){ | |
304 | cerr<<"Error occured while l"<<endl; | |
305 | } | |
306 | Int_t nall = fLoader->GetNumberOfEvents(); | |
307 | if (nEvents==0) { | |
308 | nEvents =nall; | |
309 | fNEvents=nall; | |
310 | fFirstEventNr=0; | |
311 | } | |
312 | ||
313 | if (nall<=0){ | |
314 | cerr<<"no events available"<<endl; | |
315 | fEventNr = 0; | |
316 | return; | |
317 | } | |
318 | if (firstEvent+nEvents>nall) { | |
319 | fEventNr = nall-firstEvent; | |
320 | cerr<<"restricted number of events availaible"<<endl; | |
321 | } | |
322 | AliMagF * magf = gAlice->Field(); | |
323 | AliTracker::SetFieldMap(magf,0); | |
d390cc7e | 324 | TGeoManager::Import("geometry.root"); |
325 | ||
06d06fbb | 326 | |
327 | } | |
cd875161 | 328 | //////////////////////////////////////////////////////////////////////// |
329 | AliRecInfoMaker::AliRecInfoMaker(const AliRecInfoMaker& /*info*/): | |
330 | ||
331 | fEventNr(0), //! current event number | |
332 | fNEvents(0), //! number of events to process | |
333 | fFirstEventNr(0), //! first event to process | |
334 | fFileCmp(0), //! output file with cmp tracks | |
335 | fTreeCmp(0), //! output tree with cmp tracks | |
336 | fTreeCmpKinks(0), //! output tree with cmp Kinks | |
337 | fTreeCmpV0(0), //! output tree with cmp V0 | |
338 | // | |
339 | fFileGenTracks(0), //! input files with generated tracks | |
340 | fTreeGenTracks(0), //! tree with generated tracks | |
341 | fTreeGenKinks(0), // tree with gen kinks | |
342 | fTreeGenV0(0), // tree with gen V0 | |
343 | // | |
344 | fLoader(0), //! pointer to the run loader | |
345 | // | |
346 | fIndexRecTracks(0), //! index of particle label in the TreeT_ESD | |
347 | fFakeRecTracks(0), //! number of fake tracks | |
348 | fMultiRecTracks(0), //! number of multiple reconstructions | |
349 | // | |
350 | fIndexRecKinks(0), //! index of particle label in treeesd | |
351 | fMultiRecKinks(0), //! number of multiple reconstructions | |
352 | fSignedKinks(0), //! indicator that kink was not fake | |
353 | // | |
354 | fIndexRecV0(0), //! index of particle label in treeesd | |
355 | fMultiRecV0(0), //! number of multiple reconstructions | |
356 | fSignedV0(0), //! indicator that kink was not fake | |
357 | // | |
358 | fRecArray(0), // container with rec infos | |
359 | fEvent(0), //!event | |
360 | fESDfriend(0), //!event friend | |
361 | // | |
362 | fParamTPC(0), //! AliTPCParam | |
363 | fNParticles(0), //! number of particles in the input tree genTracks | |
364 | fDebug(0), //! debug flag | |
365 | fNextTreeGenEntryToRead(0), //! last entry already read from genTracks tree | |
366 | fNextKinkToRead(0), //! last entry already read from genKinks tree | |
367 | fNextV0ToRead(0), //! last entry already read from genV0 tree | |
368 | // | |
369 | fMCInfo(0), //! MC information writen per particle | |
370 | fGenKinkInfo(0), //! MC information writen per Kink | |
371 | fGenV0Info(0), //! MC information writen per Kink | |
372 | fRecInfo(0), //! Rec. information writen per particle | |
373 | fFriend(0), //! friend track | |
374 | fRecKinkInfo(0), //! reconstructed kink info | |
375 | fRecV0Info(0) //! reconstructed kink info | |
376 | { | |
377 | // | |
378 | // Dummy copu constructor | |
379 | // | |
380 | } | |
381 | ||
382 | ||
06d06fbb | 383 | |
384 | ||
385 | //////////////////////////////////////////////////////////////////////// | |
386 | AliRecInfoMaker::~AliRecInfoMaker() | |
387 | { | |
9f3282ed | 388 | // |
389 | // Destructor | |
390 | // | |
06d06fbb | 391 | if (fLoader) { |
392 | delete fLoader; | |
393 | } | |
394 | } | |
395 | ||
396 | ////////////////////////////////////////////////////////////// | |
397 | Int_t AliRecInfoMaker::SetIO() | |
398 | { | |
399 | // | |
9f3282ed | 400 | // SetIO - Create the input trees |
401 | // | |
06d06fbb | 402 | CreateTreeCmp(); |
403 | if (!fTreeCmp) return 1; | |
404 | fParamTPC = GetTPCParam(); | |
405 | // | |
406 | if (!ConnectGenTree()) { | |
407 | cerr<<"Cannot connect tree with generated tracks"<<endl; | |
408 | return 1; | |
409 | } | |
410 | return 0; | |
411 | } | |
412 | ||
413 | ////////////////////////////////////////////////////////////// | |
414 | ||
415 | Int_t AliRecInfoMaker::SetIO(Int_t eventNr) | |
416 | { | |
417 | // | |
418 | // | |
419 | // SET INPUT | |
420 | // | |
421 | TFile f("AliESDs.root"); | |
422 | // | |
423 | ||
424 | TTree* tree = (TTree*) f.Get("esdTree"); | |
1a167b98 | 425 | tree->SetBranchStatus("*",1); |
426 | fEvent = new AliESDEvent; | |
427 | ||
428 | if (tree->GetBranch("ESD")){ | |
429 | // tree->SetBranchAddress("ESD", &fEvent); | |
430 | // tree->SetBranchAddress("ESDfriend.",&fESDfriend); | |
431 | // tree->GetEntry(eventNr); | |
432 | // fEvent->SetESDfriend(fESDfriend); | |
433 | }else{ | |
434 | fEvent->ReadFromTree(tree); | |
435 | fESDfriend = (AliESDfriend*)fEvent->FindListObject("AliESDfriend"); | |
06d06fbb | 436 | tree->GetEntry(eventNr); |
1a167b98 | 437 | fEvent->SetESDfriend(fESDfriend); |
438 | } | |
06d06fbb | 439 | |
1a167b98 | 440 | |
06d06fbb | 441 | |
442 | if (!fEvent) return 1; | |
443 | ||
444 | return 0; | |
445 | } | |
446 | ||
447 | ||
448 | ||
449 | //////////////////////////////////////////////////////////////////////// | |
450 | void AliRecInfoMaker::Reset() | |
451 | { | |
1a167b98 | 452 | // |
453 | // Reset the class | |
454 | // | |
06d06fbb | 455 | fEventNr = 0; |
456 | fNEvents = 0; | |
457 | fTreeCmp = 0; | |
458 | fTreeCmpKinks =0; | |
459 | fTreeCmpV0 =0; | |
460 | // fFnCmp = "cmpTracks.root"; | |
461 | fFileGenTracks = 0; | |
462 | fDebug = 0; | |
463 | // | |
464 | fParamTPC = 0; | |
465 | fEvent =0; | |
466 | } | |
467 | ||
468 | //////////////////////////////////////////////////////////////////////// | |
469 | Int_t AliRecInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr) | |
470 | { | |
1a167b98 | 471 | // |
472 | // Exec comparison for subrange of events | |
473 | // | |
06d06fbb | 474 | fNEvents = nEvents; |
475 | fFirstEventNr = firstEventNr; | |
476 | return Exec(); | |
477 | } | |
478 | ||
479 | //////////////////////////////////////////////////////////////////////// | |
480 | Int_t AliRecInfoMaker::Exec() | |
481 | { | |
1a167b98 | 482 | // |
483 | // Exec comparison | |
484 | // | |
06d06fbb | 485 | TStopwatch timer; |
486 | timer.Start(); | |
487 | ||
488 | if (SetIO()==1) | |
489 | return 1; | |
490 | ||
491 | fNextTreeGenEntryToRead = 0; | |
492 | fNextKinkToRead = 0; | |
493 | fNextV0ToRead =0; | |
494 | cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl; | |
495 | for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents; | |
496 | eventNr++) { | |
497 | fEventNr = eventNr; | |
498 | SetIO(fEventNr); | |
499 | fNParticles = gAlice->GetEvent(fEventNr); | |
500 | ||
501 | fIndexRecTracks = new Short_t[fNParticles*20]; //write at maximum 4 tracks corresponding to particle | |
502 | fIndexRecKinks = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle | |
503 | fIndexRecV0 = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle | |
504 | ||
505 | fFakeRecTracks = new Short_t[fNParticles]; | |
506 | fMultiRecTracks = new Short_t[fNParticles]; | |
507 | fMultiRecKinks = new Short_t[fNParticles]; | |
508 | fMultiRecV0 = new Short_t[fNParticles]; | |
509 | ||
510 | for (Int_t i = 0; i<fNParticles; i++) { | |
511 | for (Int_t j=0;j<20;j++){ | |
512 | fIndexRecTracks[20*i+j] = -1; | |
513 | fIndexRecKinks[20*i+j] = -1; | |
514 | fIndexRecV0[20*i+j] = -1; | |
515 | } | |
516 | fFakeRecTracks[i] = 0; | |
517 | fMultiRecTracks[i] = 0; | |
518 | fMultiRecKinks[i] = 0; | |
519 | fMultiRecV0[i] = 0; | |
520 | } | |
521 | ||
522 | cout<<"Start to process event "<<fEventNr<<endl; | |
523 | cout<<"\tfNParticles = "<<fNParticles<<endl; | |
524 | if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl; | |
525 | if (TreeTLoop()>0) return 1; | |
526 | ||
527 | if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl; | |
528 | if (TreeGenLoop(eventNr)>0) return 1; | |
d390cc7e | 529 | //BuildKinkInfo0(eventNr); |
a1e6aa99 | 530 | BuildV0Info(eventNr); // no V0 info for a moment |
06d06fbb | 531 | fRecArray->Delete(); |
532 | ||
533 | if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl; | |
534 | ||
535 | delete [] fIndexRecTracks; | |
536 | delete [] fIndexRecKinks; | |
537 | delete [] fIndexRecV0; | |
538 | delete [] fFakeRecTracks; | |
539 | delete [] fMultiRecTracks; | |
540 | delete [] fMultiRecKinks; | |
541 | delete [] fMultiRecV0; | |
542 | } | |
543 | ||
544 | CloseOutputFile(); | |
545 | ||
546 | cerr<<"Exec finished"<<endl; | |
547 | timer.Stop(); | |
548 | timer.Print(); | |
549 | return 0; | |
550 | ||
551 | } | |
552 | //////////////////////////////////////////////////////////////////////// | |
553 | Bool_t AliRecInfoMaker::ConnectGenTree() | |
554 | { | |
555 | // | |
556 | // connect all branches from the genTracksTree | |
557 | // use the same variables as for the new cmp tree, it may work | |
558 | // | |
559 | fFileGenTracks = TFile::Open(fFnGenTracks,"READ"); | |
560 | if (!fFileGenTracks) { | |
561 | cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl; | |
562 | return kFALSE; | |
563 | } | |
564 | fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree"); | |
565 | if (!fTreeGenTracks) { | |
566 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
567 | <<fFnGenTracks<<endl; | |
568 | return kFALSE; | |
569 | } | |
570 | // | |
571 | fMCInfo = new AliMCInfo; | |
572 | fTreeGenTracks->SetBranchAddress("MC",&fMCInfo); | |
573 | // | |
574 | // | |
575 | fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree"); | |
576 | if (!fTreeGenKinks) { | |
577 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
578 | <<fFnGenTracks<<endl; | |
579 | //return kFALSE; | |
580 | } | |
581 | else{ | |
582 | fGenKinkInfo = new AliGenKinkInfo; | |
583 | fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo); | |
584 | } | |
585 | ||
586 | fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree"); | |
587 | if (!fTreeGenV0) { | |
588 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
589 | <<fFnGenTracks<<endl; | |
590 | //return kFALSE; | |
591 | } | |
592 | else{ | |
593 | fGenV0Info = new AliGenV0Info; | |
594 | fTreeGenV0->SetBranchAddress("MC",&fGenV0Info); | |
595 | } | |
596 | // | |
597 | if (fDebug > 1) { | |
598 | cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl; | |
599 | } | |
600 | return kTRUE; | |
601 | } | |
602 | ||
603 | ||
604 | //////////////////////////////////////////////////////////////////////// | |
605 | void AliRecInfoMaker::CreateTreeCmp() | |
606 | { | |
1a167b98 | 607 | // |
608 | // Create file and tree with comparison information | |
609 | // | |
06d06fbb | 610 | fFileCmp = TFile::Open(fFnCmp,"RECREATE"); |
611 | if (!fFileCmp) { | |
612 | cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl; | |
613 | return; | |
614 | } | |
615 | // | |
616 | // | |
617 | fTreeCmp = new TTree("ESDcmpTracks","ESDcmpTracks"); | |
618 | fMCInfo = new AliMCInfo; | |
619 | fRecInfo = new AliESDRecInfo; | |
620 | AliESDtrack * esdTrack = new AliESDtrack; | |
621 | // AliITStrackMI * itsTrack = new AliITStrackMI; | |
622 | fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000); | |
623 | fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000); | |
624 | // fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack); | |
625 | delete esdTrack; | |
626 | // | |
627 | // | |
628 | fTreeCmpKinks = new TTree("ESDcmpKinks","ESDcmpKinks"); | |
629 | fGenKinkInfo = new AliGenKinkInfo; | |
630 | fRecKinkInfo = new AliESDRecKinkInfo; | |
631 | fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000); | |
632 | fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000); | |
633 | // | |
634 | // | |
635 | fTreeCmpV0 = new TTree("ESDcmpV0","ESDcmpV0"); | |
636 | fGenV0Info = new AliGenV0Info; | |
637 | fRecV0Info = new AliESDRecV0Info; | |
638 | fTreeCmpV0->Branch("MC.","AliGenV0Info", &fGenV0Info,256000); | |
639 | fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000); | |
640 | // | |
641 | fTreeCmp->AutoSave(); | |
642 | fTreeCmpKinks->AutoSave(); | |
643 | fTreeCmpV0->AutoSave(); | |
644 | } | |
1a167b98 | 645 | |
06d06fbb | 646 | //////////////////////////////////////////////////////////////////////// |
647 | void AliRecInfoMaker::CloseOutputFile() | |
648 | { | |
1a167b98 | 649 | // |
650 | // Close output file | |
651 | // | |
652 | ||
06d06fbb | 653 | if (!fFileCmp) { |
654 | cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl; | |
655 | return; | |
656 | } | |
657 | fFileCmp->cd(); | |
658 | fTreeCmp->Write(); | |
659 | delete fTreeCmp; | |
660 | ||
661 | fFileCmp->Close(); | |
662 | delete fFileCmp; | |
663 | return; | |
664 | } | |
665 | //////////////////////////////////////////////////////////////////////// | |
666 | ||
667 | TVector3 AliRecInfoMaker::TR2Local(AliTrackReference *trackRef, | |
668 | AliTPCParam *paramTPC) { | |
669 | ||
1a167b98 | 670 | // |
671 | // Transform position to the local coord frame | |
672 | // | |
673 | ||
06d06fbb | 674 | Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; |
675 | Int_t index[4]; | |
676 | paramTPC->Transform0to1(x,index); | |
677 | paramTPC->Transform1to2Ideal(x,index); | |
678 | return TVector3(x); | |
679 | } | |
680 | //////////////////////////////////////////////////////////////////////// | |
681 | ||
682 | Int_t AliRecInfoMaker::TreeTLoop() | |
683 | { | |
684 | // | |
685 | // loop over all ESD reconstructed tracks and store info in memory | |
686 | // | |
687 | // + loop over all reconstructed kinks | |
688 | TStopwatch timer; | |
689 | timer.Start(); | |
690 | // | |
691 | Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks(); | |
692 | Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks(); | |
693 | Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s(); | |
694 | fSignedKinks = new Short_t[nKinks]; | |
695 | fSignedV0 = new Short_t[nV0MIs]; | |
696 | // | |
697 | // load kinks to the memory | |
698 | for (Int_t i=0; i<nKinks;i++){ | |
9f3282ed | 699 | // AliESDkink * kink = |
700 | fEvent->GetKink(i); | |
06d06fbb | 701 | fSignedKinks[i]=0; |
702 | } | |
703 | // | |
704 | for (Int_t i=0; i<nV0MIs;i++){ | |
9f3282ed | 705 | //AliV0 * v0MI = |
706 | (AliV0*)fEvent->GetV0(i); | |
06d06fbb | 707 | fSignedV0[i]=0; |
708 | } | |
709 | ||
710 | // | |
711 | // | |
712 | AliESDtrack * track=0; | |
713 | for (Int_t iEntry=0; iEntry<nEntries;iEntry++){ | |
06d06fbb | 714 | track = (AliESDtrack*)fEvent->GetTrack(iEntry); |
715 | // | |
716 | Int_t label = track->GetLabel(); | |
717 | Int_t absLabel = abs(label); | |
718 | if (absLabel < fNParticles) { | |
719 | // fIndexRecTracks[absLabel] = iEntry; | |
720 | if (label < 0) fFakeRecTracks[absLabel]++; | |
721 | if (fMultiRecTracks[absLabel]>0){ | |
722 | if (fMultiRecTracks[absLabel]<20) | |
723 | fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] = iEntry; | |
724 | } | |
725 | else | |
726 | fIndexRecTracks[absLabel*20] = iEntry; | |
727 | fMultiRecTracks[absLabel]++; | |
728 | } | |
729 | } | |
730 | // sort reconstructed kinks | |
731 | // | |
732 | AliESDkink * kink=0; | |
733 | for (Int_t iEntry=0; iEntry<nKinks;iEntry++){ | |
734 | kink = (AliESDkink*)fEvent->GetKink(iEntry); | |
735 | if (!kink) continue; | |
736 | // | |
737 | Int_t label0 = TMath::Abs(kink->GetLabel(0)); | |
738 | Int_t label1 = TMath::Abs(kink->GetLabel(1)); | |
739 | Int_t absLabel = TMath::Min(label0,label1); | |
740 | if (absLabel < fNParticles) { | |
741 | if (fMultiRecKinks[absLabel]>0){ | |
742 | if (fMultiRecKinks[absLabel]<20) | |
743 | fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] = iEntry; | |
744 | } | |
745 | else | |
746 | fIndexRecKinks[absLabel*20] = iEntry; | |
747 | fMultiRecKinks[absLabel]++; | |
748 | } | |
749 | } | |
750 | // --sort reconstructed V0 | |
751 | // | |
f16481f0 | 752 | AliV0 * v0MI=0; |
753 | for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){ | |
754 | v0MI = (AliV0*)fEvent->GetV0(iEntry); | |
755 | if (!v0MI) continue; | |
756 | // | |
757 | // | |
758 | // | |
759 | //Int_t label0 = TMath::Abs(v0MI->GetLabel(0)); | |
760 | //Int_t label1 = TMath::Abs(v0MI->GetLabel(1)); | |
761 | AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex())); | |
762 | AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex())); | |
763 | Int_t labels[2]={-1,-1}; | |
8880c8c1 | 764 | labels[0] = (trackn==0) ? -1 : TMath::Abs(trackn->GetLabel()); |
765 | labels[1] = (trackp==0) ? -1 : TMath::Abs(trackp->GetLabel()); | |
f16481f0 | 766 | // |
767 | for (Int_t i=0;i<2;i++){ | |
8880c8c1 | 768 | Int_t absLabel = TMath::Abs(labels[i]); |
f16481f0 | 769 | if (absLabel < fNParticles) { |
770 | if (fMultiRecV0[absLabel]>0){ | |
771 | if (fMultiRecV0[absLabel]<20) | |
772 | fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] = iEntry; | |
773 | } | |
774 | else | |
775 | fIndexRecV0[absLabel*20] = iEntry; | |
776 | fMultiRecV0[absLabel]++; | |
777 | } | |
778 | } | |
779 | } | |
06d06fbb | 780 | |
781 | ||
782 | printf("Time spended in TreeTLoop\n"); | |
783 | timer.Print(); | |
784 | ||
785 | if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl; | |
786 | return 0; | |
787 | } | |
788 | ||
789 | //////////////////////////////////////////////////////////////////////// | |
790 | Int_t AliRecInfoMaker::TreeGenLoop(Int_t eventNr) | |
791 | { | |
792 | // | |
793 | // loop over all entries for a given event, find corresponding | |
794 | // rec. track and store in the fTreeCmp | |
795 | // | |
796 | TStopwatch timer; | |
797 | timer.Start(); | |
798 | Int_t entry = fNextTreeGenEntryToRead; | |
799 | Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast(); | |
800 | cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" " | |
801 | <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl; | |
802 | TBranch * branch = fTreeCmp->GetBranch("RC"); | |
1c91c693 | 803 | // TBranch * branchF = fTreeCmp->GetBranch("F"); |
06d06fbb | 804 | |
805 | branch->SetAddress(&fRecInfo); // set all pointers | |
806 | fRecArray = new TObjArray(fNParticles); | |
807 | AliESDtrack dummytrack; // | |
808 | AliESDfriendTrack dummytrackF; // | |
809 | ||
810 | while (entry < nParticlesTR) { | |
811 | fTreeGenTracks->GetEntry(entry); | |
812 | entry++; | |
813 | if (eventNr < fMCInfo->fEventNr) continue; | |
d390cc7e | 814 | if (eventNr > fMCInfo->fEventNr) continue; |
815 | if (fMCInfo->GetCharge()==0) continue; | |
06d06fbb | 816 | // |
817 | fNextTreeGenEntryToRead = entry-1; | |
818 | if (fDebug > 2 && fMCInfo->fLabel < 10) { | |
819 | cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl; | |
820 | } | |
821 | // if (fMCInfo->fNTPCRef<1) continue; // not TPCref | |
822 | // | |
823 | fRecInfo->Reset(); | |
824 | AliESDtrack * track=0; | |
825 | fRecInfo->fReconstructed =0; | |
826 | TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC); | |
827 | local.GetXYZ(fRecInfo->fTRLocalCoord); | |
828 | // | |
829 | if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) { | |
06d06fbb | 830 | track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]); |
831 | // | |
832 | // | |
833 | // find nearest track if multifound | |
834 | //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge); | |
835 | // | |
836 | Int_t status = 0; | |
837 | if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++; | |
838 | if ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++; | |
839 | if ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++; | |
840 | ||
841 | // | |
842 | if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){ | |
843 | // | |
844 | Double_t p[3]; | |
845 | track->GetInnerPxPyPz(p); | |
846 | Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; | |
847 | // | |
848 | for (Int_t i=1;i<20;i++){ | |
849 | if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){ | |
850 | AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]); | |
851 | if (!track2) continue; | |
852 | //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; // | |
853 | //if (sign2<0) continue; | |
854 | track2->GetInnerPxPyPz(p); | |
855 | Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; | |
856 | /* | |
857 | if (sign<0){ | |
858 | sign = sign2; | |
859 | track = track2; | |
860 | maxp = mom; | |
861 | continue; | |
862 | } | |
863 | */ | |
864 | // | |
865 | Int_t status2 = 0; | |
866 | if ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++; | |
867 | if ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++; | |
868 | if ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++; | |
869 | if (status2<status) continue; | |
870 | // | |
871 | if (mom<maxp) continue; | |
872 | maxp = mom; | |
873 | track = track2; | |
874 | // | |
875 | } | |
876 | } | |
877 | } | |
878 | // | |
879 | if (track) { | |
880 | fRecInfo->SetESDtrack(track); | |
881 | }else{ | |
882 | fRecInfo->SetESDtrack(&dummytrack); | |
883 | } | |
884 | // | |
885 | ||
886 | fRecInfo->fReconstructed = 1; | |
887 | fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel]; | |
888 | fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel]; | |
889 | // | |
1a167b98 | 890 | fRecInfo->Update(fMCInfo,fParamTPC,kTRUE); |
06d06fbb | 891 | } |
892 | else{ | |
893 | fRecInfo->SetESDtrack(&dummytrack); | |
1a167b98 | 894 | fRecInfo->Update(fMCInfo,fParamTPC,kFALSE); |
06d06fbb | 895 | } |
896 | fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel); | |
897 | fTreeCmp->Fill(); | |
898 | } | |
899 | fTreeCmp->AutoSave(); | |
06d06fbb | 900 | printf("Time spended in TreeGenLoop\n"); |
901 | timer.Print(); | |
902 | if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl; | |
903 | ||
904 | return 0; | |
905 | } | |
906 | ||
907 | ||
908 | ||
909 | //////////////////////////////////////////////////////////////////////// | |
910 | //////////////////////////////////////////////////////////////////////// | |
911 | //////////////////////////////////////////////////////////////////////// | |
912 | Int_t AliRecInfoMaker::BuildKinkInfo0(Int_t eventNr) | |
913 | { | |
914 | // | |
915 | // loop over all entries for a given event, find corresponding | |
916 | // rec. track and store in the fTreeCmp | |
917 | // | |
918 | TStopwatch timer; | |
919 | timer.Start(); | |
920 | Int_t entry = fNextKinkToRead; | |
921 | Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast(); | |
922 | cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" " | |
923 | <<nParticlesTR<<" "<<fNextKinkToRead<<endl; | |
924 | // | |
925 | TBranch * branch = fTreeCmpKinks->GetBranch("RC."); | |
926 | branch->SetAddress(&fRecKinkInfo); // set all pointers | |
927 | ||
928 | // | |
929 | while (entry < nParticlesTR) { | |
930 | fTreeGenKinks->GetEntry(entry); | |
931 | entry++; | |
932 | if (eventNr < fGenKinkInfo->GetMinus().fEventNr) continue; | |
933 | if (eventNr > fGenKinkInfo->GetMinus().fEventNr) continue;; | |
934 | // | |
935 | fNextKinkToRead = entry-1; | |
936 | // | |
937 | // | |
938 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetMinus().fLabel); | |
939 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetPlus().fLabel); | |
940 | fRecKinkInfo->fT1 = (*fRecInfo1); | |
941 | fRecKinkInfo->fT2 = (*fRecInfo2); | |
942 | fRecKinkInfo->fStatus =0; | |
943 | if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1; | |
944 | if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2; | |
945 | if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1; | |
946 | ||
947 | if (fRecKinkInfo->fStatus==3){ | |
948 | fRecKinkInfo->Update(); | |
949 | } | |
950 | Int_t label = TMath::Min(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel); | |
951 | Int_t label2 = TMath::Max(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel); | |
952 | ||
953 | AliESDkink *kink=0; | |
954 | fRecKinkInfo->fRecStatus =0; | |
955 | fRecKinkInfo->fMultiple = fMultiRecKinks[label]; | |
956 | fRecKinkInfo->fKinkMultiple=0; | |
957 | // | |
958 | if (fMultiRecKinks[label]>0){ | |
959 | ||
960 | // for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){ | |
961 | for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){ | |
962 | Int_t index = fIndexRecKinks[label*20+j]; | |
963 | //AliESDkink *kink2 = (AliESDkink*)fKinks->At(index); | |
964 | AliESDkink *kink2 = (AliESDkink*)fEvent->GetKink(index); | |
965 | if (TMath::Abs(kink2->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(1))==label2) { | |
966 | fRecKinkInfo->fKinkMultiple++; | |
967 | fSignedKinks[index]=1; | |
968 | Int_t c0=0; | |
969 | if (kink){ | |
970 | // if (kink->fTRDOn) c0++; | |
971 | //if (kink->fITSOn) c0++; | |
972 | if (kink->GetStatus(2)>0) c0++; | |
973 | if (kink->GetStatus(0)>0) c0++; | |
974 | } | |
975 | Int_t c2=0; | |
976 | // if (kink2->fTRDOn) c2++; | |
977 | //if (kink2->fITSOn) c2++; | |
978 | if (kink2->GetStatus(2)>0) c2++; | |
979 | if (kink2->GetStatus(0)>0) c2++; | |
980 | ||
981 | if (c2<c0) continue; | |
982 | kink =kink2; | |
983 | } | |
984 | if (TMath::Abs(kink2->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(0))==label2) { | |
985 | fRecKinkInfo->fKinkMultiple++; | |
986 | fSignedKinks[index]=1; | |
987 | Int_t c0=0; | |
988 | if (kink){ | |
989 | //if (kink->fTRDOn) c0++; | |
990 | //if (kink->fITSOn) c0++; | |
991 | if (kink->GetStatus(2)>0) c0++; | |
992 | if (kink->GetStatus(0)>0) c0++; | |
993 | ||
994 | } | |
995 | Int_t c2=0; | |
996 | // if (kink2->fTRDOn) c2++; | |
997 | //if (kink2->fITSOn) c2++; | |
998 | if (kink2->GetStatus(2)>0) c2++; | |
999 | if (kink2->GetStatus(0)>0) c2++; | |
1000 | ||
1001 | if (c2<c0) continue; | |
1002 | kink =kink2; | |
1003 | } | |
1004 | } | |
1005 | } | |
1006 | if (kink){ | |
1007 | fRecKinkInfo->fKink = *kink; | |
1008 | fRecKinkInfo->fRecStatus=1; | |
1009 | } | |
1010 | fTreeCmpKinks->Fill(); | |
1011 | } | |
1012 | // Int_t nkinks = fKinks->GetEntriesFast(); | |
1013 | Int_t nkinks = fEvent->GetNumberOfKinks(); | |
1014 | for (Int_t i=0;i<nkinks;i++){ | |
1015 | if (fSignedKinks[i]==0){ | |
1016 | // AliESDkink *kink = (AliESDkink*)fKinks->At(i); | |
1017 | AliESDkink *kink = (AliESDkink*)fEvent->GetKink(i); | |
1018 | if (!kink) continue; | |
1019 | // | |
1020 | fRecKinkInfo->fKink = *kink; | |
1021 | fRecKinkInfo->fRecStatus =-2; | |
1022 | // | |
1023 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(0))); | |
1024 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(1))); | |
1025 | if (fRecInfo1 && fRecInfo2){ | |
1026 | fRecKinkInfo->fT1 = (*fRecInfo1); | |
1027 | fRecKinkInfo->fT2 = (*fRecInfo2); | |
1028 | fRecKinkInfo->fRecStatus =-1; | |
1029 | } | |
1030 | fTreeCmpKinks->Fill(); | |
1031 | } | |
1032 | } | |
1033 | ||
1034 | ||
1035 | fTreeCmpKinks->AutoSave(); | |
1036 | printf("Time spended in BuilKinkInfo Loop\n"); | |
1037 | timer.Print(); | |
1038 | if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl; | |
1039 | return 0; | |
1040 | } | |
1041 | ||
1042 | ||
1043 | ||
1044 | //////////////////////////////////////////////////////////////////////// | |
1045 | //////////////////////////////////////////////////////////////////////// | |
1046 | //////////////////////////////////////////////////////////////////////// | |
1047 | ||
1048 | ||
1049 | ||
1050 | Int_t AliRecInfoMaker::BuildV0Info(Int_t eventNr) | |
1051 | { | |
1052 | // | |
1053 | // loop over all entries for a given event, find corresponding | |
1054 | // rec. track and store in the fTreeCmp | |
1055 | // | |
f16481f0 | 1056 | static TDatabasePDG pdgtable; |
1057 | ||
06d06fbb | 1058 | TStopwatch timer; |
1059 | timer.Start(); | |
1060 | Int_t entry = fNextV0ToRead; | |
1061 | Double_t nParticlesTR = fTreeGenV0->GetEntriesFast(); | |
1062 | cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" " | |
1063 | <<nParticlesTR<<" "<<fNextV0ToRead<<endl; | |
1064 | // | |
1065 | TBranch * branch = fTreeCmpV0->GetBranch("RC."); | |
1066 | branch->SetAddress(&fRecV0Info); // set all pointers | |
1067 | const AliESDVertex * esdvertex = fEvent->GetVertex(); | |
1068 | Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()}; | |
1069 | ||
1070 | // | |
1071 | while (entry < nParticlesTR) { | |
1072 | fTreeGenV0->GetEntry(entry); | |
1073 | entry++; | |
34737198 | 1074 | fRecV0Info->Reset(); //reset all variables |
06d06fbb | 1075 | if (eventNr < fGenV0Info->GetMinus().fEventNr) continue; |
1076 | if (eventNr > fGenV0Info->GetMinus().fEventNr) continue;; | |
1077 | // | |
1078 | fNextV0ToRead = entry-1; | |
1079 | // | |
1080 | // | |
1081 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetMinus().fLabel); | |
1082 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetPlus().fLabel); | |
1083 | if (fGenV0Info->GetMinus().fCharge*fGenV0Info->GetPlus().fCharge>0) continue; // interactions | |
1084 | if (!fRecInfo1 || !fRecInfo2) continue; | |
1085 | fRecV0Info->fT1 = (*fRecInfo1); | |
1086 | fRecV0Info->fT2 = (*fRecInfo2); | |
1087 | fRecV0Info->fV0Status =0; | |
1088 | if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1; | |
1089 | if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2; | |
1090 | ||
1091 | if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1; | |
1092 | ||
1093 | ||
1094 | if (abs(fRecV0Info->fV0Status)==3){ | |
1095 | fRecV0Info->Update(vertex); | |
1096 | { | |
1097 | // | |
1098 | // TPC V0 Info | |
1099 | Double_t x,alpha, param[5],cov[15]; | |
1100 | if ( fRecV0Info->fT1.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){ | |
1101 | fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param); | |
1102 | fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov); | |
1103 | AliExternalTrackParam paramP(x,alpha,param,cov); | |
1104 | // | |
1105 | fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param); | |
1106 | fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov); | |
1107 | AliExternalTrackParam paramM(x,alpha,param,cov); | |
1108 | // | |
1109 | fRecV0Info->fV0tpc->SetParamN(paramM); | |
1110 | fRecV0Info->fV0tpc->SetParamP(paramP); | |
1111 | Double_t pid1[5],pid2[5]; | |
1112 | fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1); | |
1113 | fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2); | |
1114 | // | |
1115 | //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2); | |
1116 | fRecV0Info->fV0tpc->Update(vertex); | |
1117 | ||
1118 | // | |
1119 | // | |
1120 | fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param); | |
1121 | fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov); | |
1122 | alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha(); | |
1123 | new (¶mP) AliExternalTrackParam(x,alpha,param,cov); | |
1124 | // | |
1125 | fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param); | |
1126 | fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov); | |
1127 | alpha = fRecV0Info->fT2.GetESDtrack()->GetAlpha(); | |
1128 | new (¶mM) AliExternalTrackParam(x,alpha,param,cov); | |
1129 | // | |
1130 | fRecV0Info->fV0its->SetParamN(paramM); | |
1131 | fRecV0Info->fV0its->SetParamP(paramP); | |
1132 | // fRecV0Info->fV0its.UpdatePID(pid1,pid2); | |
1133 | fRecV0Info->fV0its->Update(vertex); | |
1134 | } | |
1135 | } | |
f16481f0 | 1136 | // |
1137 | // ???? | |
1138 | // | |
06d06fbb | 1139 | if (TMath::Abs(fGenV0Info->GetMinus().fPdg)==11 &&TMath::Abs(fGenV0Info->GetPlus().fPdg)==11){ |
1140 | if (fRecV0Info->fDist2>10){ | |
1141 | fRecV0Info->Update(vertex); | |
1142 | } | |
1143 | if (fRecV0Info->fDist2>10){ | |
1144 | fRecV0Info->Update(vertex); | |
1145 | } | |
1146 | } | |
1147 | } | |
1148 | // | |
1149 | // take the V0 from reconstruction | |
1150 | ||
1151 | Int_t label = TMath::Min(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel); | |
1152 | Int_t label2 = TMath::Max(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel); | |
1153 | AliV0 *v0MI=0; | |
f16481f0 | 1154 | AliV0 *v0MIOff=0; |
06d06fbb | 1155 | fRecV0Info->fRecStatus =0; |
1156 | fRecV0Info->fMultiple = fMultiRecV0[label]; | |
f16481f0 | 1157 | fRecV0Info->fV0MultipleOn=0; |
1158 | fRecV0Info->fV0MultipleOff=0; | |
06d06fbb | 1159 | // |
1160 | if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){ | |
1161 | ||
1162 | // for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){ | |
a1e6aa99 | 1163 | for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){ |
1164 | Int_t index = fIndexRecV0[label*20+j]; | |
1165 | if (index<0) continue; | |
1166 | AliV0 *v0MI2 = (AliV0*)fEvent->GetV0(index); | |
1167 | // get track labels | |
1168 | AliESDtrack * trackn = fEvent->GetTrack((v0MI2->GetNindex())); | |
1169 | AliESDtrack * trackp = fEvent->GetTrack((v0MI2->GetPindex())); | |
1170 | Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel(); | |
1171 | Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel(); | |
34737198 | 1172 | fRecV0Info->fLab[0]=TMath::Abs(vlabelp); |
1173 | fRecV0Info->fLab[1]=TMath::Abs(vlabeln); | |
a1e6aa99 | 1174 | // |
1175 | if (TMath::Abs(vlabeln)==label &&TMath::Abs(vlabelp)==label2) { | |
f16481f0 | 1176 | if (v0MI2->GetOnFlyStatus()) { |
1177 | v0MI =v0MI2; | |
1178 | fRecV0Info->fV0MultipleOn++; | |
1179 | }else { | |
1180 | v0MIOff = v0MI2; | |
1181 | fRecV0Info->fV0MultipleOff++; | |
1182 | } | |
a1e6aa99 | 1183 | fSignedV0[index]=1; |
1184 | } | |
1185 | if (TMath::Abs(vlabelp)==label &&TMath::Abs(vlabeln)==label2) { | |
f16481f0 | 1186 | if (v0MI2->GetOnFlyStatus()){ |
1187 | v0MI =v0MI2; | |
1188 | fRecV0Info->fV0MultipleOn++; | |
1189 | }else { | |
1190 | v0MIOff = v0MI2; | |
1191 | fRecV0Info->fV0MultipleOff++; | |
1192 | } | |
a1e6aa99 | 1193 | fSignedV0[index]=1; |
1194 | } | |
1195 | } | |
06d06fbb | 1196 | } |
1197 | if (v0MI){ | |
f16481f0 | 1198 | new (fRecV0Info->fV0rec) AliV0(*v0MI); |
06d06fbb | 1199 | fRecV0Info->fRecStatus=1; |
1200 | } | |
f16481f0 | 1201 | if (v0MIOff){ |
1202 | new (fRecV0Info->fV0recOff) AliV0(*v0MIOff); | |
1203 | fRecV0Info->fRecStatus=1; | |
1204 | } | |
1205 | Int_t mpdg = fGenV0Info->GetMother().GetPdgCode(); | |
1206 | Float_t mass = ( pdgtable.GetParticle(mpdg)==0) ? 0 :pdgtable.GetParticle(mpdg)->Mass(); | |
1207 | fRecV0Info->UpdateKF(*esdvertex, | |
1208 | fGenV0Info->GetPlus().GetPdg(), | |
1209 | fGenV0Info->GetMinus().GetPdg(), | |
1210 | mass); | |
06d06fbb | 1211 | fTreeCmpV0->Fill(); |
1212 | } | |
1213 | // | |
1214 | // write fake v0s | |
a1e6aa99 | 1215 | // |
f16481f0 | 1216 | Int_t nV0MIs = fEvent->GetNumberOfV0s(); |
1217 | for (Int_t i=0;i<nV0MIs;i++){ | |
1218 | if (fSignedV0[i]==0){ | |
1219 | AliV0 *v0MI = (AliV0*)fEvent->GetV0(i); | |
1220 | if (!v0MI) continue; | |
34737198 | 1221 | fRecV0Info->Reset(); //reset all variables |
f16481f0 | 1222 | // |
1223 | new (fRecV0Info->fV0rec) AliV0(*v0MI); | |
1224 | fRecV0Info->fV0Status =-10; | |
1225 | fRecV0Info->fRecStatus =-2; | |
1226 | // | |
1227 | AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex())); | |
1228 | AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex())); | |
1229 | Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel(); | |
1230 | Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel(); | |
34737198 | 1231 | fRecV0Info->fLab[0]=TMath::Abs(vlabelp); |
1232 | fRecV0Info->fLab[1]=TMath::Abs(vlabeln); | |
1233 | if (TMath::Abs(fRecV0Info->fLab[0] - fRecV0Info->fLab[1])<2) continue; | |
f16481f0 | 1234 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabeln)); |
1235 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabelp)); | |
1236 | if (fRecInfo1 && fRecInfo2){ | |
1237 | fRecV0Info->fT1 = (*fRecInfo1); | |
1238 | fRecV0Info->fT2 = (*fRecInfo2); | |
1239 | fRecV0Info->fRecStatus =-1; | |
1240 | } | |
1241 | fRecV0Info->Update(vertex); | |
1242 | fRecV0Info->UpdateKF(*esdvertex,211,211,0.49767); | |
1243 | fTreeCmpV0->Fill(); | |
1244 | } | |
1245 | } | |
06d06fbb | 1246 | |
1247 | ||
1248 | ||
1249 | fTreeCmpV0->AutoSave(); | |
1250 | printf("Time spended in BuilV0Info Loop\n"); | |
1251 | timer.Print(); | |
1252 | if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl; | |
1253 | return 0; | |
1254 | } | |
1255 | //////////////////////////////////////////////////////////////////////// | |
1256 | //////////////////////////////////////////////////////////////////////// | |
1257 | ||
1258 |