AliPhysics  63e47e1 (63e47e1)
AliESDtools.cxx
Go to the documentation of this file.
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  author marian.ivanov@cern.ch
18 */
19 
20 
21 /*
22  TFile * f = TFile::Open("/lustre/nyx/alice/users/marsland/alice-tpc-notes/JIRA/ATO-436/data/LHC15o_pass1_esds/alice/data/2015/LHC15o/000246272/pass1/15000246272021.7701/AliESDs.root");
23  TTree * tree = (TTree*)f->Get("esdTree");
24  // tree = AliXRDPROOFtoolkit::MakeChainRandom("esd.list","esdTree",0,10)
25  .L $AliPhysics_SRC/PWGPP/AliESDtools.cxx+
26  AliESDtools tools;
27  tools.Init(tree);
28  //tree->GetEntry(3);
29  //tree->Scan("Tracks@.GetEntries():SPDVertex.fNContributors:Entry$","SPDVertex.fNContributors>100&&Tracks@.GetEntries()/PrimaryVertex.fNContributors>10")
30  tools->CacheTPCEventInformation()
31 
33  tree->Draw(">>entryList","SPDVertex.fNContributors>100&&Tracks@.GetEntries()/PrimaryVertex.fNContributors>10","entrylist");
34  tree->SetEntryList(entryList)
35  tree->Scan("AliESDtools::GetTrackMatchEff(0,0):AliESDtools::GetTrackCounters(0,0):AliESDtools::GetTrackCounters(4,0):AliESDtools::GetMeanHisTPCVertexA():AliESDtools::GetMeanHisTPCVertexC():Entry$","AliESDtools::SCalculateEventVariables(Entry$)")
36 */
37 
38 
39 #include "TStopwatch.h"
40 #include "TTree.h"
41 #include "TChain.h"
42 #include "TVectorF.h"
43 #include "AliStack.h"
44 #include "TDatabasePDG.h"
45 #include "TParticle.h"
46 #include "TTreeStream.h"
47 #include "AliRunLoader.h"
48 #include "AliTrackReference.h"
49 #include "AliExternalTrackParam.h"
50 #include "AliHelix.h"
51 #include "TCut.h"
52 #include "AliTreePlayer.h"
53 #include "THn.h"
54 #include "TF3.h"
55 #include "TStatToolkit.h"
56 #include <stdarg.h>
57 #include "AliNDLocalRegression.h"
58 #include "AliESDEvent.h"
59 #include "AliESDtools.h"
60 
61 ClassImp(AliESDtools)
62 AliESDtools* AliESDtools::fgInstance;
63 
65  fVerbose(0),
66  fESDtree(NULL),
67  fEvent(NULL),
68  fHisTPCVertexA(nullptr),
69  fHisTPCVertexC(nullptr),
70  fHisTPCVertexACut(nullptr),
71  fHisTPCVertexCCut(nullptr),
72  fHisTPCVertex(nullptr),
73  fHistPhiTPCcounterA(nullptr), // helper histogram phi counteres
74  fHistPhiTPCcounterC(nullptr), // helper histogram for TIdentity tree
75  fHistPhiTPCcounterAITS(nullptr), // helper histogram for TIdentity tree
76  fHistPhiTPCcounterCITS(nullptr), // helper histogram for TIdentity tree
77  fHistPhiITScounterA(nullptr), // helper histogram for TIdentity tree
78  fHistPhiITScounterC(nullptr), // helper histogram for TIdentity tree
79  fCacheTrackCounters(nullptr), // track counter
80  fCacheTrackTPCCountersZ(nullptr), // track counter
81  fCacheTrackdEdxRatio(nullptr), // dedx info counter
82  fCacheTrackNcl(nullptr), // ncl counter
83  fCacheTrackChi2(nullptr), // chi2 counter
84  fCacheTrackMatchEff(nullptr), // matchEff counter
85  fLumiGraph(nullptr) // graph for the interaction rate info for a run
86 {
87  fgInstance=this;
88 }
89 
90 void AliESDtools::Init(TTree *tree) {
91  AliESDtools & tools = *this;
92  if (tools.fESDtree) delete tools.fESDtree;
93  if (!tools.fEvent) tools.fEvent = new AliESDEvent();
94  tools.fESDtree = tree;
95  tools.fEvent->ReadFromTree(tree);
96  if (fHisTPCVertexA == nullptr) {
97  tools.fHisTPCVertexA = new TH1F("hisTPCZA", "hisTPCZA", 1000, -250, 250);
98  tools.fHisTPCVertexC = new TH1F("hisTPCZC", "hisTPCZC", 1000, -250, 250);
99  tools.fHisTPCVertex = new TH1F("hisTPCZ", "hisTPCZ", 1000, -250, 250);
100  tools.fHisTPCVertexACut = new TH1F("hisTPCZACut", "hisTPCZACut", 1000, -250, 250);
101  tools.fHisTPCVertexCCut = new TH1F("hisTPCZCCut", "hisTPCZCCut", 1000, -250, 250);
102  tools.fHisTPCVertex->SetLineColor(1);
103  tools.fHisTPCVertexA->SetLineColor(2);
104  tools.fHisTPCVertexC->SetLineColor(4);
105  tools.fHisTPCVertexACut->SetLineColor(3);
106  tools.fHisTPCVertexCCut->SetLineColor(6);
107  tools.fCacheTrackCounters = new TVectorF(20);
108  tools.fCacheTrackTPCCountersZ = new TVectorF(8);
109  tools.fCacheTrackdEdxRatio = new TVectorF(20);
110  tools.fCacheTrackNcl = new TVectorF(20);
111  tools.fCacheTrackChi2 = new TVectorF(20);
112  tools.fCacheTrackMatchEff = new TVectorF(20);
113  // **************************** Event histograms **************************
114  fHistPhiTPCcounterA = new TH1F("hPhiTPCcounterC", "control histogram to count tracks on the A side in phi ", 36, 0., 18.);
115  fHistPhiTPCcounterC = new TH1F("hPhiTPCcounterA", "control histogram to count tracks on the C side in phi ", 36, 0., 18.);
116  fHistPhiTPCcounterAITS = new TH1F("hPhiTPCcounterAITS", "control histogram to count tracks on the A side in phi ", 36, 0., 18.);
117  fHistPhiTPCcounterCITS = new TH1F("hPhiTPCcounterCITS", "control histogram to count tracks on the C side in phi ", 36, 0., 18.);
118  fHistPhiITScounterA = new TH1F("hPhiITScounterA", "control histogram to count tracks on the A side in phi ", 36, 0., 18.);
119  fHistPhiITScounterC = new TH1F("hPhiITScounterC", "control histogram to count tracks on the C side in phi ", 36, 0., 18.);
120  }
121 }
122 
126  AliESDtools &tools=*this;
127  const Int_t kNCRCut=80;
128  const Double_t kDCACut=5;
129  const Float_t knTrackletCut=1.5;
130  // FILL DCA histograms
131  tools.fHisTPCVertexA->Reset();
132  tools.fHisTPCVertexC->Reset();
133  tools.fHisTPCVertexACut->Reset();
134  tools.fHisTPCVertexCCut->Reset();
135  tools.fHisTPCVertex->Reset();
136  Int_t nTracks=tools.fEvent->GetNumberOfTracks();
137  Int_t selected=0;
138  for (Int_t iTrack=0; iTrack<nTracks; iTrack++){
139  AliESDtrack * pTrack = tools.fEvent->GetTrack(iTrack);
140  Float_t dcaxy,dcaz;
141  if (pTrack== nullptr) continue;
142  if (pTrack->IsOn(AliVTrack::kTPCin)==0) continue;
143  if (pTrack->GetTPCClusterInfo(3,1)<kNCRCut) continue;
144  pTrack->GetImpactParameters(dcaxy,dcaz);
145  if (TMath::Abs(dcaxy)>kDCACut) continue;
146  pTrack->SetESDEvent(fEvent);
147  selected++;
148  if ((pTrack->GetNumberOfTRDClusters()/20.+pTrack->GetNumberOfITSClusters())>knTrackletCut){
149  tools.fHisTPCVertex->Fill(pTrack->GetTPCInnerParam()->GetZ());
150  if (pTrack->GetTgl()>0) tools.fHisTPCVertexACut->Fill(pTrack->GetTPCInnerParam()->GetZ());
151  if (pTrack->GetTgl()<0) tools.fHisTPCVertexCCut->Fill(pTrack->GetTPCInnerParam()->GetZ());
152  }else{
153  if (pTrack->GetTgl()>0) tools.fHisTPCVertexA->Fill(pTrack->GetTPCInnerParam()->GetZ());
154  if (pTrack->GetTgl()<0) tools.fHisTPCVertexC->Fill(pTrack->GetTPCInnerParam()->GetZ());
155 
156  }
157  }
158  /*
159  TPCVertexFit(tools.fHisTPCVertex);
160  TPCVertexFit(tools.fHisTPCVertexA);
161  TPCVertexFit(tools.fHisTPCVertexC);
162  TPCVertexFit(tools.fHisTPCVertexACut);
163  TPCVertexFit(tools.fHisTPCVertexCCut);
164  */
165  if (fVerbose&0x10) printf("%d\n",selected);//acheTPCEventInformation()
166  //
167  return selected;
168 }
169 
170 void AliESDtools::TPCVertexFit(TH1F *hisVertex){
171  //0.5 cm 1 bin
172  // hisVertex=tools.fHisTPCVertexACut;
173  TAxis * axis = hisVertex->GetXaxis();
174  for (Int_t iBin=5; iBin<axis->GetNbins()-5; iBin++){
175  Double_t median10=TMath::Median(10,&(hisVertex->GetArray()[iBin-5]));
176  Double_t rms10=TMath::RMS(10,&(hisVertex->GetArray()[iBin-5]));
177  Double_t val0=TMath::Mean(3,&(hisVertex->GetArray()[iBin-2]));
178  Double_t val1=TMath::Mean(3,&(hisVertex->GetArray()[iBin-1]));
179  Double_t val2=TMath::Mean(3,&(hisVertex->GetArray()[iBin+0]));
180  if (val1>=val0 && val1>=val2 && val1>3+1.5*median10){
181  Double_t xbin=axis->GetBinCenter(iBin);
182  printf("Ibin %d\t%f\t%f\t%f\t%f\n", iBin,xbin,val1,median10,rms10);
183  hisVertex->Fit("gaus","qnrsame+","qnr",xbin-2,xbin+2);
184  }
185  }
186 }
187 
188 
189 
190 //
191 //
193  //AliVEvent *event=InputEvent();
195  fCacheTrackCounters->Zero(); // track counter
196  fCacheTrackCounters->Zero(); // track counter
197  fCacheTrackdEdxRatio->Zero(); // dedx info counter
198  fCacheTrackNcl->Zero();; // ncl counter
199  fCacheTrackChi2->Zero();; // chi2 counter
200  fCacheTrackMatchEff->Zero(); // matchEff counter
201  //
202  //
209  //
210  //
211  const Int_t kNclTPCcut=60;
212  const Int_t kTglCut=1.5;
213  const Int_t kDCACut=5; // 5 cm primary cut
214  const Int_t kMindEdxClustersRegion=15;
215  const Int_t kPtCut=0.100;
216  //
217  //
218  AliTPCdEdxInfo tpcdEdxInfo;
219  for (Int_t itrack=0;itrack<fEvent->GetNumberOfTracks();++itrack)
220  { // Track loop
221 
222  AliESDtrack *track = fEvent->GetTrack(itrack);
223  // AliESDPid *pid = track->GetDetPid();
224  Double_t eta=-100., phiTPC=0.,sectorNumber=0.;
225  Double_t tgl = track->Pz()/track->Pt();
226  Double_t phi = track->Phi()-TMath::Pi();
227  phi = track->GetParameterAtRadius(85,5,7);
228  Double_t sectorNumbertmp = (9*phi/TMath::Pi()+18*(phi<0));
229  if (track == NULL) continue;
230  eta = track->Eta();
231  if (TMath::Abs(eta)>0.9) continue;
232  Bool_t isOnITS = track->IsOn(AliESDtrack::kITSrefit);
233  Bool_t isOnTRD = track->IsOn(AliESDtrack::kTRDrefit);
234  Bool_t isOnTPC = track->IsOn(AliESDtrack::kTPCrefit);
235  //
236  //
237  if (track->GetInnerParam()) {
238  phiTPC = track->GetInnerParam()->GetParameterAtRadius(85,5,7);
239  sectorNumber = (9*phiTPC/TMath::Pi()+18*(phiTPC<0));
240  }
241  //
242  // Count only ITS tracks
243  if ( isOnITS && !isOnTPC ) {
244  if (TMath::Abs(phi)>1e-10){
245  if (tgl>0) fHistPhiITScounterA->Fill(sectorNumbertmp);
246  if (tgl<0) fHistPhiITScounterC->Fill(sectorNumbertmp);
247  }
248  }
249  //
250  if (!track->GetInnerParam()) continue;
251  if (track->IsOn(AliVTrack::kTPCout)==kFALSE) continue;
252  (*fCacheTrackCounters)[4]++; // all TPC track with out flag
253  // TPC track counters with DCAZ
254  for (Int_t izCut=1; izCut<4; izCut++){
255  Float_t impactParam[2];
256  track->GetImpactParameters(impactParam[0],impactParam[1]);
257  if (TMath::Abs(impactParam[0])>kDCACut) continue;
258  if (TMath::Abs(track->GetInnerParam()->GetParameter()[1])<10.*(izCut+1.)) (*fCacheTrackTPCCountersZ)[izCut]++;
259  if (TMath::Abs(impactParam[1])<10.*(izCut+1.)) (*fCacheTrackTPCCountersZ)[izCut+4]++;
260  }
261  //
262  //
263  Float_t dcaRPhi, dcaZ;
264  track->GetImpactParameters(dcaRPhi, dcaZ);
265  Int_t nclTPC = track->GetTPCncls();
266  Int_t nclITS = track->GetITSNcls();
267  Int_t nclTRD = track->GetTRDncls();
268  Int_t nclTOF = track->IsOn(AliVTrack::kTOFout);
269  if (nclTRD<1) nclTRD=-1;
270  if (nclITS<1) nclITS=-1;
271  //if (fNcl<1) fNcl=-1;
272  Double_t chi2TPC = TMath::Sqrt(TMath::Abs(track->GetTPCchi2()/nclTPC));
273  Double_t chi2ITS = TMath::Sqrt(TMath::Abs(track->GetITSchi2()));
274  Double_t chi2TRD = TMath::Sqrt(TMath::Abs(track->GetTRDchi2()));
275  Double_t ptot0 = track->GetP();
276  Double_t qP = track->Charge()/track->P();
277  //
278  //
279  if (nclTPC<kNclTPCcut) continue;
280  if (TMath::Abs(tgl)>kTglCut) continue;
281  if (track->Pt()<kPtCut) continue;
282  if (TMath::Abs(dcaRPhi)>kDCACut || TMath::Abs(dcaZ)>kDCACut) continue;
283  (*fCacheTrackCounters)[5]++;
284  if (TMath::Abs(phiTPC)>1e-10){
285  if (tgl>0) fHistPhiTPCcounterA->Fill(sectorNumber);
286  if (tgl<0) fHistPhiTPCcounterC->Fill(sectorNumber);
287  if(isOnITS){
288  if (tgl>0) fHistPhiTPCcounterAITS->Fill(sectorNumber);
289  if (tgl<0) fHistPhiTPCcounterCITS->Fill(sectorNumber);
290  }
291  }
292  //
293  //
294  Bool_t pileUpCut= ( (nclITS>2) || (nclTRD>40));
295  if (pileUpCut==kFALSE) continue;
296  if (TMath::Min(chi2TPC,100.)<0) continue;
297  (*fCacheTrackCounters)[1]++;
298 
299  //
300  Bool_t itsOK=track->IsOn(AliVTrack::kITSout) && nclITS>2 && chi2ITS>0;
301  Bool_t trdOK=track->IsOn(AliVTrack::kTRDout) && nclTRD>35 && chi2TRD>0;
302  Bool_t tofOK=track->IsOn(AliVTrack::kTOFout);
303  //
304  (*fCacheTrackNcl)[4]+=track->GetTPCncls(0, 63);
305  (*fCacheTrackNcl)[5]+=track->GetTPCncls(64, 127);
306  (*fCacheTrackNcl)[6]+=track->GetTPCncls(128, 159);
307  (*fCacheTrackNcl)[1] += nclTPC;
308  (*fCacheTrackChi2)[1]+= (chi2TPC>0) ? TMath::Sqrt(chi2TPC):2; // sometimes negative chi2?
309 
310  if (itsOK && track->GetTPCdEdxInfo(tpcdEdxInfo)){
311 
312  Bool_t isOK=(tpcdEdxInfo.GetNumberOfCrossedRows(0)>kMindEdxClustersRegion);
313  isOK&=(tpcdEdxInfo.GetNumberOfCrossedRows(1)>kMindEdxClustersRegion);
314  isOK&=(tpcdEdxInfo.GetNumberOfCrossedRows(2)>kMindEdxClustersRegion);
315  isOK&=((tpcdEdxInfo.GetSignalMax(0)>0) && (tpcdEdxInfo.GetSignalMax(1)>0) && (tpcdEdxInfo.GetSignalMax(2)>0));
316  isOK&=((tpcdEdxInfo.GetSignalTot(0)>0) && (tpcdEdxInfo.GetSignalTot(1)>0) && (tpcdEdxInfo.GetSignalTot(2)>0));
317  isOK&=(itsOK||trdOK); // stronger pile-up cut requiring ITS or TRD
318 
319  if (isOK) {
320  (*fCacheTrackCounters)[6]+=1; // Counter with accepted TPC dEdx info
321  (*fCacheTrackdEdxRatio)[0]+=TMath::Log(tpcdEdxInfo.GetSignalMax(3));
322  (*fCacheTrackdEdxRatio)[1]+=TMath::Log(tpcdEdxInfo.GetSignalTot(3));
323  (*fCacheTrackdEdxRatio)[2]+=TMath::Log(tpcdEdxInfo.GetSignalMax(0)/tpcdEdxInfo.GetSignalTot(0));
324  (*fCacheTrackdEdxRatio)[3]+=TMath::Log(tpcdEdxInfo.GetSignalMax(1)/tpcdEdxInfo.GetSignalTot(1));
325  (*fCacheTrackdEdxRatio)[4]+=TMath::Log(tpcdEdxInfo.GetSignalMax(2)/tpcdEdxInfo.GetSignalTot(2));
326  (*fCacheTrackdEdxRatio)[5]+=TMath::Log(tpcdEdxInfo.GetSignalMax(3)/tpcdEdxInfo.GetSignalTot(3));
327  (*fCacheTrackdEdxRatio)[6]+=TMath::Log(tpcdEdxInfo.GetSignalMax(1)/tpcdEdxInfo.GetSignalMax(0));
328  (*fCacheTrackdEdxRatio)[7]+=TMath::Log(tpcdEdxInfo.GetSignalMax(1)/tpcdEdxInfo.GetSignalMax(2));
329  (*fCacheTrackdEdxRatio)[8]+=TMath::Log(tpcdEdxInfo.GetSignalTot(1)/tpcdEdxInfo.GetSignalTot(0));
330  (*fCacheTrackdEdxRatio)[9]+=TMath::Log(tpcdEdxInfo.GetSignalTot(1)/tpcdEdxInfo.GetSignalTot(2));
331  }
332  }
333 
334  if (itsOK) { // ITS
335  (*fCacheTrackCounters)[0]++;
336  (*fCacheTrackNcl)[0] += nclITS;
337  (*fCacheTrackChi2)[0] += TMath::Min(TMath::Sqrt(chi2ITS),10.); // cutoff chi2 10
338  (*fCacheTrackMatchEff)[2]+=trdOK;
339  (*fCacheTrackMatchEff)[3]+=tofOK;
340  (*fCacheTrackChi2)[4]+= (chi2TPC>0) ? TMath::Sqrt(chi2TPC):2; // TPC chi2 in case prolongation to ITS
341  // long tracks properties
342  if (nclITS>4){
343  (*fCacheTrackCounters)[7]++;
344  (*fCacheTrackNcl)[7] += nclITS;
345  (*fCacheTrackChi2)[7]+=TMath::Min(TMath::Sqrt(chi2ITS),10.);
346  }
347  }
348  if (trdOK) {// TRD ///TODO - why chi2TRD could be smaller than 0?
349  (*fCacheTrackCounters)[2]++;
350  (*fCacheTrackNcl)[2] += nclTRD;
351  (*fCacheTrackChi2)[2] += TMath::Sqrt(chi2TRD);
352  (*fCacheTrackMatchEff)[0]+=itsOK;
353  (*fCacheTrackChi2)[5]+= (chi2TPC>0) ? TMath::Sqrt(chi2TPC):2; // TPC chi2 in case prolongation to TRD
354  if (nclTRD>80){
355  (*fCacheTrackCounters)[8]++;
356  (*fCacheTrackNcl)[8] += nclTRD;
357  (*fCacheTrackChi2)[8]+=TMath::Min(TMath::Sqrt(chi2TRD),10.);
358  }
359  }
360  if (tofOK) { // TOF
361  (*fCacheTrackCounters)[3]++;
362  (*fCacheTrackNcl)[3] += 1; // dummy for the moment
363  (*fCacheTrackChi2)[3]+= 1; //
364  }
365  } // end of track LOOP
366 
367  for (Int_t i=0; i<9; i++) if ((*fCacheTrackCounters)[i]>0) (*fCacheTrackNcl)[i]/=(*fCacheTrackCounters)[i];
368  for (Int_t i=0; i<4; i++) if ((*fCacheTrackCounters)[i]>0) (*fCacheTrackChi2)[i]/=(*fCacheTrackCounters)[i];
369 
370  for (Int_t i=4; i<7; i++) if ((*fCacheTrackCounters)[1]>0) (*fCacheTrackNcl)[i]/=(*fCacheTrackCounters)[1];
371  if ((*fCacheTrackCounters)[6]>0){
372  for (Int_t i=0; i<10; i++) (*fCacheTrackdEdxRatio)[i]/=(*fCacheTrackCounters)[6];
373  }
374  //
375  // conditional matching efficiency and chi2
376  if ((*fCacheTrackCounters)[0]>0){
377  (*fCacheTrackMatchEff)[2]/=(*fCacheTrackCounters)[0]; // TRD if ITS
378  (*fCacheTrackMatchEff)[3]/=(*fCacheTrackCounters)[0]; // TOF if ITS
379  (*fCacheTrackChi2)[4]/=(*fCacheTrackCounters)[0];
380  }
381  if ((*fCacheTrackCounters)[2]>0) {
382  (*fCacheTrackMatchEff)[0]/=(*fCacheTrackCounters)[2];
383  (*fCacheTrackChi2)[5]/=(*fCacheTrackCounters)[2];
384  } //ITS if TRD
385  (*fCacheTrackCounters)[9]=fEvent->GetNumberOfTracks(); // original number of ESDtracks
386  return 1;
387 }
388 
389 
400 Int_t AliESDtools::GetNearestTrack(const AliExternalTrackParam * trackMatch, Int_t indexSkip, AliESDEvent*event, Int_t trackType, Int_t paramType, AliExternalTrackParam & paramNearest){
401  //
402  // Find track with closest chi2 distance (assume all track ae propagated to the DCA)
403 
404  // paramType = 0 - global track
405  // 1 - track at inner wall of TPC
406  if (trackMatch==NULL){
407  ::Error("AliAnalysisTaskFilteredTree::GetNearestTrack","invalid track pointer");
408  return -1;
409  }
410  Int_t ntracks=event->GetNumberOfTracks();
411  const Double_t ktglCut=0.1;
412  const Double_t kqptCut=0.4;
413  const Double_t kAlphaCut=0.2;
414  //
415  Double_t chi2Min=100000;
416  Int_t indexMin=-1;
417  for (Int_t itrack=0; itrack<ntracks; itrack++){
418  if (itrack==indexSkip) continue;
419  AliESDtrack *ptrack=event->GetTrack(itrack);
420  if (ptrack==NULL) continue;
421  if (trackType==0 && (ptrack->IsOn(0x1)==kFALSE || ptrack->IsOn(0x10)==kTRUE)) continue; // looks for track without TPC information
422  if (trackType==1 && (ptrack->IsOn(0x10)==kFALSE)) continue; // looks for tracks with TPC information
423  if (trackType==2 && (ptrack->IsOn(0x1)==kFALSE || ptrack->IsOn(0x10)==kFALSE)) continue; // looks for tracks with TPC+ITS information
424 
425  if (ptrack->GetKinkIndex(0)<0) continue; // skip kink daughters
426  const AliExternalTrackParam * track=0; //
427  if (paramType==0) track=ptrack; // Global track
428  if (paramType==1) track=ptrack->GetInnerParam(); // TPC only track at inner wall of TPC
429  if (track==NULL) {
430  continue;
431  }
432  // first rough cuts
433  // fP3 cut
434  if (TMath::Abs((track->GetTgl()-trackMatch->GetTgl()))>ktglCut) continue;
435  // fP4 cut
436  if (TMath::Abs((track->GetSigned1Pt()-trackMatch->GetSigned1Pt()))>kqptCut) continue;
437  // fAlpha cut
438  //Double_t alphaDist=TMath::Abs((track->GetAlpha()-trackMatch->GetAlpha()));
439  Double_t alphaDist=TMath::Abs(TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(trackMatch->Py(),trackMatch->Py()));
440  if (alphaDist>TMath::Pi()) alphaDist-=TMath::TwoPi();
441  if (alphaDist>kAlphaCut) continue;
442  // calculate and extract track with smallest chi2 distance
443  AliExternalTrackParam param(*track);
444  if (param.Rotate(trackMatch->GetAlpha())==kFALSE) continue;
445  if (param.PropagateTo(trackMatch->GetX(),trackMatch->GetBz())==kFALSE) continue;
446  Double_t chi2=trackMatch->GetPredictedChi2(&param);
447  if (chi2<chi2Min){
448  indexMin=itrack;
449  chi2Min=chi2;
450  paramNearest=param;
451  }
452  }
453  return indexMin;
454 
455 }
456 
457 
462 void AliESDtools::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, TTreeStream *pcstream){
463  //
464  // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
465  // marian.ivanov@cern.ch
466  // 0.) Init variables
467  // 1.) GetTrack parameters at TPC inner wall
468  // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
469  //
470  // Logic to be used in reco:
471  // 1.) Find matching ITSalone->TPCalone
472  // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
473  // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
474  const Double_t radiusMatch=84.; // redius to propagate
475  //
476  const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
477  const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
478  const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
479  const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
480  const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
481  const Double_t chi2Cut=100; // chi2 matching cut
482  //
483  if (!esdFriend) return; // not ITS standalone track
484  if (esdFriend->TestSkipBit()) return; // friends tracks not stored
485  Int_t ntracks=esdEvent->GetNumberOfTracks();
486  Float_t bz = esdEvent->GetMagneticField();
487  //
488  // 0.) Get parameters in reference radius TPC Inner wall
489  //
490  //
491  TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
492  TMatrixD vecMomR0(ntracks,6); //
493  TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
494  TMatrixD vecMomR1(ntracks,6); //
495  Double_t xyz[3], pxyz[3]; //
496  for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
497  AliESDtrack *track = esdEvent->GetTrack(iTrack);
498  if(!track) continue;
499  if (track->GetInnerParam()){
500  const AliExternalTrackParam *trackTPC=track->GetInnerParam();
501  trackTPC->GetXYZAt(radiusMatch,bz,xyz);
502  trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
503  for (Int_t i=0; i<3; i++){
504  vecPosR1(iTrack,i)=xyz[i];
505  vecMomR1(iTrack,i)=pxyz[i];
506  }
507  vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
508  vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
509  vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
510  vecMomR1(iTrack,5)= trackTPC->GetTgl();;
511  }
512  AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
513  if(!friendTrack) continue;
514  if (friendTrack->GetITSOut()){
515  const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
516  trackITS->GetXYZAt(radiusMatch,bz,xyz);
517  trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
518  for (Int_t i=0; i<3; i++){
519  vecPosR0(iTrack,i)=xyz[i];
520  vecMomR0(iTrack,i)=pxyz[i];
521  }
522  vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
523  vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
524  vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
525  vecMomR0(iTrack,5)= trackITS->GetTgl();;
526  }
527  }
528  //
529  // 1.) Find closest matching tracks, between the ITS standalone track
530  // and the all other tracks
531  // a.) caltegory - All
532  // b.) category - without ITS
533  //
534  //
535  Int_t ntracksPropagated=0;
536  AliExternalTrackParam extTrackDummy;
537  AliESDtrack esdTrackDummy;
538  AliExternalTrackParam itsAtTPC;
539  AliExternalTrackParam itsAtITSTPC;
540  for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
541  AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
542  if(!track0) continue;
543  if (track0->IsOn(AliVTrack::kTPCin)) continue;
544  AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
545  if (!friendTrack0) continue;
546  //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
547  //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
548  ntracksPropagated++;
549  //
550  // 2.) find clostest TPCtrack
551  // a.) all tracks
552  Double_t minChi2All=10000000;
553  Double_t minChi2TPC=10000000;
554  Double_t minChi2TPCITS=10000000;
555  Int_t indexAll=-1;
556  Int_t indexTPC=-1;
557  Int_t indexTPCITS=-1;
558  Int_t ncandidates0=0; // n candidates - rough cut
559  Int_t ncandidates1=0; // n candidates - rough + chi2 cut
560  itsAtTPC=*(friendTrack0->GetITSOut());
561  itsAtITSTPC=*(friendTrack0->GetITSOut());
562  for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
563  AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
564  if(!track1) continue;
565  if (!track1->IsOn(AliVTrack::kTPCin)) continue;
566  // fast checks
567  //
568  if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
569  if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
570  if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
571  if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
572  if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
573  ncandidates0++;
574  //
575  const AliExternalTrackParam * param1= track1->GetInnerParam();
576  if (!friendTrack0->GetITSOut()) continue;
577  AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
578  if (!outerITS.Rotate(param1->GetAlpha())) continue;
579  if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
580  Double_t chi2 = outerITS.GetPredictedChi2(param1);
581  if (chi2>chi2Cut) continue;
582  ncandidates1++;
583  if (chi2<minChi2All){
584  minChi2All=chi2;
585  indexAll=iTrack1;
586  }
587  if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
588  minChi2TPC=chi2;
589  indexTPC=iTrack1;
590  itsAtTPC=outerITS;
591  }
592  if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
593  minChi2TPCITS=chi2;
594  indexTPCITS=iTrack1;
595  itsAtITSTPC=outerITS;
596  }
597  }
598  //
599  AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
600  AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
601  AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
602  (*pcstream)<<"itsTPC"<<
603  "indexAll="<<indexAll<< // index of closest track (chi2)
604  "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
605  "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
606  "ncandidates0="<<ncandidates0<< // number of candidates
607  "ncandidates1="<<ncandidates1<<
608  //
609  "chi2All="<<minChi2All<< // chi2 of closest tracks
610  "chi2TPC="<<minChi2TPC<<
611  "chi2TPCITS="<<minChi2TPCITS<<
612  //
613  "track0.="<<track0<< // ITS standalone tracks
614  "trackAll.="<<trackAll<< // Closets other track
615  "trackTPC.="<<trackTPC<< // Closest TPC only track
616  "trackTPCITS.="<<trackTPCITS<< // closest combined track
617  //
618  "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
619  "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
620  "\n";
621  }
622 }
623 
624 
625 
TH1F * fHistPhiTPCcounterA
Definition: AliESDtools.h:35
TH1F * fHistPhiITScounterC
Definition: AliESDtools.h:40
double Double_t
Definition: External.C:58
TVectorF * fCacheTrackMatchEff
Definition: AliESDtools.h:46
TVectorF * fCacheTrackTPCCountersZ
Definition: AliESDtools.h:42
TVectorF * fCacheTrackChi2
Definition: AliESDtools.h:45
TH1F * fHistPhiTPCcounterCITS
Definition: AliESDtools.h:38
Double_t bz
Int_t GetNearestTrack(const AliExternalTrackParam *trackMatch, Int_t indexSkip, AliESDEvent *event, Int_t trackType, Int_t paramType, AliExternalTrackParam &paramNearest)
TH1F * fHisTPCVertexC
Definition: AliESDtools.h:31
TH1F * fHisTPCVertexCCut
Definition: AliESDtools.h:34
void Init(TTree *tree)
Definition: AliESDtools.cxx:90
void TPCVertexFit(TH1F *hisVertex)
int Int_t
Definition: External.C:63
TH1F * fHistPhiTPCcounterC
Definition: AliESDtools.h:36
float Float_t
Definition: External.C:68
TH1F * fHistPhiTPCcounterAITS
Definition: AliESDtools.h:37
Int_t CalculateEventVariables()
TH1F * fHisTPCVertexACut
Definition: AliESDtools.h:33
TTree * fESDtree
Definition: AliESDtools.h:28
Int_t fVerbose
Definition: AliESDtools.h:27
TH1F * fHisTPCVertexA
Definition: AliESDtools.h:30
TH1F * fHistPhiITScounterA
Definition: AliESDtools.h:39
void ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, TTreeStream *pcstream)
TVectorF * fCacheTrackNcl
Definition: AliESDtools.h:44
Int_t CacheTPCEventInformation()
caching
TH1F * fHisTPCVertex
Definition: AliESDtools.h:32
AliESDEvent * fEvent
Definition: AliESDtools.h:29
bool Bool_t
Definition: External.C:53
TVectorF * fCacheTrackdEdxRatio
Definition: AliESDtools.h:43
TVectorF * fCacheTrackCounters
Definition: AliESDtools.h:41