]>
Commit | Line | Data |
---|---|---|
ae7d73d2 | 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 | .L $ALICE_ROOT/STEER/AliGenInfo.C+ | |
32 | //be sure you created genTracks file before | |
33 | .L $ALICE_ROOT/STEER/AliESDComparisonMI.C+ | |
34 | ||
35 | // | |
36 | ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,1,0); | |
37 | t2->Exec(); | |
38 | ||
39 | // | |
40 | //some cuts definition | |
51ad6848 | 41 | TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01") |
ae7d73d2 | 42 | //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5") |
51ad6848 | 43 | //TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9"); |
f007cb5f | 44 | TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5"); |
ae7d73d2 | 45 | TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5"); |
46 | ||
47 | ||
f007cb5f | 48 | |
49 | ||
50 | ||
ae7d73d2 | 51 | TCut crec("crec","fReconstructed==1"); |
51ad6848 | 52 | TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25"); |
f007cb5f | 53 | TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1"); |
54 | ||
ae7d73d2 | 55 | TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1"); |
56 | TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50"); | |
57 | TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)"); | |
58 | TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.") | |
59 | ||
60 | AliESDComparisonDraw comp; | |
61 | comp.SetIO(); | |
51ad6848 | 62 | TFile f("genHits.root"); |
63 | TTree * treel = (TTree*)f.Get("HitLines"); | |
64 | if (treel) comp->fTree->AddFriend(treel,"L"); | |
ae7d73d2 | 65 | |
66 | // | |
67 | //example | |
68 | comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); | |
51ad6848 | 69 | comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); |
70 | comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); | |
71 | comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()"); | |
72 | comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); | |
51ad6848 | 73 | comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); |
f007cb5f | 74 | comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); |
75 | ||
ae7d73d2 | 76 | |
77 | TH1F his("his","his",100,0,20); | |
78 | TH1F hpools("hpools","hpools",100,-7,7); | |
79 | TH1F hfake("hfake","hfake",1000,0,150); | |
80 | TProfile profp0("profp0","profp0",20,-0.4,0.9) | |
81 | ||
f007cb5f | 82 | comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06) |
ae7d73d2 | 83 | comp.fRes->Draw(); |
84 | comp.fMean->Draw(); | |
51ad6848 | 85 | |
f007cb5f | 86 | comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06) |
51ad6848 | 87 | comp.fRes->Draw(); |
88 | ||
89 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5) | |
90 | comp.fRes->Draw(); | |
91 | ||
92 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5) | |
93 | comp.fRes->Draw(); | |
94 | comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5) | |
ae7d73d2 | 95 | comp.fRes->Draw(); |
96 | ||
97 | comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") | |
98 | ||
99 | TH1F his("his","his",100,0,20); | |
100 | TH1F hpools("hpools","hpools",100,-7,7); | |
101 | ||
f007cb5f | 102 | TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1); |
103 | TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4); | |
104 | TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3); | |
105 | TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2); | |
106 | ||
107 | comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) | |
108 | comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) | |
109 | comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) | |
110 | comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) | |
51ad6848 | 111 | |
51ad6848 | 112 | |
f007cb5f | 113 | comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") |
114 | comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") | |
115 | comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") | |
116 | comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") | |
51ad6848 | 117 | |
f007cb5f | 118 | comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") |
119 | comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") | |
120 | comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") | |
121 | comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") | |
51ad6848 | 122 | |
123 | ||
f007cb5f | 124 | hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same"); |
ae7d73d2 | 125 | |
126 | comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8) | |
127 | ||
f007cb5f | 128 | TProfile prof("prof","prof",10,0.5,5); |
129 | ||
130 | ||
131 | ||
132 | ||
ae7d73d2 | 133 | */ |
134 | ||
135 | ||
136 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
137 | #include <stdio.h> | |
138 | #include <string.h> | |
139 | //ROOT includes | |
140 | #include "Rtypes.h" | |
141 | #include "TFile.h" | |
142 | #include "TTree.h" | |
143 | #include "TChain.h" | |
144 | #include "TCut.h" | |
145 | #include "TString.h" | |
146 | #include "TBenchmark.h" | |
147 | #include "TStopwatch.h" | |
148 | #include "TParticle.h" | |
149 | #include "TSystem.h" | |
150 | #include "TTimer.h" | |
151 | #include "TVector3.h" | |
152 | #include "TPad.h" | |
153 | #include "TCanvas.h" | |
154 | #include "TH1F.h" | |
155 | #include "TH2F.h" | |
156 | #include "TF1.h" | |
157 | #include "TText.h" | |
158 | #include "Getline.h" | |
159 | #include "TStyle.h" | |
160 | ||
161 | //ALIROOT includes | |
162 | #include "AliRun.h" | |
163 | #include "AliStack.h" | |
164 | #include "AliESDtrack.h" | |
165 | #include "AliSimDigits.h" | |
166 | #include "AliTPCParam.h" | |
167 | #include "AliTPC.h" | |
168 | #include "AliTPCLoader.h" | |
169 | #include "AliDetector.h" | |
170 | #include "AliTrackReference.h" | |
171 | #include "AliRun.h" | |
172 | #include "AliTPCParamSR.h" | |
173 | #include "AliTracker.h" | |
174 | #include "AliComplexCluster.h" | |
175 | #include "AliMagF.h" | |
176 | #include "AliESD.h" | |
177 | #include "AliESDtrack.h" | |
51ad6848 | 178 | #include "AliITStrackMI.h" |
ae7d73d2 | 179 | #include "AliTRDtrack.h" |
51ad6848 | 180 | #include "AliHelix.h" |
181 | #include "AliESDVertex.h" | |
182 | #include "AliExternalTrackParam.h" | |
183 | #include "AliESDkink.h" | |
184 | #include "AliESDV0MI.h" | |
185 | ||
ae7d73d2 | 186 | #endif |
187 | #include "AliGenInfo.h" | |
188 | #include "AliESDComparisonMI.h" | |
189 | ||
190 | ||
191 | ||
81e97e0d | 192 | TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01"); |
193 | TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5"); | |
194 | TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5"); | |
195 | TCut crec("crec","fReconstructed==1"); | |
196 | TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25"); | |
197 | TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1"); | |
198 | ||
199 | TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1"); | |
200 | TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50"); | |
201 | TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)"); | |
202 | TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6."); | |
203 | ||
204 | ||
205 | ||
206 | void MakeAliases(AliESDComparisonDraw&comp) | |
207 | { | |
208 | // | |
209 | // aliases definition | |
210 | // | |
211 | comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); | |
212 | comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); | |
213 | comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); | |
214 | comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()"); | |
215 | comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); | |
216 | comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
217 | comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); | |
218 | ||
219 | comp.fTree->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."); | |
220 | ||
221 | comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)"); | |
222 | comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)"); | |
223 | ||
224 | comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]"); | |
225 | comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum"); | |
226 | comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum"); | |
227 | comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum"); | |
228 | comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum"); | |
229 | comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum"); | |
230 | comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)"); | |
231 | } | |
232 | ||
233 | ||
234 | ||
51ad6848 | 235 | |
236 | void AliESDRecInfo::UpdatePoints(AliESDtrack*track) | |
237 | { | |
238 | // | |
239 | // | |
240 | Int_t iclusters[200]; | |
241 | Float_t density[160]; | |
242 | for (Int_t i=0;i<160;i++) density[i]=-1.; | |
243 | fTPCPoints[0]= 160; | |
244 | fTPCPoints[1] = -1; | |
245 | // | |
246 | if (fTPCPoints[0]<fTPCPoints[1]) return; | |
247 | // Int_t nclusters=track->GetTPCclusters(iclusters); | |
248 | ||
249 | Int_t ngood=0; | |
250 | Int_t undeff=0; | |
251 | Int_t nall =0; | |
252 | Int_t range=20; | |
253 | for (Int_t i=0;i<160;i++){ | |
254 | Int_t last = i-range; | |
255 | if (nall<range) nall++; | |
256 | if (last>=0){ | |
257 | if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--; | |
258 | if (iclusters[last]==-1) undeff--; | |
259 | } | |
260 | if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++; | |
261 | if (iclusters[i]==-1) undeff++; | |
262 | if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff); | |
263 | } | |
264 | Float_t maxdens=0; | |
265 | Int_t indexmax =0; | |
266 | for (Int_t i=0;i<160;i++){ | |
267 | if (density[i]<0) continue; | |
268 | if (density[i]>maxdens){ | |
269 | maxdens=density[i]; | |
270 | indexmax=i; | |
271 | } | |
272 | } | |
273 | // | |
274 | //max dens point | |
275 | fTPCPoints[3] = maxdens; | |
276 | fTPCPoints[1] = indexmax; | |
277 | // | |
278 | // last point | |
279 | for (Int_t i=indexmax;i<160;i++){ | |
280 | if (density[i]<0) continue; | |
281 | if (density[i]<maxdens/2.) { | |
282 | break; | |
283 | } | |
284 | fTPCPoints[2]=i; | |
285 | } | |
286 | // | |
287 | // first point | |
288 | for (Int_t i=indexmax;i>0;i--){ | |
289 | if (density[i]<0) continue; | |
290 | if (density[i]<maxdens/2.) { | |
291 | break; | |
292 | } | |
293 | fTPCPoints[0]=i; | |
294 | } | |
295 | // | |
f007cb5f | 296 | // Density at the last 30 padrows |
297 | // | |
298 | // | |
299 | nall = 0; | |
300 | ngood = 0; | |
301 | for (Int_t i=159;i>0;i--){ | |
302 | if (iclusters[i]==-1) continue; //dead zone | |
303 | nall++; | |
304 | if (iclusters[i]>0) ngood++; | |
305 | if (nall>20) break; | |
306 | } | |
307 | fTPCPoints[4] = Float_t(ngood)/Float_t(nall); | |
308 | // | |
51ad6848 | 309 | if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1; |
f007cb5f | 310 | // |
311 | // | |
312 | // check TRDPoints | |
313 | /* | |
314 | nclusters=track->GetTRDclusters(iclusters); | |
315 | for (Int_t i=nclusters;i>0;i--){ | |
316 | ||
317 | } | |
318 | */ | |
319 | ||
320 | ||
51ad6848 | 321 | } |
322 | ||
ae7d73d2 | 323 | // |
324 | // | |
f007cb5f | 325 | void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed) |
ae7d73d2 | 326 | { |
327 | // | |
328 | // | |
329 | //calculates derived variables | |
330 | // | |
331 | // | |
51ad6848 | 332 | UpdatePoints(&fESDTrack); |
333 | fBestTOFmatch=1000; | |
ae7d73d2 | 334 | AliTrackReference * ref = &(info->fTrackRef); |
335 | fTPCinR0[0] = info->fTrackRef.X(); | |
336 | fTPCinR0[1] = info->fTrackRef.Y(); | |
337 | fTPCinR0[2] = info->fTrackRef.Z(); | |
338 | fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]); | |
339 | fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]); | |
340 | // | |
341 | fTPCinP0[0] = ref->Px(); | |
342 | fTPCinP0[1] = ref->Py(); | |
343 | fTPCinP0[2] = ref->Pz(); | |
344 | fTPCinP0[3] = ref->Pt(); | |
345 | fTPCinP0[4] = ref->P(); | |
346 | fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P(); | |
347 | // | |
348 | // | |
349 | if (fTPCinP0[3]>0.0000001){ | |
350 | // | |
351 | fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]); | |
352 | fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]); | |
353 | } | |
354 | // | |
355 | // | |
356 | fITSinP0[0]=info->fParticle.Px(); | |
357 | fITSinP0[1]=info->fParticle.Py(); | |
358 | fITSinP0[2]=info->fParticle.Pz(); | |
359 | fITSinP0[3]=info->fParticle.Pt(); | |
360 | // | |
361 | fITSinR0[0]=info->fParticle.Vx(); | |
362 | fITSinR0[1]=info->fParticle.Vy(); | |
363 | fITSinR0[2]=info->fParticle.Vz(); | |
364 | fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]); | |
365 | fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]); | |
366 | // | |
367 | // | |
368 | if (fITSinP0[3]>0.0000001){ | |
369 | fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]); | |
370 | fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]); | |
371 | } | |
372 | // | |
51ad6848 | 373 | for (Int_t i=0;i<4;i++) fStatus[i] =0; |
374 | fReconstructed = kFALSE; | |
375 | fTPCOn = kFALSE; | |
376 | fITSOn = kFALSE; | |
377 | fTRDOn = kFALSE; | |
378 | if (reconstructed==kFALSE) return; | |
379 | ||
380 | fLabels[0] = info->fLabel; | |
381 | fLabels[1] = info->fPrimPart; | |
ae7d73d2 | 382 | fReconstructed = kTRUE; |
383 | fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE; | |
384 | fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE; | |
385 | fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE; | |
51ad6848 | 386 | // |
387 | // | |
388 | if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){ | |
389 | fStatus[1] =3; | |
390 | } | |
391 | else{ | |
392 | if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){ | |
393 | fStatus[1] =2; | |
394 | } | |
395 | else{ | |
396 | fStatus[1]=1; | |
397 | } | |
398 | } | |
f007cb5f | 399 | // |
400 | if ((fESDTrack.GetStatus()&AliESDtrack::kITSout)>0){ | |
401 | fStatus[0] =2; | |
402 | } | |
403 | else{ | |
404 | if ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0){ | |
405 | fStatus[0] =1; | |
406 | } | |
407 | else{ | |
408 | fStatus[0]=0; | |
409 | } | |
410 | } | |
411 | ||
51ad6848 | 412 | // |
413 | // | |
414 | if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){ | |
415 | fStatus[2] =2; | |
416 | } | |
417 | else{ | |
418 | if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){ | |
419 | fStatus[2] =1; | |
420 | } | |
421 | } | |
422 | if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){ | |
423 | fStatus[2] =10; | |
424 | } | |
425 | ||
426 | // | |
427 | //TOF | |
428 | // | |
429 | if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){ | |
430 | // | |
431 | // best tof match | |
432 | Double_t times[5]; | |
433 | fESDTrack.GetIntegratedTimes(times); | |
434 | for (Int_t i=0;i<5;i++){ | |
435 | if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){ | |
436 | fBestTOFmatch = fESDTrack.GetTOFsignal()-times[i]; | |
437 | } | |
438 | } | |
439 | Int_t toflabel[3]; | |
440 | fESDTrack.GetTOFLabel(toflabel); | |
441 | Bool_t toffake=kTRUE; | |
f007cb5f | 442 | Bool_t tofdaughter=kFALSE; |
51ad6848 | 443 | for (Int_t i=0;i<3;i++){ |
f007cb5f | 444 | if (toflabel[i]<0) continue; |
51ad6848 | 445 | if (toflabel[i]== TMath::Abs(fESDTrack.GetLabel())) toffake=kFALSE; |
f007cb5f | 446 | if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle |
51ad6848 | 447 | fStatus[3]=1; |
448 | } | |
f007cb5f | 449 | if (toffake) fStatus[3] =3; //total fake |
450 | if (tofdaughter) fStatus[3]=2; //fake because of decay | |
51ad6848 | 451 | }else{ |
452 | fStatus[3]=0; | |
453 | } | |
454 | ||
455 | ||
81e97e0d | 456 | if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){ |
ae7d73d2 | 457 | //TPC |
458 | fESDTrack.GetInnerXYZ(fTPCinR1); | |
459 | fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]); | |
460 | fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]); | |
461 | fESDTrack.GetInnerPxPyPz(fTPCinP1); | |
462 | fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]); | |
463 | fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]); | |
464 | // | |
465 | // | |
466 | if (fTPCinP1[3]>0.0000001){ | |
467 | fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]); | |
468 | fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]); | |
469 | } | |
470 | Double_t cov[15], param[5],x; | |
471 | fESDTrack.GetInnerExternalCovariance(cov); | |
51ad6848 | 472 | fESDTrack.GetInnerExternalParameters(x,param); |
f007cb5f | 473 | if (x<50) return ; |
ae7d73d2 | 474 | // |
475 | fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi | |
476 | fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]); | |
477 | fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z | |
478 | fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]); | |
479 | fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]); | |
480 | fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]); | |
481 | fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1])); | |
482 | fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]); | |
483 | fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]); | |
484 | Double_t sign = (param[4]>0)? 1.:-1; | |
51ad6848 | 485 | fSign =sign; |
ae7d73d2 | 486 | fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]); |
487 | } | |
488 | if (fITSOn){ | |
489 | // ITS | |
490 | Double_t param[5],x; | |
f007cb5f | 491 | fESDTrack.GetExternalParameters(x,param); |
492 | // fESDTrack.GetConstrainedExternalParameters(x,param); | |
ae7d73d2 | 493 | Double_t cov[15]; |
494 | fESDTrack.GetExternalCovariance(cov); | |
f007cb5f | 495 | //fESDTrack.GetConstrainedExternalCovariance(cov); |
ae7d73d2 | 496 | if (TMath::Abs(param[4])<0.0000000001) return; |
497 | ||
498 | fESDTrack.GetXYZ(fITSinR1); | |
499 | fESDTrack.GetPxPyPz(fITSinP1); | |
500 | fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]); | |
501 | // | |
502 | fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]); | |
503 | fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]); | |
504 | // | |
505 | // | |
506 | if (fITSinP1[3]>0.0000001){ | |
507 | fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]); | |
508 | fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]); | |
509 | } | |
510 | // | |
511 | // | |
512 | fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi | |
513 | fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]); | |
514 | fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z | |
515 | fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]); | |
516 | fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]); | |
517 | fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]); | |
518 | fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1])); | |
519 | fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]); | |
520 | fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]); | |
521 | Double_t sign = (param[4]>0) ? 1:-1; | |
51ad6848 | 522 | fSign = sign; |
ae7d73d2 | 523 | fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]); |
524 | } | |
525 | ||
526 | } | |
527 | ||
528 | ||
529 | void AliESDRecV0Info::Update(Float_t vertex[3]) | |
51ad6848 | 530 | { |
531 | ||
532 | if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){ | |
533 | Float_t distance1,distance2; | |
534 | Double_t xx[3],pp[3]; | |
535 | // | |
536 | Double_t xd[3],pd[3],signd; | |
537 | Double_t xm[3],pm[3],signm; | |
f007cb5f | 538 | // |
539 | // | |
540 | if (fT1.fITSOn&&fT2.fITSOn){ | |
541 | for (Int_t i=0;i<3;i++){ | |
542 | xd[i] = fT2.fITSinR1[i]; | |
543 | pd[i] = fT2.fITSinP1[i]; | |
544 | xm[i] = fT1.fITSinR1[i]; | |
545 | pm[i] = fT1.fITSinP1[i]; | |
546 | } | |
547 | } | |
548 | else{ | |
549 | ||
550 | for (Int_t i=0;i<3;i++){ | |
551 | xd[i] = fT2.fTPCinR1[i]; | |
552 | pd[i] = fT2.fTPCinP1[i]; | |
553 | xm[i] = fT1.fTPCinR1[i]; | |
554 | pm[i] = fT1.fTPCinP1[i]; | |
555 | } | |
51ad6848 | 556 | } |
f007cb5f | 557 | // |
558 | // | |
51ad6848 | 559 | signd = fT2.fSign<0 ? -1:1; |
560 | signm = fT1.fSign<0 ? -1:1; | |
561 | ||
562 | AliHelix dhelix1(xd,pd,signd); | |
563 | dhelix1.GetMomentum(0,pp,0); | |
564 | dhelix1.Evaluate(0,xx); | |
565 | // | |
566 | // Double_t x2[3],p2[3]; | |
567 | // | |
568 | AliHelix mhelix(xm,pm,signm); | |
569 | // | |
570 | //find intersection linear | |
571 | // | |
572 | Double_t phase[2][2],radius[2]; | |
573 | Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200); | |
574 | Double_t delta1=10000,delta2=10000; | |
575 | ||
576 | if (points==1){ | |
f007cb5f | 577 | fRs[0] = TMath::Sqrt(radius[0]); |
578 | fRs[1] = TMath::Sqrt(radius[0]); | |
51ad6848 | 579 | } |
580 | if (points==2){ | |
f007cb5f | 581 | fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1])); |
582 | fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1])); | |
51ad6848 | 583 | } |
584 | ||
585 | if (points>0){ | |
586 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
587 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
588 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
589 | } | |
590 | if (points==2){ | |
591 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
592 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
593 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
594 | } | |
595 | if (points==1){ | |
f007cb5f | 596 | fRs[0] = TMath::Sqrt(radius[0]); |
597 | fRs[1] = TMath::Sqrt(radius[0]); | |
51ad6848 | 598 | fDistMinR = delta1; |
599 | } | |
600 | if (points==2){ | |
601 | if (radius[0]<radius[1]){ | |
f007cb5f | 602 | fRs[0] = TMath::Sqrt(radius[0]); |
603 | fRs[1] = TMath::Sqrt(radius[1]); | |
51ad6848 | 604 | fDistMinR = delta1; |
605 | } | |
606 | else{ | |
f007cb5f | 607 | fRs[0] = TMath::Sqrt(radius[1]); |
608 | fRs[1] = TMath::Sqrt(radius[0]); | |
51ad6848 | 609 | fDistMinR = delta2; |
610 | } | |
611 | } | |
612 | // | |
613 | // | |
614 | distance1 = TMath::Min(delta1,delta2); | |
615 | // | |
616 | //find intersection parabolic | |
617 | // | |
618 | points = dhelix1.GetRPHIintersections(mhelix, phase, radius); | |
619 | delta1=10000,delta2=10000; | |
620 | ||
621 | if (points>0){ | |
622 | dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
623 | } | |
624 | if (points==2){ | |
625 | dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
626 | } | |
627 | ||
628 | distance2 = TMath::Min(delta1,delta2); | |
f007cb5f | 629 | if (distance2>100) fDist2 =100; |
630 | return; | |
51ad6848 | 631 | if (delta1<delta2){ |
632 | //get V0 info | |
633 | dhelix1.Evaluate(phase[0][0],fXr); | |
634 | dhelix1.GetMomentum(phase[0][0],fPdr); | |
635 | mhelix.GetMomentum(phase[0][1],fPm); | |
636 | dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle); | |
637 | fRr = TMath::Sqrt(radius[0]); | |
638 | } | |
639 | else{ | |
640 | dhelix1.Evaluate(phase[1][0],fXr); | |
641 | dhelix1.GetMomentum(phase[1][0], fPdr); | |
642 | mhelix.GetMomentum(phase[1][1], fPm); | |
643 | dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle); | |
644 | fRr = TMath::Sqrt(radius[1]); | |
645 | } | |
646 | fDist1 = TMath::Sqrt(distance1); | |
647 | fDist2 = TMath::Sqrt(distance2); | |
648 | ||
649 | if (fDist2<10.5){ | |
650 | Double_t x,alpha,param[5],cov[15]; | |
651 | // | |
652 | fT1.fESDTrack.GetInnerExternalParameters(x,param); | |
653 | fT1.fESDTrack.GetInnerExternalCovariance(cov); | |
654 | alpha = fT1.fESDTrack.GetInnerAlpha(); | |
655 | AliExternalTrackParam paramm(x,alpha,param,cov); | |
656 | // | |
657 | fT2.fESDTrack.GetInnerExternalParameters(x,param); | |
658 | fT2.fESDTrack.GetInnerExternalCovariance(cov); | |
659 | alpha = fT2.fESDTrack.GetInnerAlpha(); | |
660 | AliExternalTrackParam paramd(x,alpha,param,cov); | |
661 | } | |
662 | // | |
663 | // | |
664 | ||
665 | Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; | |
666 | Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]}; | |
667 | ||
668 | Float_t vnorm2 = v[0]*v[0]+v[1]*v[1]; | |
669 | Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2); | |
670 | vnorm2 = TMath::Sqrt(vnorm2); | |
671 | Float_t pnorm2 = p[0]*p[0]+p[1]*p[1]; | |
672 | Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2); | |
673 | pnorm2 = TMath::Sqrt(pnorm2); | |
674 | ||
675 | fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2); | |
676 | fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3); | |
677 | fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3); | |
678 | } | |
679 | } | |
680 | ||
681 | //// | |
682 | void AliESDRecKinkInfo::Update() | |
ae7d73d2 | 683 | { |
51ad6848 | 684 | |
685 | if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){ | |
686 | // | |
687 | // IF BOTH RECONSTRUCTED | |
688 | Float_t distance1,distance2; | |
689 | Double_t xx[3],pp[3]; | |
690 | // | |
691 | Double_t xd[3],pd[3],signd; | |
692 | Double_t xm[3],pm[3],signm; | |
693 | for (Int_t i=0;i<3;i++){ | |
694 | xd[i] = fT2.fTPCinR1[i]; | |
695 | pd[i] = fT2.fTPCinP1[i]; | |
696 | xm[i] = fT1.fTPCinR1[i]; | |
697 | pm[i] = fT1.fTPCinP1[i]; | |
698 | } | |
699 | signd = fT2.fSign<0 ? -1:1; | |
700 | signm = fT1.fSign<0 ? -1:1; | |
701 | ||
702 | AliHelix dhelix1(xd,pd,signd); | |
703 | dhelix1.GetMomentum(0,pp,0); | |
704 | dhelix1.Evaluate(0,xx); | |
705 | // | |
706 | // Double_t x2[3],p2[3]; | |
707 | // | |
708 | AliHelix mhelix(xm,pm,signm); | |
709 | // | |
710 | //find intersection linear | |
711 | // | |
712 | Double_t phase[2][2],radius[2]; | |
713 | Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200); | |
714 | Double_t delta1=10000,delta2=10000; | |
715 | ||
716 | if (points==1){ | |
717 | fMinR = TMath::Sqrt(radius[0]); | |
718 | } | |
719 | if (points==2){ | |
720 | fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1])); | |
721 | } | |
722 | ||
723 | if (points>0){ | |
724 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
725 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
726 | dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
727 | } | |
728 | if (points==2){ | |
729 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
730 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
731 | dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
732 | } | |
733 | if (points==1){ | |
734 | fMinR = TMath::Sqrt(radius[0]); | |
735 | fDistMinR = delta1; | |
736 | } | |
737 | if (points==2){ | |
738 | if (radius[0]<radius[1]){ | |
739 | fMinR = TMath::Sqrt(radius[0]); | |
740 | fDistMinR = delta1; | |
741 | } | |
742 | else{ | |
743 | fMinR = TMath::Sqrt(radius[1]); | |
744 | fDistMinR = delta2; | |
745 | } | |
746 | } | |
747 | // | |
748 | // | |
749 | distance1 = TMath::Min(delta1,delta2); | |
750 | // | |
751 | //find intersection parabolic | |
752 | // | |
753 | points = dhelix1.GetRPHIintersections(mhelix, phase, radius); | |
754 | delta1=10000,delta2=10000; | |
755 | ||
756 | if (points>0){ | |
757 | dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
758 | } | |
759 | if (points==2){ | |
760 | dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
761 | } | |
762 | ||
763 | distance2 = TMath::Min(delta1,delta2); | |
764 | if (delta1<delta2){ | |
765 | //get V0 info | |
766 | dhelix1.Evaluate(phase[0][0],fXr); | |
767 | dhelix1.GetMomentum(phase[0][0],fPdr); | |
768 | mhelix.GetMomentum(phase[0][1],fPm); | |
769 | dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle); | |
770 | fRr = TMath::Sqrt(radius[0]); | |
771 | } | |
772 | else{ | |
773 | dhelix1.Evaluate(phase[1][0],fXr); | |
774 | dhelix1.GetMomentum(phase[1][0], fPdr); | |
775 | mhelix.GetMomentum(phase[1][1], fPm); | |
776 | dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle); | |
777 | fRr = TMath::Sqrt(radius[1]); | |
778 | } | |
779 | fDist1 = TMath::Sqrt(distance1); | |
780 | fDist2 = TMath::Sqrt(distance2); | |
781 | ||
782 | if (fDist2<10.5){ | |
783 | Double_t x,alpha,param[5],cov[15]; | |
784 | // | |
785 | fT1.fESDTrack.GetInnerExternalParameters(x,param); | |
786 | fT1.fESDTrack.GetInnerExternalCovariance(cov); | |
787 | alpha = fT1.fESDTrack.GetInnerAlpha(); | |
788 | AliExternalTrackParam paramm(x,alpha,param,cov); | |
789 | // | |
790 | fT2.fESDTrack.GetInnerExternalParameters(x,param); | |
791 | fT2.fESDTrack.GetInnerExternalCovariance(cov); | |
792 | alpha = fT2.fESDTrack.GetInnerAlpha(); | |
793 | AliExternalTrackParam paramd(x,alpha,param,cov); | |
794 | /* | |
795 | AliESDkink kink; | |
796 | kink.Update(¶mm,¶md); | |
797 | // kink.Dump(); | |
798 | Double_t diff = kink.fRr-fRr; | |
799 | Double_t diff2 = kink.fDist2-fDist2; | |
800 | printf("Diff\t%f\t%f\n",diff,diff2); | |
801 | */ | |
802 | } | |
803 | ||
804 | // | |
805 | // | |
806 | } | |
ae7d73d2 | 807 | |
808 | } | |
809 | ||
810 | ||
811 | //////////////////////////////////////////////////////////////////////// | |
812 | ESDCmpTr::ESDCmpTr() | |
813 | { | |
814 | Reset(); | |
815 | } | |
816 | ||
817 | //////////////////////////////////////////////////////////////////////// | |
818 | ESDCmpTr::ESDCmpTr(const char* fnGenTracks, | |
819 | const char* fnCmp, | |
820 | const char* fnGalice, Int_t direction, | |
821 | Int_t nEvents, Int_t firstEvent) | |
822 | { | |
823 | Reset(); | |
824 | // fFnGenTracks = fnGenTracks; | |
825 | // fFnCmp = fnCmp; | |
826 | sprintf(fFnGenTracks,"%s",fnGenTracks); | |
827 | sprintf(fFnCmp,"%s",fnCmp); | |
828 | ||
829 | fFirstEventNr = firstEvent; | |
830 | fEventNr = firstEvent; | |
831 | fNEvents = nEvents; | |
832 | fDirection = direction; | |
833 | // | |
834 | fLoader = AliRunLoader::Open(fnGalice); | |
835 | if (gAlice){ | |
836 | //delete gAlice->GetRunLoader(); | |
837 | delete gAlice; | |
838 | gAlice = 0x0; | |
839 | } | |
840 | if (fLoader->LoadgAlice()){ | |
841 | cerr<<"Error occured while l"<<endl; | |
842 | } | |
843 | Int_t nall = fLoader->GetNumberOfEvents(); | |
844 | if (nEvents==0) { | |
845 | nEvents =nall; | |
846 | fNEvents=nall; | |
847 | fFirstEventNr=0; | |
848 | } | |
849 | ||
850 | if (nall<=0){ | |
851 | cerr<<"no events available"<<endl; | |
852 | fEventNr = 0; | |
853 | return; | |
854 | } | |
855 | if (firstEvent+nEvents>nall) { | |
856 | fEventNr = nall-firstEvent; | |
857 | cerr<<"restricted number of events availaible"<<endl; | |
858 | } | |
859 | AliMagF * magf = gAlice->Field(); | |
860 | AliTracker::SetFieldMap(magf); | |
861 | ||
862 | } | |
863 | ||
864 | ||
865 | //////////////////////////////////////////////////////////////////////// | |
866 | ESDCmpTr::~ESDCmpTr() | |
867 | { | |
868 | if (fLoader) { | |
869 | delete fLoader; | |
870 | } | |
871 | } | |
872 | ||
873 | ////////////////////////////////////////////////////////////// | |
874 | Int_t ESDCmpTr::SetIO() | |
875 | { | |
876 | // | |
877 | // | |
878 | CreateTreeCmp(); | |
879 | if (!fTreeCmp) return 1; | |
880 | fParamTPC = GetTPCParam(); | |
881 | // | |
882 | if (!ConnectGenTree()) { | |
883 | cerr<<"Cannot connect tree with generated tracks"<<endl; | |
884 | return 1; | |
885 | } | |
886 | return 0; | |
887 | } | |
888 | ||
889 | ////////////////////////////////////////////////////////////// | |
890 | ||
891 | Int_t ESDCmpTr::SetIO(Int_t eventNr) | |
892 | { | |
893 | // | |
894 | // | |
895 | // SET INPUT | |
896 | // | |
897 | TFile f("AliESDs.root"); | |
898 | // | |
899 | ||
900 | TTree* tree = (TTree*) f.Get("esdTree"); | |
901 | if (!tree) { | |
902 | Char_t ename[100]; | |
903 | sprintf(ename,"%d",eventNr); | |
904 | fEvent = (AliESD*)f.Get(ename); | |
905 | if (!fEvent){ | |
906 | sprintf(ename,"ESD%d",eventNr); | |
907 | fEvent = (AliESD*)f.Get(ename); | |
908 | } | |
909 | } | |
910 | else{ | |
911 | tree->SetBranchAddress("ESD", &fEvent); | |
912 | tree->GetEntry(eventNr); | |
913 | } | |
914 | ||
915 | ||
916 | /* | |
917 | Char_t ename[100]; | |
918 | sprintf(ename,"%d",eventNr); | |
919 | fEvent = (AliESD*)f.Get(ename); | |
920 | if (!fEvent){ | |
921 | sprintf(ename,"ESD%d",eventNr); | |
922 | fEvent = (AliESD*)f.Get(ename); | |
923 | } | |
924 | ||
925 | TTree* tree = (TTree*) f.Get("esdTree"); | |
926 | if (!tree) { | |
927 | Error("CheckESD", "no ESD tree found"); | |
928 | return kFALSE; | |
929 | } | |
930 | tree->SetBranchAddress("ESD", &fEvent); | |
931 | tree->GetEntry(eventNr); | |
932 | */ | |
933 | ||
934 | if (!fEvent) return 1; | |
935 | ||
936 | return 0; | |
937 | } | |
938 | ||
939 | ||
940 | ||
941 | //////////////////////////////////////////////////////////////////////// | |
942 | void ESDCmpTr::Reset() | |
943 | { | |
944 | fEventNr = 0; | |
945 | fNEvents = 0; | |
946 | fTreeCmp = 0; | |
51ad6848 | 947 | fTreeCmpKinks =0; |
948 | fTreeCmpV0 =0; | |
ae7d73d2 | 949 | // fFnCmp = "cmpTracks.root"; |
950 | fFileGenTracks = 0; | |
951 | fDebug = 0; | |
952 | // | |
953 | fParamTPC = 0; | |
954 | fEvent =0; | |
955 | } | |
956 | ||
957 | //////////////////////////////////////////////////////////////////////// | |
958 | Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr) | |
959 | { | |
960 | fNEvents = nEvents; | |
961 | fFirstEventNr = firstEventNr; | |
962 | return Exec(); | |
963 | } | |
964 | ||
965 | //////////////////////////////////////////////////////////////////////// | |
966 | Int_t ESDCmpTr::Exec() | |
967 | { | |
968 | TStopwatch timer; | |
969 | timer.Start(); | |
970 | ||
971 | if (SetIO()==1) | |
972 | return 1; | |
973 | ||
974 | fNextTreeGenEntryToRead = 0; | |
51ad6848 | 975 | fNextKinkToRead = 0; |
81e97e0d | 976 | fNextV0ToRead =0; |
ae7d73d2 | 977 | cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl; |
978 | for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents; | |
979 | eventNr++) { | |
980 | fEventNr = eventNr; | |
981 | SetIO(fEventNr); | |
982 | fNParticles = gAlice->GetEvent(fEventNr); | |
983 | ||
f007cb5f | 984 | fIndexRecTracks = new Short_t[fNParticles*20]; //write at maximum 4 tracks corresponding to particle |
985 | fIndexRecKinks = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle | |
986 | fIndexRecV0 = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle | |
51ad6848 | 987 | |
f007cb5f | 988 | fFakeRecTracks = new Short_t[fNParticles]; |
989 | fMultiRecTracks = new Short_t[fNParticles]; | |
990 | fMultiRecKinks = new Short_t[fNParticles]; | |
991 | fMultiRecV0 = new Short_t[fNParticles]; | |
ae7d73d2 | 992 | |
51ad6848 | 993 | for (Int_t i = 0; i<fNParticles; i++) { |
994 | for (Int_t j=0;j<20;j++){ | |
995 | fIndexRecTracks[20*i+j] = -1; | |
996 | fIndexRecKinks[20*i+j] = -1; | |
997 | fIndexRecV0[20*i+j] = -1; | |
998 | } | |
ae7d73d2 | 999 | fFakeRecTracks[i] = 0; |
1000 | fMultiRecTracks[i] = 0; | |
51ad6848 | 1001 | fMultiRecKinks[i] = 0; |
1002 | fMultiRecV0[i] = 0; | |
ae7d73d2 | 1003 | } |
1004 | ||
1005 | cout<<"Start to process event "<<fEventNr<<endl; | |
1006 | cout<<"\tfNParticles = "<<fNParticles<<endl; | |
1007 | if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl; | |
1008 | if (TreeTLoop()>0) return 1; | |
1009 | ||
1010 | if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl; | |
1011 | if (TreeGenLoop(eventNr)>0) return 1; | |
51ad6848 | 1012 | BuildKinkInfo0(eventNr); |
1013 | BuildV0Info(eventNr); | |
1014 | fRecArray->Delete(); | |
1015 | ||
ae7d73d2 | 1016 | if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl; |
1017 | ||
1018 | delete [] fIndexRecTracks; | |
51ad6848 | 1019 | delete [] fIndexRecKinks; |
1020 | delete [] fIndexRecV0; | |
ae7d73d2 | 1021 | delete [] fFakeRecTracks; |
1022 | delete [] fMultiRecTracks; | |
51ad6848 | 1023 | delete [] fMultiRecKinks; |
1024 | delete [] fMultiRecV0; | |
ae7d73d2 | 1025 | } |
1026 | ||
1027 | CloseOutputFile(); | |
1028 | ||
1029 | cerr<<"Exec finished"<<endl; | |
1030 | timer.Stop(); | |
1031 | timer.Print(); | |
1032 | return 0; | |
1033 | ||
1034 | } | |
1035 | //////////////////////////////////////////////////////////////////////// | |
1036 | Bool_t ESDCmpTr::ConnectGenTree() | |
1037 | { | |
1038 | // | |
1039 | // connect all branches from the genTracksTree | |
1040 | // use the same variables as for the new cmp tree, it may work | |
1041 | // | |
1042 | fFileGenTracks = TFile::Open(fFnGenTracks,"READ"); | |
1043 | if (!fFileGenTracks) { | |
1044 | cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl; | |
1045 | return kFALSE; | |
1046 | } | |
1047 | fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree"); | |
1048 | if (!fTreeGenTracks) { | |
1049 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
1050 | <<fFnGenTracks<<endl; | |
1051 | return kFALSE; | |
1052 | } | |
1053 | // | |
1054 | fMCInfo = new AliMCInfo; | |
1055 | fTreeGenTracks->SetBranchAddress("MC",&fMCInfo); | |
1056 | // | |
51ad6848 | 1057 | // |
1058 | fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree"); | |
1059 | if (!fTreeGenKinks) { | |
1060 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
1061 | <<fFnGenTracks<<endl; | |
1062 | //return kFALSE; | |
1063 | } | |
1064 | else{ | |
1065 | fGenKinkInfo = new AliGenKinkInfo; | |
1066 | fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo); | |
1067 | } | |
1068 | ||
1069 | fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree"); | |
1070 | if (!fTreeGenV0) { | |
1071 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
1072 | <<fFnGenTracks<<endl; | |
1073 | //return kFALSE; | |
1074 | } | |
1075 | else{ | |
1076 | fGenV0Info = new AliGenV0Info; | |
1077 | fTreeGenV0->SetBranchAddress("MC",&fGenV0Info); | |
1078 | } | |
1079 | // | |
ae7d73d2 | 1080 | if (fDebug > 1) { |
1081 | cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl; | |
1082 | } | |
1083 | return kTRUE; | |
1084 | } | |
1085 | ||
1086 | ||
1087 | //////////////////////////////////////////////////////////////////////// | |
1088 | void ESDCmpTr::CreateTreeCmp() | |
1089 | { | |
1090 | fFileCmp = TFile::Open(fFnCmp,"RECREATE"); | |
1091 | if (!fFileCmp) { | |
1092 | cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl; | |
1093 | return; | |
1094 | } | |
ae7d73d2 | 1095 | // |
51ad6848 | 1096 | // |
1097 | fTreeCmp = new TTree("ESDcmpTracks","ESDcmpTracks"); | |
ae7d73d2 | 1098 | fMCInfo = new AliMCInfo; |
1099 | fRecInfo = new AliESDRecInfo; | |
ae7d73d2 | 1100 | AliESDtrack * esdTrack = new AliESDtrack; |
51ad6848 | 1101 | // AliITStrackMI * itsTrack = new AliITStrackMI; |
ae7d73d2 | 1102 | fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo); |
1103 | fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo); | |
f007cb5f | 1104 | // fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack); |
51ad6848 | 1105 | // fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack); |
ae7d73d2 | 1106 | delete esdTrack; |
1107 | // | |
51ad6848 | 1108 | // |
1109 | fTreeCmpKinks = new TTree("ESDcmpKinks","ESDcmpKinks"); | |
1110 | fGenKinkInfo = new AliGenKinkInfo; | |
1111 | fRecKinkInfo = new AliESDRecKinkInfo; | |
1112 | fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo); | |
1113 | fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo); | |
1114 | // | |
1115 | // | |
1116 | fTreeCmpV0 = new TTree("ESDcmpV0","ESDcmpV0"); | |
1117 | fGenV0Info = new AliGenV0Info; | |
1118 | fRecV0Info = new AliESDRecV0Info; | |
1119 | fTreeCmpV0->Branch("MC.","AliGenV0Info", &fGenV0Info); | |
1120 | fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info); | |
1121 | // | |
ae7d73d2 | 1122 | fTreeCmp->AutoSave(); |
51ad6848 | 1123 | fTreeCmpKinks->AutoSave(); |
1124 | fTreeCmpV0->AutoSave(); | |
ae7d73d2 | 1125 | } |
1126 | //////////////////////////////////////////////////////////////////////// | |
1127 | void ESDCmpTr::CloseOutputFile() | |
1128 | { | |
1129 | if (!fFileCmp) { | |
1130 | cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl; | |
1131 | return; | |
1132 | } | |
1133 | fFileCmp->cd(); | |
1134 | fTreeCmp->Write(); | |
1135 | delete fTreeCmp; | |
1136 | ||
1137 | fFileCmp->Close(); | |
1138 | delete fFileCmp; | |
1139 | return; | |
1140 | } | |
1141 | //////////////////////////////////////////////////////////////////////// | |
1142 | ||
1143 | TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef, | |
1144 | AliTPCParam *paramTPC) { | |
1145 | ||
1146 | Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; | |
1147 | Int_t index[4]; | |
1148 | paramTPC->Transform0to1(x,index); | |
1149 | paramTPC->Transform1to2(x,index); | |
1150 | return TVector3(x); | |
1151 | } | |
1152 | //////////////////////////////////////////////////////////////////////// | |
1153 | ||
1154 | Int_t ESDCmpTr::TreeTLoop() | |
1155 | { | |
1156 | // | |
1157 | // loop over all ESD reconstructed tracks and store info in memory | |
1158 | // | |
51ad6848 | 1159 | // + loop over all reconstructed kinks |
ae7d73d2 | 1160 | TStopwatch timer; |
1161 | timer.Start(); | |
1162 | // | |
1163 | Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks(); | |
51ad6848 | 1164 | Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks(); |
1165 | Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0MIs(); | |
f007cb5f | 1166 | fSignedKinks = new Short_t[nKinks]; |
1167 | fSignedV0 = new Short_t[nV0MIs]; | |
ae7d73d2 | 1168 | // |
51ad6848 | 1169 | // load kinks to the memory |
1170 | for (Int_t i=0; i<nKinks;i++){ | |
1171 | AliESDkink * kink =fEvent->GetKink(i); | |
1172 | fSignedKinks[i]=0; | |
1173 | if (kink->fStatus<0) continue; | |
1174 | } | |
ae7d73d2 | 1175 | // |
51ad6848 | 1176 | for (Int_t i=0; i<nV0MIs;i++){ |
1177 | AliESDV0MI * v0MI =fEvent->GetV0MI(i); | |
1178 | fSignedV0[i]=0; | |
1179 | if (v0MI->fStatus<0) continue; | |
ae7d73d2 | 1180 | } |
51ad6848 | 1181 | |
ae7d73d2 | 1182 | // |
1183 | // | |
1184 | AliESDtrack * track=0; | |
1185 | for (Int_t iEntry=0; iEntry<nEntries;iEntry++){ | |
51ad6848 | 1186 | //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry); |
1187 | track = (AliESDtrack*)fEvent->GetTrack(iEntry); | |
ae7d73d2 | 1188 | // |
1189 | Int_t label = track->GetLabel(); | |
1190 | Int_t absLabel = abs(label); | |
1191 | if (absLabel < fNParticles) { | |
1192 | // fIndexRecTracks[absLabel] = iEntry; | |
1193 | if (label < 0) fFakeRecTracks[absLabel]++; | |
1194 | if (fMultiRecTracks[absLabel]>0){ | |
51ad6848 | 1195 | if (fMultiRecTracks[absLabel]<20) |
1196 | fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] = iEntry; | |
ae7d73d2 | 1197 | } |
1198 | else | |
51ad6848 | 1199 | fIndexRecTracks[absLabel*20] = iEntry; |
ae7d73d2 | 1200 | fMultiRecTracks[absLabel]++; |
1201 | } | |
51ad6848 | 1202 | } |
1203 | // sort reconstructed kinks | |
1204 | // | |
1205 | AliESDkink * kink=0; | |
1206 | for (Int_t iEntry=0; iEntry<nKinks;iEntry++){ | |
1207 | kink = (AliESDkink*)fEvent->GetKink(iEntry); | |
1208 | if (!kink) continue; | |
1209 | // | |
1210 | Int_t label0 = TMath::Abs(kink->fLab[0]); | |
1211 | Int_t label1 = TMath::Abs(kink->fLab[1]); | |
51ad6848 | 1212 | Int_t absLabel = TMath::Min(label0,label1); |
1213 | if (absLabel < fNParticles) { | |
1214 | if (fMultiRecKinks[absLabel]>0){ | |
1215 | if (fMultiRecKinks[absLabel]<20) | |
1216 | fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] = iEntry; | |
1217 | } | |
1218 | else | |
1219 | fIndexRecKinks[absLabel*20] = iEntry; | |
1220 | fMultiRecKinks[absLabel]++; | |
1221 | } | |
ae7d73d2 | 1222 | } |
51ad6848 | 1223 | // --sort reconstructed V0 |
1224 | // | |
1225 | AliESDV0MI * v0MI=0; | |
1226 | for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){ | |
1227 | v0MI = fEvent->GetV0MI(iEntry); | |
1228 | if (!v0MI) continue; | |
1229 | // | |
f007cb5f | 1230 | // Int_t label0 = TMath::Abs(v0MI->fLab[0]); |
1231 | //Int_t label1 = TMath::Abs(v0MI->fLab[1]); | |
1232 | // | |
1233 | for (Int_t i=0;i<2;i++){ | |
1234 | // Int_t absLabel = TMath::Min(label0,label1); | |
1235 | Int_t absLabel = TMath::Abs(v0MI->fLab[i]); | |
1236 | if (absLabel < fNParticles) { | |
1237 | if (fMultiRecV0[absLabel]>0){ | |
1238 | if (fMultiRecV0[absLabel]<20) | |
1239 | fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] = iEntry; | |
1240 | } | |
1241 | else | |
1242 | fIndexRecV0[absLabel*20] = iEntry; | |
1243 | fMultiRecV0[absLabel]++; | |
51ad6848 | 1244 | } |
51ad6848 | 1245 | } |
1246 | } | |
1247 | ||
1248 | ||
ae7d73d2 | 1249 | printf("Time spended in TreeTLoop\n"); |
1250 | timer.Print(); | |
1251 | ||
1252 | if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl; | |
1253 | return 0; | |
1254 | } | |
1255 | ||
1256 | //////////////////////////////////////////////////////////////////////// | |
1257 | Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr) | |
1258 | { | |
1259 | // | |
1260 | // loop over all entries for a given event, find corresponding | |
1261 | // rec. track and store in the fTreeCmp | |
1262 | // | |
1263 | TStopwatch timer; | |
1264 | timer.Start(); | |
1265 | Int_t entry = fNextTreeGenEntryToRead; | |
1266 | Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast(); | |
1267 | cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" " | |
1268 | <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl; | |
1269 | TBranch * branch = fTreeCmp->GetBranch("RC"); | |
1270 | branch->SetAddress(&fRecInfo); // set all pointers | |
51ad6848 | 1271 | fRecArray = new TObjArray(fNParticles); |
ae7d73d2 | 1272 | |
1273 | while (entry < nParticlesTR) { | |
1274 | fTreeGenTracks->GetEntry(entry); | |
1275 | entry++; | |
1276 | if (eventNr < fMCInfo->fEventNr) continue; | |
1277 | if (eventNr > fMCInfo->fEventNr) continue;; | |
1278 | // | |
1279 | fNextTreeGenEntryToRead = entry-1; | |
1280 | if (fDebug > 2 && fMCInfo->fLabel < 10) { | |
1281 | cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl; | |
1282 | } | |
51ad6848 | 1283 | // if (fMCInfo->fNTPCRef<1) continue; // not TPCref |
ae7d73d2 | 1284 | // |
1285 | fRecInfo->Reset(); | |
1286 | AliESDtrack * track=0; | |
1287 | fRecInfo->fReconstructed =0; | |
1288 | TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC); | |
1289 | local.GetXYZ(fRecInfo->fTRLocalCoord); | |
1290 | // | |
51ad6848 | 1291 | if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) { |
1292 | //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]); | |
1293 | track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]); | |
ae7d73d2 | 1294 | // |
1295 | // | |
1296 | // find nearest track if multifound | |
51ad6848 | 1297 | //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge); |
1298 | // | |
1299 | Int_t status = 0; | |
1300 | if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++; | |
1301 | if ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++; | |
1302 | if ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++; | |
1303 | ||
1304 | // | |
1305 | if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){ | |
1306 | // | |
1307 | Double_t p[3]; | |
1308 | track->GetInnerPxPyPz(p); | |
1309 | Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; | |
1310 | // | |
1311 | for (Int_t i=1;i<20;i++){ | |
1312 | if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){ | |
1313 | AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]); | |
1314 | if (!track2) continue; | |
1315 | //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; // | |
1316 | //if (sign2<0) continue; | |
1317 | track2->GetInnerPxPyPz(p); | |
1318 | Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; | |
1319 | /* | |
1320 | if (sign<0){ | |
1321 | sign = sign2; | |
1322 | track = track2; | |
1323 | maxp = mom; | |
1324 | continue; | |
1325 | } | |
1326 | */ | |
1327 | // | |
1328 | Int_t status2 = 0; | |
1329 | if ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++; | |
1330 | if ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++; | |
1331 | if ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++; | |
1332 | if (status2<status) continue; | |
1333 | // | |
1334 | if (mom<maxp) continue; | |
1335 | maxp = mom; | |
1336 | track = track2; | |
1337 | // | |
ae7d73d2 | 1338 | } |
1339 | } | |
1340 | } | |
1341 | // | |
1342 | fRecInfo->fESDTrack =*track; | |
1343 | if (track->GetITStrack()) | |
51ad6848 | 1344 | fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack()); |
81e97e0d | 1345 | else{ |
1346 | fRecInfo->fITStrack = *track; | |
1347 | } | |
ae7d73d2 | 1348 | if (track->GetTRDtrack()){ |
1349 | fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack()); | |
1350 | } | |
1351 | else{ | |
1352 | fRecInfo->fTRDtrack.SetdEdx(-1); | |
1353 | } | |
1354 | fRecInfo->fReconstructed = 1; | |
1355 | fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel]; | |
1356 | fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel]; | |
1357 | // | |
1358 | fRecInfo->Update(fMCInfo,fParamTPC,kTRUE); | |
1359 | } | |
1360 | else{ | |
1361 | fRecInfo->Update(fMCInfo,fParamTPC,kFALSE); | |
1362 | } | |
51ad6848 | 1363 | fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel); |
ae7d73d2 | 1364 | fTreeCmp->Fill(); |
1365 | } | |
1366 | fTreeCmp->AutoSave(); | |
51ad6848 | 1367 | //fTracks->Delete(); |
ae7d73d2 | 1368 | printf("Time spended in TreeGenLoop\n"); |
1369 | timer.Print(); | |
1370 | if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl; | |
1371 | ||
1372 | return 0; | |
1373 | } | |
51ad6848 | 1374 | |
1375 | ||
1376 | ||
1377 | //////////////////////////////////////////////////////////////////////// | |
1378 | //////////////////////////////////////////////////////////////////////// | |
1379 | //////////////////////////////////////////////////////////////////////// | |
1380 | Int_t ESDCmpTr::BuildKinkInfo0(Int_t eventNr) | |
1381 | { | |
1382 | // | |
1383 | // loop over all entries for a given event, find corresponding | |
1384 | // rec. track and store in the fTreeCmp | |
1385 | // | |
1386 | TStopwatch timer; | |
1387 | timer.Start(); | |
1388 | Int_t entry = fNextKinkToRead; | |
1389 | Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast(); | |
1390 | cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" " | |
1391 | <<nParticlesTR<<" "<<fNextKinkToRead<<endl; | |
1392 | // | |
1393 | TBranch * branch = fTreeCmpKinks->GetBranch("RC."); | |
1394 | branch->SetAddress(&fRecKinkInfo); // set all pointers | |
1395 | ||
1396 | // | |
1397 | while (entry < nParticlesTR) { | |
1398 | fTreeGenKinks->GetEntry(entry); | |
1399 | entry++; | |
1400 | if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue; | |
1401 | if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;; | |
1402 | // | |
1403 | fNextKinkToRead = entry-1; | |
1404 | // | |
1405 | // | |
1406 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel); | |
1407 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel); | |
1408 | fRecKinkInfo->fT1 = (*fRecInfo1); | |
1409 | fRecKinkInfo->fT2 = (*fRecInfo2); | |
1410 | fRecKinkInfo->fStatus =0; | |
1411 | if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1; | |
1412 | if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2; | |
1413 | if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1; | |
1414 | ||
1415 | if (fRecKinkInfo->fStatus==3){ | |
1416 | fRecKinkInfo->Update(); | |
1417 | } | |
1418 | Int_t label = TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel); | |
1419 | Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel); | |
1420 | ||
1421 | AliESDkink *kink=0; | |
1422 | fRecKinkInfo->fRecStatus =0; | |
1423 | fRecKinkInfo->fMultiple = fMultiRecKinks[label]; | |
1424 | fRecKinkInfo->fKinkMultiple=0; | |
1425 | // | |
1426 | if (fMultiRecKinks[label]>0){ | |
1427 | ||
1428 | // for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){ | |
f007cb5f | 1429 | for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){ |
51ad6848 | 1430 | Int_t index = fIndexRecKinks[label*20+j]; |
1431 | //AliESDkink *kink2 = (AliESDkink*)fKinks->At(index); | |
1432 | AliESDkink *kink2 = (AliESDkink*)fEvent->GetKink(index); | |
1433 | if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) { | |
51ad6848 | 1434 | fRecKinkInfo->fKinkMultiple++; |
1435 | fSignedKinks[index]=1; | |
f007cb5f | 1436 | Int_t c0=0; |
1437 | if (kink){ | |
1438 | // if (kink->fTRDOn) c0++; | |
1439 | //if (kink->fITSOn) c0++; | |
1440 | if (kink->GetStatus(2)>0) c0++; | |
1441 | if (kink->GetStatus(0)>0) c0++; | |
1442 | } | |
1443 | Int_t c2=0; | |
1444 | // if (kink2->fTRDOn) c2++; | |
1445 | //if (kink2->fITSOn) c2++; | |
1446 | if (kink2->GetStatus(2)>0) c2++; | |
1447 | if (kink2->GetStatus(0)>0) c2++; | |
1448 | ||
1449 | if (c2<c0) continue; | |
1450 | kink =kink2; | |
51ad6848 | 1451 | } |
1452 | if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) { | |
51ad6848 | 1453 | fRecKinkInfo->fKinkMultiple++; |
1454 | fSignedKinks[index]=1; | |
f007cb5f | 1455 | Int_t c0=0; |
1456 | if (kink){ | |
1457 | //if (kink->fTRDOn) c0++; | |
1458 | //if (kink->fITSOn) c0++; | |
1459 | if (kink->GetStatus(2)>0) c0++; | |
1460 | if (kink->GetStatus(0)>0) c0++; | |
1461 | ||
1462 | } | |
1463 | Int_t c2=0; | |
1464 | // if (kink2->fTRDOn) c2++; | |
1465 | //if (kink2->fITSOn) c2++; | |
1466 | if (kink2->GetStatus(2)>0) c2++; | |
1467 | if (kink2->GetStatus(0)>0) c2++; | |
1468 | ||
1469 | if (c2<c0) continue; | |
1470 | kink =kink2; | |
51ad6848 | 1471 | } |
1472 | } | |
1473 | } | |
1474 | if (kink){ | |
1475 | fRecKinkInfo->fKink = *kink; | |
1476 | fRecKinkInfo->fRecStatus=1; | |
1477 | } | |
1478 | fTreeCmpKinks->Fill(); | |
1479 | } | |
1480 | // Int_t nkinks = fKinks->GetEntriesFast(); | |
1481 | Int_t nkinks = fEvent->GetNumberOfKinks(); | |
1482 | for (Int_t i=0;i<nkinks;i++){ | |
1483 | if (fSignedKinks[i]==0){ | |
1484 | // AliESDkink *kink = (AliESDkink*)fKinks->At(i); | |
1485 | AliESDkink *kink = (AliESDkink*)fEvent->GetKink(i); | |
1486 | if (!kink) continue; | |
1487 | // | |
1488 | fRecKinkInfo->fKink = *kink; | |
1489 | fRecKinkInfo->fRecStatus =-2; | |
1490 | // | |
1491 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0])); | |
1492 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1])); | |
1493 | if (fRecInfo1 && fRecInfo2){ | |
1494 | fRecKinkInfo->fT1 = (*fRecInfo1); | |
1495 | fRecKinkInfo->fT2 = (*fRecInfo2); | |
1496 | fRecKinkInfo->fRecStatus =-1; | |
1497 | } | |
1498 | fTreeCmpKinks->Fill(); | |
1499 | } | |
1500 | } | |
1501 | ||
1502 | ||
1503 | fTreeCmpKinks->AutoSave(); | |
1504 | printf("Time spended in BuilKinkInfo Loop\n"); | |
1505 | timer.Print(); | |
1506 | if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl; | |
1507 | return 0; | |
1508 | } | |
1509 | ||
f007cb5f | 1510 | |
1511 | ||
1512 | void ESDCmpTr::MakePoints(AliESDtrack * track, AliPointsMI &points) | |
1513 | { | |
1514 | // | |
1515 | // make points in global coordinate frame | |
1516 | // | |
1517 | return; | |
1518 | /* | |
1519 | UInt_t itscl[10]; | |
1520 | Int_t nits = track->GetITSclusters(itscl); | |
1521 | Int_t tpccl[1000]; | |
1522 | Int_t ntpc = track->GetTPcclusters(tpccl); | |
1523 | UInt_t trdcl[1000]; | |
1524 | Int_t ntrd = track->GetTRDclusters(trdcl); | |
1525 | // | |
1526 | AliLoader *itsloader = fLoader->GetLoader("ITSLoader"); | |
1527 | AliLoader *tpcloader = fLoader->GetLoader("TPCLoader"); | |
1528 | AliLoader *trdloader = fLoader->GetLoader("TRDLoader"); | |
1529 | // | |
1530 | AliITStrackerMI itstracker(); | |
1531 | */ | |
1532 | } | |
1533 | ||
1534 | ||
51ad6848 | 1535 | //////////////////////////////////////////////////////////////////////// |
1536 | //////////////////////////////////////////////////////////////////////// | |
1537 | //////////////////////////////////////////////////////////////////////// | |
1538 | ||
f007cb5f | 1539 | |
1540 | ||
51ad6848 | 1541 | Int_t ESDCmpTr::BuildV0Info(Int_t eventNr) |
1542 | { | |
1543 | // | |
1544 | // loop over all entries for a given event, find corresponding | |
1545 | // rec. track and store in the fTreeCmp | |
1546 | // | |
1547 | TStopwatch timer; | |
1548 | timer.Start(); | |
1549 | Int_t entry = fNextV0ToRead; | |
1550 | Double_t nParticlesTR = fTreeGenV0->GetEntriesFast(); | |
1551 | cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" " | |
1552 | <<nParticlesTR<<" "<<fNextV0ToRead<<endl; | |
1553 | // | |
1554 | TBranch * branch = fTreeCmpV0->GetBranch("RC."); | |
1555 | branch->SetAddress(&fRecV0Info); // set all pointers | |
1556 | const AliESDVertex * esdvertex = fEvent->GetVertex(); | |
1557 | Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()}; | |
1558 | ||
1559 | // | |
1560 | while (entry < nParticlesTR) { | |
1561 | fTreeGenV0->GetEntry(entry); | |
1562 | entry++; | |
1563 | if (eventNr < fGenV0Info->fMCm.fEventNr) continue; | |
1564 | if (eventNr > fGenV0Info->fMCm.fEventNr) continue;; | |
1565 | // | |
1566 | fNextV0ToRead = entry-1; | |
1567 | // | |
1568 | // | |
1569 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel); | |
1570 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel); | |
1571 | if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue; // interactions | |
1572 | if (!fRecInfo1 || !fRecInfo2) continue; | |
1573 | fRecV0Info->fT1 = (*fRecInfo1); | |
1574 | fRecV0Info->fT2 = (*fRecInfo2); | |
1575 | fRecV0Info->fV0Status =0; | |
1576 | if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1; | |
1577 | if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2; | |
1578 | ||
1579 | if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1; | |
1580 | ||
1581 | ||
1582 | if (abs(fRecV0Info->fV0Status)==3){ | |
1583 | fRecV0Info->Update(vertex); | |
1584 | { | |
1585 | // | |
1586 | // TPC V0 Info | |
1587 | Double_t x,alpha, param[5],cov[15]; | |
1588 | fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(x,param); | |
1589 | fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov); | |
1590 | alpha = fRecV0Info->fT1.fESDTrack.GetInnerAlpha(); | |
1591 | AliExternalTrackParam paramP(x,alpha,param,cov); | |
1592 | // | |
1593 | fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(x,param); | |
1594 | fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov); | |
1595 | alpha = fRecV0Info->fT2.fESDTrack.GetInnerAlpha(); | |
1596 | AliExternalTrackParam paramM(x,alpha,param,cov); | |
1597 | // | |
1598 | fRecV0Info->fV0tpc.SetM(paramM); | |
1599 | fRecV0Info->fV0tpc.SetP(paramP); | |
1600 | Double_t pid1[5],pid2[5]; | |
1601 | fRecV0Info->fT1.fESDTrack.GetESDpid(pid1); | |
1602 | fRecV0Info->fT1.fESDTrack.GetESDpid(pid2); | |
1603 | // | |
1604 | fRecV0Info->fV0tpc.UpdatePID(pid1,pid2); | |
1605 | fRecV0Info->fV0tpc.Update(vertex); | |
1606 | // | |
1607 | // | |
1608 | fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param); | |
1609 | fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov); | |
1610 | alpha = fRecV0Info->fT1.fESDTrack.GetAlpha(); | |
1611 | new (¶mP) AliExternalTrackParam(x,alpha,param,cov); | |
1612 | // | |
1613 | fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param); | |
1614 | fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov); | |
1615 | alpha = fRecV0Info->fT2.fESDTrack.GetAlpha(); | |
1616 | new (¶mM) AliExternalTrackParam(x,alpha,param,cov); | |
1617 | // | |
1618 | fRecV0Info->fV0its.SetM(paramM); | |
1619 | fRecV0Info->fV0its.SetP(paramP); | |
1620 | fRecV0Info->fV0its.UpdatePID(pid1,pid2); | |
1621 | fRecV0Info->fV0its.Update(vertex); | |
1622 | ||
1623 | } | |
1624 | if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){ | |
1625 | if (fRecV0Info->fDist2>10){ | |
1626 | fRecV0Info->Update(vertex); | |
1627 | } | |
1628 | if (fRecV0Info->fDist2>10){ | |
1629 | fRecV0Info->Update(vertex); | |
1630 | } | |
1631 | } | |
1632 | } | |
1633 | // | |
1634 | // take the V0 from reconstruction | |
1635 | ||
1636 | Int_t label = TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel); | |
1637 | Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel); | |
1638 | AliESDV0MI *v0MI=0; | |
1639 | fRecV0Info->fRecStatus =0; | |
1640 | fRecV0Info->fMultiple = fMultiRecV0[label]; | |
1641 | fRecV0Info->fV0Multiple=0; | |
1642 | // | |
f007cb5f | 1643 | if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){ |
51ad6848 | 1644 | |
1645 | // for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){ | |
f007cb5f | 1646 | for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){ |
51ad6848 | 1647 | Int_t index = fIndexRecV0[label*20+j]; |
1648 | AliESDV0MI *v0MI2 = fEvent->GetV0MI(index); | |
1649 | if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) { | |
1650 | v0MI =v0MI2; | |
1651 | fRecV0Info->fV0Multiple++; | |
1652 | fSignedV0[index]=1; | |
1653 | } | |
1654 | if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) { | |
1655 | v0MI =v0MI2; | |
1656 | fRecV0Info->fV0Multiple++; | |
1657 | fSignedV0[index]=1; | |
1658 | } | |
1659 | } | |
1660 | } | |
1661 | if (v0MI){ | |
1662 | fRecV0Info->fV0rec = *v0MI; | |
1663 | fRecV0Info->fRecStatus=1; | |
1664 | } | |
1665 | ||
1666 | fTreeCmpV0->Fill(); | |
1667 | } | |
1668 | // | |
1669 | // write fake v0s | |
1670 | ||
1671 | Int_t nV0MIs = fEvent->GetNumberOfV0MIs(); | |
1672 | for (Int_t i=0;i<nV0MIs;i++){ | |
1673 | if (fSignedV0[i]==0){ | |
1674 | AliESDV0MI *v0MI = fEvent->GetV0MI(i); | |
1675 | if (!v0MI) continue; | |
1676 | // | |
1677 | fRecV0Info->fV0rec = *v0MI; | |
1678 | fRecV0Info->fV0Status =-10; | |
1679 | fRecV0Info->fRecStatus =-2; | |
1680 | // | |
1681 | AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0])); | |
1682 | AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1])); | |
1683 | if (fRecInfo1 && fRecInfo2){ | |
1684 | fRecV0Info->fT1 = (*fRecInfo1); | |
1685 | fRecV0Info->fT2 = (*fRecInfo2); | |
1686 | fRecV0Info->fRecStatus =-1; | |
1687 | } | |
f007cb5f | 1688 | fRecV0Info->Update(vertex); |
51ad6848 | 1689 | fTreeCmpV0->Fill(); |
1690 | } | |
1691 | } | |
1692 | ||
1693 | ||
1694 | ||
1695 | fTreeCmpV0->AutoSave(); | |
1696 | fRecArray->Delete(); | |
1697 | printf("Time spended in BuilV0Info Loop\n"); | |
1698 | timer.Print(); | |
1699 | if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl; | |
1700 | return 0; | |
1701 | } | |
ae7d73d2 | 1702 | //////////////////////////////////////////////////////////////////////// |
1703 | //////////////////////////////////////////////////////////////////////// | |
1704 | ||
1705 | void AliESDComparisonDraw::SetIO(const char *fname) | |
1706 | { | |
1707 | // | |
1708 | TFile* file = TFile::Open(fname); | |
1709 | // | |
1710 | fTree = (TTree*) file->Get("ESDcmpTracks"); | |
1711 | if (!fTree) { | |
1712 | printf("no track comparison tree found\n"); | |
1713 | file->Close(); | |
1714 | delete file; | |
1715 | } | |
1716 | } | |
1717 | ||
1718 |