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