1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Time Projection Chamber //
20 // Comparison macro for ESD //
22 // marian.ivanov@cern.ch //
32 gSystem->Load("libPWG1.so");
34 AliESDComparisonMI *t2 = new AliESDComparisonMI("genTracks.root","cmpESDTracks.root","galice.root",-1,0,0);
38 TFile f("cmpESDTracks.root");
40 TTree * tree = (TTree*)f.Get("ESDcmpTracks");
46 //some cuts definition
47 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
48 //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
49 //TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
50 TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
51 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
54 TCut crec("crec","fReconstructed==1");
55 TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
56 TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
58 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
59 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
60 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
61 TCut cchi2("cchi2","fESDtrack.fITSchi2MIP[0]<7.&&fESDtrack.fITSchi2MIP[1]<5.&&fESDtrack.fITSchi2MIP[2]<7.&&fESDtrack.fITSchi2MIP[3]<7.5&&fESDtrack.fITSchi2MIP[4]<6.")
66 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
67 comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
68 comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
69 comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
70 comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
71 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
72 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
75 TH1F his("his","his",100,0,20);
76 TH1F hpools("hpools","hpools",100,-7,7);
77 TH1F hfake("hfake","hfake",1000,0,150);
78 TProfile profp0("profp0","profp0",20,-0.4,0.9)
80 comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
84 comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
87 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
90 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio<0.1",10,0.2,1.5)
92 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio>0.1",10,0.2,1.5)
95 comp.fTree->Draw("fESDtrack.fITSsignal/fESDtrack.fTPCsignal","fITSOn&&fTPCOn&&fESDtrack.fITSFakeRatio==0")
97 TH1F his("his","his",100,0,20);
98 TH1F hpools("hpools","hpools",100,-7,7);
100 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
101 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
102 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
103 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
105 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
106 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
107 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
108 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
111 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1")
112 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1")
113 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1")
114 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1")
116 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
117 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
118 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
119 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
121 hdedx3->SetXTitle("P(GeV/c)");
122 hdedx3->SetYTitle("dEdx(unit)");
123 hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
125 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
127 TProfile prof("prof","prof",10,0.5,5);
144 #include "TBenchmark.h"
145 #include "TStopwatch.h"
146 #include "TParticle.h"
149 #include "TVector3.h"
160 #include "AliStack.h"
161 #include "AliESDtrack.h"
162 #include "AliSimDigits.h"
163 #include "AliTPCParam.h"
165 #include "AliTPCLoader.h"
166 #include "AliDetector.h"
167 #include "AliTrackReference.h"
169 #include "AliTPCParamSR.h"
170 #include "AliTracker.h"
171 #include "AliComplexCluster.h"
174 #include "AliESDfriend.h"
175 #include "AliESDtrack.h"
176 #include "AliTPCseed.h"
177 #include "AliITStrackMI.h"
178 #include "AliTRDtrack.h"
179 #include "AliHelix.h"
180 #include "AliESDVertex.h"
181 #include "AliExternalTrackParam.h"
182 #include "AliESDkink.h"
183 #include "AliESDv0.h"
186 #include "AliTreeDraw.h"
187 #include "AliGenInfo.h"
188 #include "AliESDComparisonMI.h"
192 ClassImp(AliESDRecInfo)
193 ClassImp(AliESDRecV0Info)
194 ClassImp(AliESDRecKinkInfo)
195 ClassImp(AliESDComparisonMI)
201 void MakeAliases(AliTreeDraw&comp)
204 // aliases definition
206 comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
207 comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
208 comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
209 comp.T()->SetAlias("theta","MC.fTrackRef.Theta()");
210 comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
211 comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
212 comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
214 comp.T()->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.");
216 comp.T()->SetAlias("dtofmc2","fESDtrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
217 comp.T()->SetAlias("dtofrc2","(fESDtrack.fTrackTime[2]-fESDtrack.fTOFsignal)");
219 comp.T()->SetAlias("psum","fESDtrack.fTOFr[4]+fESDtrack.fTOFr[3]+fESDtrack.fTOFr[2]+fESDtrack.fTOFr[1]+fESDtrack.fTOFr[0]");
220 comp.T()->SetAlias("P0","fESDtrack.fTOFr[0]/psum");
221 comp.T()->SetAlias("P1","fESDtrack.fTOFr[1]/psum");
222 comp.T()->SetAlias("P2","fESDtrack.fTOFr[2]/psum");
223 comp.T()->SetAlias("P3","fESDtrack.fTOFr[3]/psum");
224 comp.T()->SetAlias("P4","fESDtrack.fTOFr[4]/psum");
225 comp.T()->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
229 AliESDRecInfo::AliESDRecInfo():
230 fITSOn(0), // ITS refitted inward
231 fTRDOn(0), // ITS refitted inward
232 fDeltaP(0), //delta of momenta
234 fReconstructed(0), //flag if track was reconstructed
235 fFake(0), // fake track
236 fMultiple(0), // number of reconstructions
237 fTPCOn(0), // TPC refitted inward
238 fBestTOFmatch(0), //best matching between times
239 fESDtrack(0), // esd track
240 fTrackF(0), // friend track
241 fTPCtrack(0), // tpc track
242 fITStrack(0), // its track
243 fTRDtrack(0) // trd track
246 // default constructor
251 AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
257 memcpy(this,&recinfo, sizeof(recinfo));
258 fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
259 SetESDtrack(recinfo.GetESDtrack());
263 AliESDRecInfo::~AliESDRecInfo()
269 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
270 if (fTrackF) { delete fTrackF; fTrackF=0;}
271 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
272 if (fITStrack) { delete fITStrack; fITStrack=0;}
273 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
279 void AliESDRecInfo::Reset()
287 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
288 if (fTrackF) { delete fTrackF; fTrackF=0;}
289 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
290 if (fITStrack) { delete fITStrack; fITStrack=0;}
291 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
294 void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
298 if (fESDtrack) delete fESDtrack;
299 fESDtrack = (AliESDtrack*)track->Clone();
300 if (track->GetFriendTrack()){
301 if (fTrackF) delete fTrackF;
302 fTrackF = (AliESDfriendTrack*)track->GetFriendTrack()->Clone();
307 void AliESDRecInfo::UpdatePoints(AliESDtrack*track)
311 Int_t iclusters[200];
312 Float_t density[160];
313 for (Int_t i=0;i<160;i++) density[i]=-1.;
317 if (fTPCPoints[0]<fTPCPoints[1]) return;
318 // Int_t nclusters=track->GetTPCclusters(iclusters);
324 for (Int_t i=0;i<160;i++){
325 Int_t last = i-range;
326 if (nall<range) nall++;
328 if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
329 if (iclusters[last]==-1) undeff--;
331 if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++;
332 if (iclusters[i]==-1) undeff++;
333 if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
337 for (Int_t i=0;i<160;i++){
338 if (density[i]<0) continue;
339 if (density[i]>maxdens){
346 fTPCPoints[3] = maxdens;
347 fTPCPoints[1] = indexmax;
350 for (Int_t i=indexmax;i<160;i++){
351 if (density[i]<0) continue;
352 if (density[i]<maxdens/2.) {
359 for (Int_t i=indexmax;i>0;i--){
360 if (density[i]<0) continue;
361 if (density[i]<maxdens/2.) {
367 // Density at the last 30 padrows
372 for (Int_t i=159;i>0;i--){
373 if (iclusters[i]==-1) continue; //dead zone
375 if (iclusters[i]>0) ngood++;
378 fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
380 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
387 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESD */*event*/)
391 //calculates derived variables
394 UpdatePoints(fESDtrack);
396 AliTrackReference * ref = &(info->fTrackRef);
397 fTPCinR0[0] = info->fTrackRef.X();
398 fTPCinR0[1] = info->fTrackRef.Y();
399 fTPCinR0[2] = info->fTrackRef.Z();
400 fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
401 fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
403 fTPCinP0[0] = ref->Px();
404 fTPCinP0[1] = ref->Py();
405 fTPCinP0[2] = ref->Pz();
406 fTPCinP0[3] = ref->Pt();
407 fTPCinP0[4] = ref->P();
408 fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P();
411 if (fTPCinP0[3]>0.0000001){
413 fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
414 fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
418 fITSinP0[0]=info->fParticle.Px();
419 fITSinP0[1]=info->fParticle.Py();
420 fITSinP0[2]=info->fParticle.Pz();
421 fITSinP0[3]=info->fParticle.Pt();
423 fITSinR0[0]=info->fParticle.Vx();
424 fITSinR0[1]=info->fParticle.Vy();
425 fITSinR0[2]=info->fParticle.Vz();
426 fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
427 fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
430 if (fITSinP0[3]>0.0000001){
431 fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
432 fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
435 for (Int_t i=0;i<4;i++) fStatus[i] =0;
436 fReconstructed = kFALSE;
440 if (reconstructed==kFALSE) return;
442 fLabels[0] = info->fLabel;
443 fLabels[1] = info->fPrimPart;
444 fReconstructed = kTRUE;
445 fTPCOn = ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
446 fITSOn = ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
447 fTRDOn = ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
450 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0){
454 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCout)>0){
458 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCin)>0)
463 if ((fESDtrack->GetStatus()&AliESDtrack::kITSout)>0){
467 if ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0){
477 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0){
481 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDout)>0){
485 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDStop)>0){
492 if (((fESDtrack->GetStatus()&AliESDtrack::kTOFout)>0)){
496 fESDtrack->GetIntegratedTimes(times);
497 for (Int_t i=0;i<5;i++){
498 if ( TMath::Abs(fESDtrack->GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
499 fBestTOFmatch = fESDtrack->GetTOFsignal()-times[i];
503 fESDtrack->GetTOFLabel(toflabel);
504 Bool_t toffake=kTRUE;
505 Bool_t tofdaughter=kFALSE;
506 for (Int_t i=0;i<3;i++){
507 if (toflabel[i]<0) continue;
508 if (toflabel[i]== TMath::Abs(fESDtrack->GetLabel())) toffake=kFALSE;
509 if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle
512 if (toffake) fStatus[3] =3; //total fake
513 if (tofdaughter) fStatus[3]=2; //fake because of decay
519 if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
521 fESDtrack->GetInnerXYZ(fTPCinR1);
522 fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
523 fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
524 fESDtrack->GetInnerPxPyPz(fTPCinP1);
525 fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
526 fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
529 if (fTPCinP1[3]>0.000000000000001){
530 fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
531 fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
533 Double_t cov[15], param[5],x, alpha;
534 fESDtrack->GetInnerExternalCovariance(cov);
535 fESDtrack->GetInnerExternalParameters(alpha, x,param);
538 fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
539 fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
540 fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
541 fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
542 fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
543 fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
544 fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
545 fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
546 fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
547 Double_t sign = (param[4]>0)? 1.:-1;
549 fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
554 fESDtrack->GetExternalParameters(x,param);
555 // fESDtrack->GetConstrainedExternalParameters(x,param);
557 fESDtrack->GetExternalCovariance(cov);
558 //fESDtrack->GetConstrainedExternalCovariance(cov);
559 if (TMath::Abs(param[4])<0.0000000001) return;
561 fESDtrack->GetXYZ(fITSinR1);
562 fESDtrack->GetPxPyPz(fITSinP1);
563 fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
565 fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
566 fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
569 if (fITSinP1[3]>0.0000001){
570 fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
571 fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
575 fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi
576 fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
577 fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z
578 fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
579 fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
580 fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
581 fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
582 fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
583 fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
584 Double_t sign = (param[4]>0) ? 1:-1;
586 fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
592 void AliESDRecV0Info::Update(Float_t vertex[3])
595 if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
596 Float_t distance1,distance2;
597 Double_t xx[3],pp[3];
599 Double_t xd[3],pd[3],signd;
600 Double_t xm[3],pm[3],signm;
603 if (fT1.fITSOn&&fT2.fITSOn){
604 for (Int_t i=0;i<3;i++){
605 xd[i] = fT2.fITSinR1[i];
606 pd[i] = fT2.fITSinP1[i];
607 xm[i] = fT1.fITSinR1[i];
608 pm[i] = fT1.fITSinP1[i];
613 for (Int_t i=0;i<3;i++){
614 xd[i] = fT2.fTPCinR1[i];
615 pd[i] = fT2.fTPCinP1[i];
616 xm[i] = fT1.fTPCinR1[i];
617 pm[i] = fT1.fTPCinP1[i];
622 signd = fT2.fSign<0 ? -1:1;
623 signm = fT1.fSign<0 ? -1:1;
625 AliHelix dhelix1(xd,pd,signd);
626 dhelix1.GetMomentum(0,pp,0);
627 dhelix1.Evaluate(0,xx);
629 // Double_t x2[3],p2[3];
631 AliHelix mhelix(xm,pm,signm);
633 //find intersection linear
635 Double_t phase[2][2],radius[2];
636 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
637 Double_t delta1=10000,delta2=10000;
640 fRs[0] = TMath::Sqrt(radius[0]);
641 fRs[1] = TMath::Sqrt(radius[0]);
644 fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
645 fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
649 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
650 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
651 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
654 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
655 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
656 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
659 fRs[0] = TMath::Sqrt(radius[0]);
660 fRs[1] = TMath::Sqrt(radius[0]);
664 if (radius[0]<radius[1]){
665 fRs[0] = TMath::Sqrt(radius[0]);
666 fRs[1] = TMath::Sqrt(radius[1]);
670 fRs[0] = TMath::Sqrt(radius[1]);
671 fRs[1] = TMath::Sqrt(radius[0]);
677 distance1 = TMath::Min(delta1,delta2);
679 //find intersection parabolic
681 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
682 delta1=10000,delta2=10000;
685 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
688 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
691 distance2 = TMath::Min(delta1,delta2);
692 if (distance2>100) fDist2 =100;
696 dhelix1.Evaluate(phase[0][0],fXr);
697 dhelix1.GetMomentum(phase[0][0],fPdr);
698 mhelix.GetMomentum(phase[0][1],fPm);
699 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
700 fRr = TMath::Sqrt(radius[0]);
703 dhelix1.Evaluate(phase[1][0],fXr);
704 dhelix1.GetMomentum(phase[1][0], fPdr);
705 mhelix.GetMomentum(phase[1][1], fPm);
706 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
707 fRr = TMath::Sqrt(radius[1]);
709 fDist1 = TMath::Sqrt(distance1);
710 fDist2 = TMath::Sqrt(distance2);
713 Double_t x,alpha,param[5],cov[15];
715 fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
716 fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
717 AliExternalTrackParam paramm(x,alpha,param,cov);
719 fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
720 fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
721 AliExternalTrackParam paramd(x,alpha,param,cov);
726 Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
727 Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
729 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
730 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
731 vnorm2 = TMath::Sqrt(vnorm2);
732 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
733 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
734 pnorm2 = TMath::Sqrt(pnorm2);
736 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
737 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
738 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
743 void AliESDRecKinkInfo::Update()
746 if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
748 // IF BOTH RECONSTRUCTED
749 Float_t distance1,distance2;
750 Double_t xx[3],pp[3];
752 Double_t xd[3],pd[3],signd;
753 Double_t xm[3],pm[3],signm;
754 for (Int_t i=0;i<3;i++){
755 xd[i] = fT2.fTPCinR1[i];
756 pd[i] = fT2.fTPCinP1[i];
757 xm[i] = fT1.fTPCinR1[i];
758 pm[i] = fT1.fTPCinP1[i];
760 signd = fT2.fSign<0 ? -1:1;
761 signm = fT1.fSign<0 ? -1:1;
763 AliHelix dhelix1(xd,pd,signd);
764 dhelix1.GetMomentum(0,pp,0);
765 dhelix1.Evaluate(0,xx);
767 // Double_t x2[3],p2[3];
769 AliHelix mhelix(xm,pm,signm);
771 //find intersection linear
773 Double_t phase[2][2],radius[2];
774 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
775 Double_t delta1=10000,delta2=10000;
778 fMinR = TMath::Sqrt(radius[0]);
781 fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
785 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
786 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
787 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
790 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
791 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
792 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
795 fMinR = TMath::Sqrt(radius[0]);
799 if (radius[0]<radius[1]){
800 fMinR = TMath::Sqrt(radius[0]);
804 fMinR = TMath::Sqrt(radius[1]);
810 distance1 = TMath::Min(delta1,delta2);
812 //find intersection parabolic
814 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
815 delta1=10000,delta2=10000;
818 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
821 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
824 distance2 = TMath::Min(delta1,delta2);
827 dhelix1.Evaluate(phase[0][0],fXr);
828 dhelix1.GetMomentum(phase[0][0],fPdr);
829 mhelix.GetMomentum(phase[0][1],fPm);
830 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
831 fRr = TMath::Sqrt(radius[0]);
834 dhelix1.Evaluate(phase[1][0],fXr);
835 dhelix1.GetMomentum(phase[1][0], fPdr);
836 mhelix.GetMomentum(phase[1][1], fPm);
837 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
838 fRr = TMath::Sqrt(radius[1]);
840 fDist1 = TMath::Sqrt(distance1);
841 fDist2 = TMath::Sqrt(distance2);
844 Double_t x,alpha,param[5],cov[15];
846 fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
847 fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
848 AliExternalTrackParam paramm(x,alpha,param,cov);
850 fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
851 fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
852 AliExternalTrackParam paramd(x,alpha,param,cov);
855 kink.Update(¶mm,¶md);
857 Double_t diff = kink.fRr-fRr;
858 Double_t diff2 = kink.fDist2-fDist2;
859 printf("Diff\t%f\t%f\n",diff,diff2);
870 ////////////////////////////////////////////////////////////////////////
871 AliESDComparisonMI::AliESDComparisonMI()
876 ////////////////////////////////////////////////////////////////////////
877 AliESDComparisonMI::AliESDComparisonMI(const char* fnGenTracks,
879 const char* fnGalice, Int_t direction,
880 Int_t nEvents, Int_t firstEvent)
883 // fFnGenTracks = fnGenTracks;
885 sprintf(fFnGenTracks,"%s",fnGenTracks);
886 sprintf(fFnCmp,"%s",fnCmp);
888 fFirstEventNr = firstEvent;
889 fEventNr = firstEvent;
891 fDirection = direction;
893 fLoader = AliRunLoader::Open(fnGalice);
895 //delete gAlice->GetRunLoader();
899 if (fLoader->LoadgAlice()){
900 cerr<<"Error occured while l"<<endl;
902 Int_t nall = fLoader->GetNumberOfEvents();
910 cerr<<"no events available"<<endl;
914 if (firstEvent+nEvents>nall) {
915 fEventNr = nall-firstEvent;
916 cerr<<"restricted number of events availaible"<<endl;
918 AliMagF * magf = gAlice->Field();
919 AliTracker::SetFieldMap(magf,0);
924 ////////////////////////////////////////////////////////////////////////
925 AliESDComparisonMI::~AliESDComparisonMI()
932 //////////////////////////////////////////////////////////////
933 Int_t AliESDComparisonMI::SetIO()
938 if (!fTreeCmp) return 1;
939 fParamTPC = GetTPCParam();
941 if (!ConnectGenTree()) {
942 cerr<<"Cannot connect tree with generated tracks"<<endl;
948 //////////////////////////////////////////////////////////////
950 Int_t AliESDComparisonMI::SetIO(Int_t eventNr)
956 TFile f("AliESDs.root");
959 TTree* tree = (TTree*) f.Get("esdTree");
962 sprintf(ename,"%d",eventNr);
963 fEvent = (AliESD*)f.Get(ename);
965 sprintf(ename,"ESD%d",eventNr);
966 fEvent = (AliESD*)f.Get(ename);
970 tree->SetBranchStatus("*",1);
971 tree->SetBranchAddress("ESD", &fEvent);
972 tree->SetBranchAddress("ESDfriend.",&fESDfriend);
973 tree->GetEntry(eventNr);
974 fEvent->SetESDfriend(fESDfriend);
980 sprintf(ename,"%d",eventNr);
981 fEvent = (AliESD*)f.Get(ename);
983 sprintf(ename,"ESD%d",eventNr);
984 fEvent = (AliESD*)f.Get(ename);
987 TTree* tree = (TTree*) f.Get("esdTree");
989 Error("CheckESD", "no ESD tree found");
992 tree->SetBranchAddress("ESD", &fEvent);
993 tree->GetEntry(eventNr);
996 if (!fEvent) return 1;
1003 ////////////////////////////////////////////////////////////////////////
1004 void AliESDComparisonMI::Reset()
1011 // fFnCmp = "cmpTracks.root";
1019 ////////////////////////////////////////////////////////////////////////
1020 Int_t AliESDComparisonMI::Exec(Int_t nEvents, Int_t firstEventNr)
1023 fFirstEventNr = firstEventNr;
1027 ////////////////////////////////////////////////////////////////////////
1028 Int_t AliESDComparisonMI::Exec()
1036 fNextTreeGenEntryToRead = 0;
1037 fNextKinkToRead = 0;
1039 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
1040 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
1044 fNParticles = gAlice->GetEvent(fEventNr);
1046 fIndexRecTracks = new Short_t[fNParticles*20]; //write at maximum 4 tracks corresponding to particle
1047 fIndexRecKinks = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle
1048 fIndexRecV0 = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle
1050 fFakeRecTracks = new Short_t[fNParticles];
1051 fMultiRecTracks = new Short_t[fNParticles];
1052 fMultiRecKinks = new Short_t[fNParticles];
1053 fMultiRecV0 = new Short_t[fNParticles];
1055 for (Int_t i = 0; i<fNParticles; i++) {
1056 for (Int_t j=0;j<20;j++){
1057 fIndexRecTracks[20*i+j] = -1;
1058 fIndexRecKinks[20*i+j] = -1;
1059 fIndexRecV0[20*i+j] = -1;
1061 fFakeRecTracks[i] = 0;
1062 fMultiRecTracks[i] = 0;
1063 fMultiRecKinks[i] = 0;
1067 cout<<"Start to process event "<<fEventNr<<endl;
1068 cout<<"\tfNParticles = "<<fNParticles<<endl;
1069 if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
1070 if (TreeTLoop()>0) return 1;
1072 if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
1073 if (TreeGenLoop(eventNr)>0) return 1;
1074 BuildKinkInfo0(eventNr);
1075 //BuildV0Info(eventNr); // no V0 info for a moment
1076 fRecArray->Delete();
1078 if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
1080 delete [] fIndexRecTracks;
1081 delete [] fIndexRecKinks;
1082 delete [] fIndexRecV0;
1083 delete [] fFakeRecTracks;
1084 delete [] fMultiRecTracks;
1085 delete [] fMultiRecKinks;
1086 delete [] fMultiRecV0;
1091 cerr<<"Exec finished"<<endl;
1097 ////////////////////////////////////////////////////////////////////////
1098 Bool_t AliESDComparisonMI::ConnectGenTree()
1101 // connect all branches from the genTracksTree
1102 // use the same variables as for the new cmp tree, it may work
1104 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1105 if (!fFileGenTracks) {
1106 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1109 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1110 if (!fTreeGenTracks) {
1111 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1112 <<fFnGenTracks<<endl;
1116 fMCInfo = new AliMCInfo;
1117 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1120 fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
1121 if (!fTreeGenKinks) {
1122 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1123 <<fFnGenTracks<<endl;
1127 fGenKinkInfo = new AliGenKinkInfo;
1128 fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
1131 fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
1133 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1134 <<fFnGenTracks<<endl;
1138 fGenV0Info = new AliGenV0Info;
1139 fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
1143 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1149 ////////////////////////////////////////////////////////////////////////
1150 void AliESDComparisonMI::CreateTreeCmp()
1152 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1154 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1159 fTreeCmp = new TTree("ESDcmpTracks","ESDcmpTracks");
1160 fMCInfo = new AliMCInfo;
1161 fRecInfo = new AliESDRecInfo;
1162 AliESDtrack * esdTrack = new AliESDtrack;
1163 // AliITStrackMI * itsTrack = new AliITStrackMI;
1164 fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
1165 fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
1166 // fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
1170 fTreeCmpKinks = new TTree("ESDcmpKinks","ESDcmpKinks");
1171 fGenKinkInfo = new AliGenKinkInfo;
1172 fRecKinkInfo = new AliESDRecKinkInfo;
1173 fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
1174 fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
1177 fTreeCmpV0 = new TTree("ESDcmpV0","ESDcmpV0");
1178 fGenV0Info = new AliGenV0Info;
1179 fRecV0Info = new AliESDRecV0Info;
1180 fTreeCmpV0->Branch("MC.","AliGenV0Info", &fGenV0Info,256000);
1181 fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
1183 fTreeCmp->AutoSave();
1184 fTreeCmpKinks->AutoSave();
1185 fTreeCmpV0->AutoSave();
1187 ////////////////////////////////////////////////////////////////////////
1188 void AliESDComparisonMI::CloseOutputFile()
1191 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1202 ////////////////////////////////////////////////////////////////////////
1204 TVector3 AliESDComparisonMI::TR2Local(AliTrackReference *trackRef,
1205 AliTPCParam *paramTPC) {
1207 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1209 paramTPC->Transform0to1(x,index);
1210 paramTPC->Transform1to2Ideal(x,index);
1213 ////////////////////////////////////////////////////////////////////////
1215 Int_t AliESDComparisonMI::TreeTLoop()
1218 // loop over all ESD reconstructed tracks and store info in memory
1220 // + loop over all reconstructed kinks
1224 Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
1225 Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
1226 Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
1227 fSignedKinks = new Short_t[nKinks];
1228 fSignedV0 = new Short_t[nV0MIs];
1230 // load kinks to the memory
1231 for (Int_t i=0; i<nKinks;i++){
1232 AliESDkink * kink =fEvent->GetKink(i);
1236 for (Int_t i=0; i<nV0MIs;i++){
1237 AliV0 * v0MI = (AliV0*)fEvent->GetV0(i);
1243 AliESDtrack * track=0;
1244 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1245 //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
1246 track = (AliESDtrack*)fEvent->GetTrack(iEntry);
1248 Int_t label = track->GetLabel();
1249 Int_t absLabel = abs(label);
1250 if (absLabel < fNParticles) {
1251 // fIndexRecTracks[absLabel] = iEntry;
1252 if (label < 0) fFakeRecTracks[absLabel]++;
1253 if (fMultiRecTracks[absLabel]>0){
1254 if (fMultiRecTracks[absLabel]<20)
1255 fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] = iEntry;
1258 fIndexRecTracks[absLabel*20] = iEntry;
1259 fMultiRecTracks[absLabel]++;
1262 // sort reconstructed kinks
1264 AliESDkink * kink=0;
1265 for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
1266 kink = (AliESDkink*)fEvent->GetKink(iEntry);
1267 if (!kink) continue;
1269 Int_t label0 = TMath::Abs(kink->GetLabel(0));
1270 Int_t label1 = TMath::Abs(kink->GetLabel(1));
1271 Int_t absLabel = TMath::Min(label0,label1);
1272 if (absLabel < fNParticles) {
1273 if (fMultiRecKinks[absLabel]>0){
1274 if (fMultiRecKinks[absLabel]<20)
1275 fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] = iEntry;
1278 fIndexRecKinks[absLabel*20] = iEntry;
1279 fMultiRecKinks[absLabel]++;
1282 // --sort reconstructed V0
1285 for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
1286 v0MI = (AliV0*)fEvent->GetV0(iEntry);
1287 if (!v0MI) continue;
1289 Int_t label0 = TMath::Abs(v0MI->GetLabel(0));
1290 Int_t label1 = TMath::Abs(v0MI->GetLabel(1));
1292 for (Int_t i=0;i<2;i++){
1293 Int_t absLabel = TMath::Abs(v0MI->GetLabel(i));
1294 if (absLabel < fNParticles) {
1295 if (fMultiRecV0[absLabel]>0){
1296 if (fMultiRecV0[absLabel]<20)
1297 fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] = iEntry;
1300 fIndexRecV0[absLabel*20] = iEntry;
1301 fMultiRecV0[absLabel]++;
1307 printf("Time spended in TreeTLoop\n");
1310 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1314 ////////////////////////////////////////////////////////////////////////
1315 Int_t AliESDComparisonMI::TreeGenLoop(Int_t eventNr)
1318 // loop over all entries for a given event, find corresponding
1319 // rec. track and store in the fTreeCmp
1323 Int_t entry = fNextTreeGenEntryToRead;
1324 Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1325 cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1326 <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1327 TBranch * branch = fTreeCmp->GetBranch("RC");
1328 TBranch * branchF = fTreeCmp->GetBranch("F");
1330 branch->SetAddress(&fRecInfo); // set all pointers
1331 fRecArray = new TObjArray(fNParticles);
1332 AliESDtrack dummytrack; //
1333 AliESDfriendTrack dummytrackF; //
1335 while (entry < nParticlesTR) {
1336 fTreeGenTracks->GetEntry(entry);
1338 if (eventNr < fMCInfo->fEventNr) continue;
1339 if (eventNr > fMCInfo->fEventNr) continue;;
1341 fNextTreeGenEntryToRead = entry-1;
1342 if (fDebug > 2 && fMCInfo->fLabel < 10) {
1343 cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1345 // if (fMCInfo->fNTPCRef<1) continue; // not TPCref
1348 AliESDtrack * track=0;
1349 fRecInfo->fReconstructed =0;
1350 TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1351 local.GetXYZ(fRecInfo->fTRLocalCoord);
1353 if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
1354 //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1355 track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
1358 // find nearest track if multifound
1359 //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
1362 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
1363 if ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
1364 if ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
1367 if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
1370 track->GetInnerPxPyPz(p);
1371 Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1373 for (Int_t i=1;i<20;i++){
1374 if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
1375 AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
1376 if (!track2) continue;
1377 //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //
1378 //if (sign2<0) continue;
1379 track2->GetInnerPxPyPz(p);
1380 Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1391 if ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
1392 if ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
1393 if ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
1394 if (status2<status) continue;
1396 if (mom<maxp) continue;
1405 fRecInfo->SetESDtrack(track);
1407 fRecInfo->SetESDtrack(&dummytrack);
1411 fRecInfo->fReconstructed = 1;
1412 fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1413 fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1415 fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent);
1418 fRecInfo->SetESDtrack(&dummytrack);
1419 fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent);
1421 fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
1424 fTreeCmp->AutoSave();
1425 //fTracks->Delete();
1426 printf("Time spended in TreeGenLoop\n");
1428 if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
1435 ////////////////////////////////////////////////////////////////////////
1436 ////////////////////////////////////////////////////////////////////////
1437 ////////////////////////////////////////////////////////////////////////
1438 Int_t AliESDComparisonMI::BuildKinkInfo0(Int_t eventNr)
1441 // loop over all entries for a given event, find corresponding
1442 // rec. track and store in the fTreeCmp
1446 Int_t entry = fNextKinkToRead;
1447 Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
1448 cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
1449 <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
1451 TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
1452 branch->SetAddress(&fRecKinkInfo); // set all pointers
1455 while (entry < nParticlesTR) {
1456 fTreeGenKinks->GetEntry(entry);
1458 if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
1459 if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
1461 fNextKinkToRead = entry-1;
1464 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
1465 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
1466 fRecKinkInfo->fT1 = (*fRecInfo1);
1467 fRecKinkInfo->fT2 = (*fRecInfo2);
1468 fRecKinkInfo->fStatus =0;
1469 if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
1470 if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
1471 if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
1473 if (fRecKinkInfo->fStatus==3){
1474 fRecKinkInfo->Update();
1476 Int_t label = TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1477 Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1480 fRecKinkInfo->fRecStatus =0;
1481 fRecKinkInfo->fMultiple = fMultiRecKinks[label];
1482 fRecKinkInfo->fKinkMultiple=0;
1484 if (fMultiRecKinks[label]>0){
1486 // for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
1487 for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
1488 Int_t index = fIndexRecKinks[label*20+j];
1489 //AliESDkink *kink2 = (AliESDkink*)fKinks->At(index);
1490 AliESDkink *kink2 = (AliESDkink*)fEvent->GetKink(index);
1491 if (TMath::Abs(kink2->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(1))==label2) {
1492 fRecKinkInfo->fKinkMultiple++;
1493 fSignedKinks[index]=1;
1496 // if (kink->fTRDOn) c0++;
1497 //if (kink->fITSOn) c0++;
1498 if (kink->GetStatus(2)>0) c0++;
1499 if (kink->GetStatus(0)>0) c0++;
1502 // if (kink2->fTRDOn) c2++;
1503 //if (kink2->fITSOn) c2++;
1504 if (kink2->GetStatus(2)>0) c2++;
1505 if (kink2->GetStatus(0)>0) c2++;
1507 if (c2<c0) continue;
1510 if (TMath::Abs(kink2->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(0))==label2) {
1511 fRecKinkInfo->fKinkMultiple++;
1512 fSignedKinks[index]=1;
1515 //if (kink->fTRDOn) c0++;
1516 //if (kink->fITSOn) c0++;
1517 if (kink->GetStatus(2)>0) c0++;
1518 if (kink->GetStatus(0)>0) c0++;
1522 // if (kink2->fTRDOn) c2++;
1523 //if (kink2->fITSOn) c2++;
1524 if (kink2->GetStatus(2)>0) c2++;
1525 if (kink2->GetStatus(0)>0) c2++;
1527 if (c2<c0) continue;
1533 fRecKinkInfo->fKink = *kink;
1534 fRecKinkInfo->fRecStatus=1;
1536 fTreeCmpKinks->Fill();
1538 // Int_t nkinks = fKinks->GetEntriesFast();
1539 Int_t nkinks = fEvent->GetNumberOfKinks();
1540 for (Int_t i=0;i<nkinks;i++){
1541 if (fSignedKinks[i]==0){
1542 // AliESDkink *kink = (AliESDkink*)fKinks->At(i);
1543 AliESDkink *kink = (AliESDkink*)fEvent->GetKink(i);
1544 if (!kink) continue;
1546 fRecKinkInfo->fKink = *kink;
1547 fRecKinkInfo->fRecStatus =-2;
1549 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(0)));
1550 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(1)));
1551 if (fRecInfo1 && fRecInfo2){
1552 fRecKinkInfo->fT1 = (*fRecInfo1);
1553 fRecKinkInfo->fT2 = (*fRecInfo2);
1554 fRecKinkInfo->fRecStatus =-1;
1556 fTreeCmpKinks->Fill();
1561 fTreeCmpKinks->AutoSave();
1562 printf("Time spended in BuilKinkInfo Loop\n");
1564 if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1570 ////////////////////////////////////////////////////////////////////////
1571 ////////////////////////////////////////////////////////////////////////
1572 ////////////////////////////////////////////////////////////////////////
1576 Int_t AliESDComparisonMI::BuildV0Info(Int_t eventNr)
1579 // loop over all entries for a given event, find corresponding
1580 // rec. track and store in the fTreeCmp
1584 Int_t entry = fNextV0ToRead;
1585 Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1586 cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1587 <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1589 TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1590 branch->SetAddress(&fRecV0Info); // set all pointers
1591 const AliESDVertex * esdvertex = fEvent->GetVertex();
1592 Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1595 while (entry < nParticlesTR) {
1596 fTreeGenV0->GetEntry(entry);
1598 if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
1599 if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
1601 fNextV0ToRead = entry-1;
1604 AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
1605 AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
1606 if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue; // interactions
1607 if (!fRecInfo1 || !fRecInfo2) continue;
1608 fRecV0Info->fT1 = (*fRecInfo1);
1609 fRecV0Info->fT2 = (*fRecInfo2);
1610 fRecV0Info->fV0Status =0;
1611 if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1612 if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1614 if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1617 if (abs(fRecV0Info->fV0Status)==3){
1618 fRecV0Info->Update(vertex);
1622 Double_t x,alpha, param[5],cov[15];
1623 if ( fRecV0Info->fT1.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){
1624 fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1625 fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
1626 AliExternalTrackParam paramP(x,alpha,param,cov);
1628 fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1629 fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
1630 AliExternalTrackParam paramM(x,alpha,param,cov);
1632 fRecV0Info->fV0tpc->SetParamN(paramM);
1633 fRecV0Info->fV0tpc->SetParamP(paramP);
1634 Double_t pid1[5],pid2[5];
1635 fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1);
1636 fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2);
1638 //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1639 fRecV0Info->fV0tpc->Update(vertex);
1643 fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param);
1644 fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov);
1645 alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha();
1646 new (¶mP) AliExternalTrackParam(x,alpha,param,cov);
1648 fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param);
1649 fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov);
1650 alpha = fRecV0Info->fT2.GetESDtrack()->GetAlpha();
1651 new (¶mM) AliExternalTrackParam(x,alpha,param,cov);
1653 fRecV0Info->fV0its->SetParamN(paramM);
1654 fRecV0Info->fV0its->SetParamP(paramP);
1655 // fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1656 fRecV0Info->fV0its->Update(vertex);
1659 if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
1660 if (fRecV0Info->fDist2>10){
1661 fRecV0Info->Update(vertex);
1663 if (fRecV0Info->fDist2>10){
1664 fRecV0Info->Update(vertex);
1669 // take the V0 from reconstruction
1671 Int_t label = TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1672 Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1674 fRecV0Info->fRecStatus =0;
1675 fRecV0Info->fMultiple = fMultiRecV0[label];
1676 fRecV0Info->fV0Multiple=0;
1678 if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1680 // for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1681 for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1682 Int_t index = fIndexRecV0[label*20+j];
1683 if (index<0) continue;
1684 AliV0 *v0MI2 = (AliV0*)fEvent->GetV0(index);
1685 if (TMath::Abs(v0MI2->GetLabel(0))==label &&TMath::Abs(v0MI2->GetLabel(1))==label2) {
1687 fRecV0Info->fV0Multiple++;
1690 if (TMath::Abs(v0MI2->GetLabel(1))==label &&TMath::Abs(v0MI2->GetLabel(0))==label2) {
1692 fRecV0Info->fV0Multiple++;
1698 fRecV0Info->fV0rec = v0MI;
1699 fRecV0Info->fRecStatus=1;
1707 Int_t nV0MIs = fEvent->GetNumberOfV0s();
1708 for (Int_t i=0;i<nV0MIs;i++){
1709 if (fSignedV0[i]==0){
1710 AliV0 *v0MI = (AliV0*)fEvent->GetV0(i);
1711 if (!v0MI) continue;
1713 fRecV0Info->fV0rec = v0MI;
1714 fRecV0Info->fV0Status =-10;
1715 fRecV0Info->fRecStatus =-2;
1717 // AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->GetLabel(0)));
1718 // AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->GetLabel(1)));
1719 // if (fRecInfo1 && fRecInfo2){
1720 // fRecV0Info->fT1 = (*fRecInfo1);
1721 // fRecV0Info->fT2 = (*fRecInfo2);
1722 // fRecV0Info->fRecStatus =-1;
1724 fRecV0Info->Update(vertex);
1731 fTreeCmpV0->AutoSave();
1732 printf("Time spended in BuilV0Info Loop\n");
1734 if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1737 ////////////////////////////////////////////////////////////////////////
1738 ////////////////////////////////////////////////////////////////////////