]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliESDComparisonMI.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Replacing array of objects by array of pointers
[u/mrichter/AliRoot.git] / STEER / AliESDComparisonMI.C
CommitLineData
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/*
27marian.ivanov@cern.ch
28Usage:
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//
36ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,1,0);
37t2->Exec();
38
39//
40//some cuts definition
51ad6848 41TCut 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 44TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
ae7d73d2 45TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
46
47
f007cb5f 48
49
50
ae7d73d2 51TCut crec("crec","fReconstructed==1");
51ad6848 52TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
f007cb5f 53TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
54
ae7d73d2 55TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
56TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
57TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
58TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
59
60AliESDComparisonDraw comp;
61comp.SetIO();
51ad6848 62TFile f("genHits.root");
63TTree * treel = (TTree*)f.Get("HitLines");
64if (treel) comp->fTree->AddFriend(treel,"L");
ae7d73d2 65
66//
67//example
68comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
51ad6848 69comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
70comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
71comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
72comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
51ad6848 73comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
f007cb5f 74comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
75
ae7d73d2 76
77TH1F his("his","his",100,0,20);
78TH1F hpools("hpools","hpools",100,-7,7);
79TH1F hfake("hfake","hfake",1000,0,150);
80TProfile profp0("profp0","profp0",20,-0.4,0.9)
81
f007cb5f 82comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
ae7d73d2 83comp.fRes->Draw();
84comp.fMean->Draw();
51ad6848 85
f007cb5f 86comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
51ad6848 87comp.fRes->Draw();
88
89comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
90comp.fRes->Draw();
91
92comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5)
93comp.fRes->Draw();
94comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5)
ae7d73d2 95comp.fRes->Draw();
96
97comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0")
98
99TH1F his("his","his",100,0,20);
100TH1F hpools("hpools","hpools",100,-7,7);
101
f007cb5f 102TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
103TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
104TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
105TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
106
107comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
108comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
109comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
110comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
51ad6848 111
51ad6848 112
f007cb5f 113comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1")
114comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1")
115comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1")
116comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1")
51ad6848 117
f007cb5f 118comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)")
119comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)")
120comp.fTree->Draw("fESDTrack.fTPCsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)")
121comp.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 124hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
ae7d73d2 125
126comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
127
f007cb5f 128TProfile 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 192TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
193TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
194TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
195TCut crec("crec","fReconstructed==1");
196TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
197TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
198
199TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
200TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
201TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
202TCut 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
206void 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
236void 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 325void 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
529void 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////
682void 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(&paramm,&paramd);
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////////////////////////////////////////////////////////////////////////
812ESDCmpTr::ESDCmpTr()
813{
814 Reset();
815}
816
817////////////////////////////////////////////////////////////////////////
818ESDCmpTr::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////////////////////////////////////////////////////////////////////////
866ESDCmpTr::~ESDCmpTr()
867{
868 if (fLoader) {
869 delete fLoader;
870 }
871}
872
873//////////////////////////////////////////////////////////////
874Int_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
891Int_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////////////////////////////////////////////////////////////////////////
942void 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////////////////////////////////////////////////////////////////////////
958Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
959{
960 fNEvents = nEvents;
961 fFirstEventNr = firstEventNr;
962 return Exec();
963}
964
965////////////////////////////////////////////////////////////////////////
966Int_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////////////////////////////////////////////////////////////////////////
1036Bool_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////////////////////////////////////////////////////////////////////////
1088void 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////////////////////////////////////////////////////////////////////////
1127void 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
1143TVector3 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
1154Int_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////////////////////////////////////////////////////////////////////////
1257Int_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////////////////////////////////////////////////////////////////////////
1380Int_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
1512void 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 1541Int_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 (&paramP) 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 (&paramM) 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
1705void 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