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