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