AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMCTreeTools.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  Class to enable easy correlation of the reconstructed information and MC(true)
18  WARNING : this code is for fast visualization and prototyping NOT FOR PRODUCTION SOFTWARE
19  : Goals was flexibility not speed
20  Example usage:
21  1.) GetParticle properties at vertex (using ESd tree)
22  - e.g to compare reconstructed track momenta and MC momenta ( mentioned in PWGPP-205)
23  esdTree->Draw("Tracks[].fIp.Pt()/AliMCTreeTools::GetValueAt(Entry$,abs(Tracks[].fLabel),0,0,0,8,1,0)","(Tracks[].fTPCncls)>80&&Tracks[].fITSncls>3","")
24  2.) Find TPC space point porperty using closest MC information ( mentioned in ATO-168)
25  - used for loop finder as a ground cluster/tracklet true
26  -- fast as used togeter with MakeCacheTree(TTree * tree, TString varList, TString outFile, TString outTree, TCut selection);
27  author marian.ivanov@cern.ch
28 */
29 
30 
31 /*
32  .L $AliPhysics_SRC/PWGPP/AliMCTreeTools.cxx+
33  AliMCTreeTools::SetWDir(gSystem->pwd());
34  TFile * f = TFile::Open("AliESDs.root");
35  esdTree->Draw("Tracks[].fIp.Pt()/AliMCTreeTools::GetValueAt(Entry$,abs(Tracks[].fLabel),0,0,0,8,1,0)","(Tracks[].fTPCncls)>80&&Tracks[].fITSncls>3","")
36 
37 */
38 
39 
40 #include "TStopwatch.h"
41 #include "TTree.h"
42 #include "TChain.h"
43 #include "TVectorF.h"
44 #include "AliStack.h"
45 #include "TDatabasePDG.h"
46 #include "TParticle.h"
47 #include "TTreeStream.h"
48 #include "AliRunLoader.h"
49 #include "AliTrackReference.h"
50 #include "AliExternalTrackParam.h"
51 #include "AliHelix.h"
52 #include "TCut.h"
53 #include "AliTreePlayer.h"
54 #include "THn.h"
55 #include "TF3.h"
56 #include "TStatToolkit.h"
57 #include <stdarg.h>
58 #include "AliNDLocalRegression.h"
59 #include "AliMCTreeTools.h"
60 
61 ClassImp(AliMCTreeTools)
62 
63 
64 const Double_t kMaxRadius=300;
65 TClonesArray* trackRefs = 0;
66 map<Int_t, TClonesArray*> mapTR;
67 map<Int_t, AliTrackReference*> mapFirstTr;
68 map<Int_t, TClonesArray*> mapHelix;
69 //
70 
72 AliStack *stack=NULL;
73 AliRunLoader *rl=NULL;
78 
79 
80 void AliMCTreeTools::SetWDir(TString dir){
81  wdir=dir;
82 }
83 
84 //
86  // clear cache needed to speedup acces to MC information
87  mapTR.clear();
88  mapHelix.clear();
89  mapFirstTr.clear();
90 }
91 
92 
93 
95  if (rl==NULL){
96  rl = AliRunLoader::Open(Form("%s/galice.root", wdir.Data()));
97  }
98  if (treeTR && rl->GetEventNumber()==iEvent) return;
99  rl->GetEvent(iEvent);
100  rl->LoadTrackRefs();
101  rl->LoadKinematics();
102  //
103  treeTR=rl->TreeTR();
104  stack=rl->Stack();
106 }
107 
108 
109 Double_t AliMCTreeTools::GetValueAt(Int_t iEvent, Int_t itrack, Double_t x, Double_t y, Double_t z, Int_t returnValue, Int_t interpolationType, Int_t verbose){
110  //
111  // GetValueAt
112  // return values
113  // 0 - gx
114  // 1 - gy
115  // 2 - gz
116  // 3 - r
117  // 4 - phi
118  // 5 - deltaPhi (pos-dir)
119  // 6 - minimal distance to track ref.
120  // 7 - weighted interpolated P
121  // 8 - weighted interpolated Pt
122  // 9 - P at first tpc ref
123  // 10 - dEdx at local P
124  // 11 - pdg code
125  // 12 - P at vertex
126  // 13 - P at vertex
127 
128  Int_t index= FindNearestReference(iEvent,itrack,x,y,z,0,verbose);
129  if (returnValue==11 || returnValue==12 || returnValue==13){
130  if (stack==NULL) return -1;
131  TParticle *particle = stack->Particle(itrack);
132  Int_t pdgCode=(particle!=NULL) ? particle->GetPdgCode():-1;
133  if (returnValue==11) return pdgCode;
134  if (returnValue==12) return (particle!=NULL) ? particle->P():-1;
135  if (returnValue==13) return (particle!=NULL) ? particle->Pt():-1;
136  }
137  TClonesArray *helixArray = mapHelix[itrack];
138  TClonesArray *trArray = mapTR[itrack];
139  if (helixArray==NULL) return 0;
140  if (trArray==NULL) return 0;
141  if (helixArray->GetEntriesFast()<2) return 0;
142  AliHelix *helix0=(AliHelix*)helixArray->At(index);
143  AliHelix *helix1=(AliHelix*)helixArray->At(index+1);
144  AliTrackReference *ref0= (AliTrackReference *)trArray->At(index);
145  AliTrackReference *ref1= (AliTrackReference *)trArray->At(index+1);
146  if (ref0->GetTrack()!=itrack){
147  static Int_t wCounter=0;
148  wCounter++;
149  return -1;
150  }
151  Double_t d0=TMath::Sqrt((ref0->X()-x)*(ref0->X()-x)+(ref0->Y()-y)*(ref0->Y()-y)+(ref0->Z()-z)*(ref0->Z()-z));
152  Double_t d1=TMath::Sqrt((ref1->X()-x)*(ref1->X()-x)+(ref1->Y()-y)*(ref1->Y()-y)+(ref1->Z()-z)*(ref1->Z()-z));
153  Double_t w0=d1/(d1+d0);
154  Double_t w1=d0/(d1+d0);
155  //
156  if (helix0==NULL || helix1==NULL) return 0;
157  if (returnValue==6) {
158  Double_t value=TMath::Sqrt(TMath::Min(d0,d1));
159  return value;
160  }
161  if (returnValue==7) {
162  Double_t value =w0*ref0->P()+w1*ref1->P();
163  return value;
164  }
165  if (returnValue==8) {
166  Double_t value =w0*ref0->Pt()+w1*ref1->Pt();
167  return value;
168  }
169  //
170  if (returnValue==9){
171  AliTrackReference * ref = mapFirstTr[itrack];
172  return (ref!=NULL) ? ref->P():0;
173  }
174  if (returnValue==10){
175  TParticle *particle = stack->Particle(itrack);
176  if (particle==NULL) return 0;
177  TParticlePDG *mcParticle = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode());
178  if (mcParticle==NULL) return 0;
179  Double_t mass = mcParticle->Mass();
180  Double_t wP =w0*ref0->P()+w1*ref1->P();
181  Double_t dEdx= AliExternalTrackParam::BetheBlochAleph(wP/mass);
182  return dEdx;
183  }
184  //
185 
186  Double_t xyz[3], dxyz[3], ddxyz[3];
187  if (interpolationType==0){
188  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
189  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
190  if (TMath::Abs(zLoop)>0.5){
191  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
192  xyz[2]-=nLoops*zLoop;
193  }
194  }
195  if (interpolationType==1){
196  helix1->Evaluate(helix1->GetPhase(x,y),xyz,dxyz,ddxyz);
197  Double_t zLoop=TMath::TwoPi()/helix1->GetHelix(4);
198  if (TMath::Abs(zLoop)>0.5){
199  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
200  xyz[2]-=nLoops*zLoop;
201  }
202  }
203  if (interpolationType==2){
204  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
205  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
206  if (TMath::Abs(zLoop)>0.5){
207  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
208  xyz[2]-=nLoops*zLoop;
209  }
210  Double_t xyz1[3], dxyz1[3], ddxyz1[3];
211  helix1->Evaluate(helix1->GetPhase(x,y),xyz1,dxyz1,ddxyz1);
212  zLoop=TMath::TwoPi()/helix1->GetHelix(4);
213  if (TMath::Abs(zLoop)>0.5){
214  Int_t nLoops=TMath::Nint((xyz1[2]-z)/zLoop);
215  xyz1[2]-=nLoops*zLoop;
216  }
217  for (Int_t i=0; i<3;i++){
218  xyz[i]*=w0;
219  xyz[i]+=w1*xyz1[i];
220  }
221  }
222  if (returnValue<3) return xyz[returnValue];
223  if (returnValue==3){
224  return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
225  }
226  if (returnValue==4){
227  return TMath::ATan2(xyz[1],xyz[0]);
228  }
229  if (returnValue==5){
230  return TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(dxyz[1],dxyz[0]);
231  }
232 
233 }
234 
235 
237  //
238  TDatabasePDG *pdg = TDatabasePDG::Instance();
239  if (itrack<0) return 0;
240  InitStack(iEvent);
241  TClonesArray *trefs=mapTR[itrack];
242  if (trefs==NULL){ // cache Trackrefs if not done before
243  treeTR->SetBranchAddress("TrackReferences", &trackRefs);
244  treeTR->GetEntry(stack->TreeKEntry(itrack));
245  mapTR[itrack]=(TClonesArray*)trackRefs->Clone();
246  trefs= mapTR[itrack];
247  TClonesArray * helixArray = new TClonesArray("AliHelix",trackRefs->GetEntries());
248  helixArray->ExpandCreateFast(trackRefs->GetEntries());
249  TParticle *particle = stack->Particle(itrack);
250  TParticlePDG *mcparticle = pdg->GetParticle(particle->GetPdgCode());
251  if (mcparticle==NULL) return 0;
252  //
253  Float_t conversion = -1000/0.299792458/bz; // AliTracker::GetBz();
254  if (mcparticle->Charge()!=0){
255  for (Int_t itr=0; itr<trackRefs->GetEntriesFast(); itr++){
256  AliTrackReference *ref= (AliTrackReference *)trackRefs->At(itr);
257  if (ref->GetTrack()!=itrack) continue;
258  if ( ref->DetectorId()==AliTrackReference::kTPC && mapFirstTr[itrack]==NULL){
259  mapFirstTr[itrack]=(AliTrackReference*)ref->Clone();
260  }
261  Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
262  Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
263  if (ref->P()==0) pxyz[0]+=0.00000000001; // create dummy track reference in case of 0 moment (track disappeared)
264  new ((*helixArray)[itr]) AliHelix(xyz,pxyz,mcparticle->Charge()/3.,conversion);
265  }
266  mapHelix[itrack]=helixArray;
267  }
268  }
269 
270  Int_t nTrackRefs = trefs->GetEntriesFast();
271  Int_t nTPCRef=0;
272  AliTrackReference *refNearest=0;
273  TVectorF fdist(nTrackRefs);
274  for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR) {
275  AliTrackReference* ref = static_cast<AliTrackReference*>(trefs->UncheckedAt(itrR));
276  Double_t lDist=(ref->X()-x)*(ref->X()-x)+(ref->Y()-y)*(ref->Y()-y)+(ref->Z()-z)*(ref->Z()-z);
277  fdist[itrR]=lDist;
278  }
279  Double_t dist=250*250;
280  Int_t index0=0,index1=0;
281  for (Int_t itrR = 1; itrR < nTrackRefs-1; ++itrR){
282  if (fdist[itrR]<dist){
283  dist=fdist[itrR];
284  if (fdist[itrR-1]<fdist[itrR+1]){
285  index0=itrR-1;
286  index1=itrR;
287  }else{
288  index0=itrR;
289  index1=itrR+1;
290  }
291  }
292  refNearest=static_cast<AliTrackReference*>(trefs->UncheckedAt(index0));
293  if (verbose ) {
294  trefs->UncheckedAt(index0)->Print();
295  trefs->UncheckedAt(index1)->Print();
296  }
297  }
298  if (returnValue==0) return index0;
299  if (returnValue==1) return TMath::Sqrt(dist);
300 }
301 
Int_t pdg
map< Int_t, AliTrackReference * > mapFirstTr
double Double_t
Definition: External.C:58
Double_t mass
AliRunLoader * rl
Double_t bz
AliStack * stack
map< Int_t, TClonesArray * > mapTR
static Double_t FindNearestReference(Int_t iEvent, Int_t itrack, Double_t x, Double_t y, Double_t z, Int_t returnValue, Int_t verbose=0)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
static void ClearCache()
map< Int_t, TClonesArray * > mapHelix
TTree * treeCl
static Double_t GetValueAt(Int_t iEvent, Int_t itrack, Double_t x, Double_t y, Double_t z, Int_t returnValue, Int_t interpolationType, Int_t verbose)
static void InitStack(Int_t iEvent)
TTree * treeTR
const Double_t kMaxRadius
TTree * treeMC
TClonesArray * trackRefs
TString wdir
TDirectoryFile * dir