AliPhysics  a56b849 (a56b849)
 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  if (wdir.Contains(".zip")){
97  rl = AliRunLoader::Open(Form("%s#galice.root", wdir.Data()));
98  }else {
99  rl = AliRunLoader::Open(Form("%s/galice.root", wdir.Data()));
100  }
101  }
102  if (treeTR && rl->GetEventNumber()==iEvent) return;
103  rl->GetEvent(iEvent);
104  rl->LoadTrackRefs();
105  rl->LoadKinematics();
106  //
107  treeTR=rl->TreeTR();
108  stack=rl->Stack();
110 }
111 
122 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,Double_t rotation){
123  //
124  // GetValueAt
125  // return values
126  // 0 - gx - interpolated MC point
127  // 1 - gy -
128  // 2 - gz
129  // 3 - r
130  // 4 - phi
131  // 5 - deltaPhi (pos-dir)
132  // 6 - minimal distance to track ref.
133  // 7 - weighted interpolated P
134  // 8 - weighted interpolated Pt
135  // 9 - P at first tpc ref
136  // 10 - dEdx at local P
137  // 11 - pdg code
138  // 12 - P at vertex
139  // 13 - P at vertex
140  // 14 - lx
141  // 15 - ly
142  // 16 - gx of nearest reference
143  // 17 - gy of nearest reference
144  // 18 - gz of nearest reference
145  // rotate position vector local to global
146  {
147  Double_t rx= x*TMath::Cos(rotation)-y*TMath::Sin(rotation);
148  Double_t ry= x*TMath::Sin(rotation)+y*TMath::Cos(rotation);
149  x=rx;
150  y=ry;
151  }
152 
153  Int_t index= FindNearestReference(iEvent,itrack,x,y,z,0,verbose);
154  if (returnValue==11 || returnValue==12 || returnValue==13){
155  if (stack==NULL) return -1;
156  TParticle *particle = stack->Particle(itrack);
157  Int_t pdgCode=(particle!=NULL) ? particle->GetPdgCode():-1;
158  if (returnValue==11) return pdgCode;
159  if (returnValue==12) return (particle!=NULL) ? particle->P():-1;
160  if (returnValue==13) return (particle!=NULL) ? particle->Pt():-1;
161  }
162  TClonesArray *helixArray = mapHelix[itrack];
163  TClonesArray *trArray = mapTR[itrack];
164  if (helixArray==NULL) return 0;
165  if (trArray==NULL) return 0;
166  if (helixArray->GetEntriesFast()<2) return 0;
167  AliHelix *helix0=(AliHelix*)helixArray->At(index);
168  AliHelix *helix1=(AliHelix*)helixArray->At(index+1);
169  AliTrackReference *ref0= (AliTrackReference *)trArray->At(index);
170  AliTrackReference *ref1= (AliTrackReference *)trArray->At(index+1);
171  if (ref0->GetTrack()!=itrack){
172  static Int_t wCounter=0;
173  wCounter++;
174  return -1;
175  }
176  if (returnValue==16) return ref0->X();
177  if (returnValue==17) return ref0->Y();
178  if (returnValue==18) return ref0->Z();
179 
180  Double_t d0=TMath::Sqrt((ref0->X()-x)*(ref0->X()-x)+(ref0->Y()-y)*(ref0->Y()-y)+(ref0->Z()-z)*(ref0->Z()-z));
181  Double_t d1=TMath::Sqrt((ref1->X()-x)*(ref1->X()-x)+(ref1->Y()-y)*(ref1->Y()-y)+(ref1->Z()-z)*(ref1->Z()-z));
182  Double_t w0=d1/(d1+d0);
183  Double_t w1=d0/(d1+d0);
184  //
185  if (helix0==NULL || helix1==NULL) return 0;
186  if (returnValue==6) {
187  Double_t value=TMath::Sqrt(TMath::Min(d0,d1));
188  return value;
189  }
190  if (returnValue==7) {
191  Double_t value =w0*ref0->P()+w1*ref1->P();
192  return value;
193  }
194  if (returnValue==8) {
195  Double_t value =w0*ref0->Pt()+w1*ref1->Pt();
196  return value;
197  }
198  //
199  if (returnValue==9){
200  AliTrackReference * ref = mapFirstTr[itrack];
201  return (ref!=NULL) ? ref->P():0;
202  }
203  if (returnValue==10){
204  TParticle *particle = stack->Particle(itrack);
205  if (particle==NULL) return 0;
206  TParticlePDG *mcParticle = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode());
207  if (mcParticle==NULL) return 0;
208  Double_t mass = mcParticle->Mass();
209  Double_t wP =w0*ref0->P()+w1*ref1->P();
210  Double_t dEdx= AliExternalTrackParam::BetheBlochAleph(wP/mass);
211  return dEdx;
212  }
213  //
214 
215  Double_t xyz[3], dxyz[3], ddxyz[3];
216  if (interpolationType==0){
217  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
218  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
219  if (TMath::Abs(zLoop)>0.5){
220  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
221  xyz[2]-=nLoops*zLoop;
222  }
223  }
224  if (interpolationType==1){
225  helix1->Evaluate(helix1->GetPhase(x,y),xyz,dxyz,ddxyz);
226  Double_t zLoop=TMath::TwoPi()/helix1->GetHelix(4);
227  if (TMath::Abs(zLoop)>0.5){
228  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
229  xyz[2]-=nLoops*zLoop;
230  }
231  }
232  if (interpolationType==2){
233  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
234  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
235  if (TMath::Abs(zLoop)>0.5){
236  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
237  xyz[2]-=nLoops*zLoop;
238  }
239  Double_t xyz1[3], dxyz1[3], ddxyz1[3];
240  helix1->Evaluate(helix1->GetPhase(x,y),xyz1,dxyz1,ddxyz1);
241  zLoop=TMath::TwoPi()/helix1->GetHelix(4);
242  if (TMath::Abs(zLoop)>0.5){
243  Int_t nLoops=TMath::Nint((xyz1[2]-z)/zLoop);
244  xyz1[2]-=nLoops*zLoop;
245  }
246  for (Int_t i=0; i<3;i++){
247  xyz[i]*=w0;
248  xyz[i]+=w1*xyz1[i];
249  }
250  }
251  if (returnValue<3) return xyz[returnValue];
252  if (returnValue==14 || returnValue==15 ){
253  // global to local
254  Double_t lx= xyz[0]*TMath::Cos(rotation)+xyz[1]*TMath::Sin(rotation);
255  Double_t ly= -xyz[0]*TMath::Sin(rotation)+xyz[1]*TMath::Cos(rotation);
256  if (returnValue==14) return lx;
257  if (returnValue==15) return ly;
258 
259  }
260  if (returnValue==3){
261  return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
262  }
263  if (returnValue==4){
264  return TMath::ATan2(xyz[1],xyz[0]);
265  }
266  if (returnValue==5){
267  return TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(dxyz[1],dxyz[0]);
268  }
269 
270 }
271 
282  Int_t verbose) {
283  //
284  TDatabasePDG *pdg = TDatabasePDG::Instance();
285  if (itrack < 0) return 0;
286  InitStack(iEvent);
287  TClonesArray *trefs = mapTR[itrack];
288  if (trefs == NULL) { // cache Trackrefs if not done before
289  treeTR->SetBranchAddress("TrackReferences", &trackRefs);
290  treeTR->GetEntry(stack->TreeKEntry(itrack));
291  mapTR[itrack] = (TClonesArray *) trackRefs->Clone();
292  trefs = mapTR[itrack];
293  TClonesArray *helixArray = new TClonesArray("AliHelix", trackRefs->GetEntries());
294  helixArray->ExpandCreateFast(trackRefs->GetEntries());
295  TParticle *particle = stack->Particle(itrack);
296  TParticlePDG *mcparticle = pdg->GetParticle(particle->GetPdgCode());
297  if (mcparticle == NULL) return 0;
298  //
299  Float_t conversion = -1000 / 0.299792458 / bz; // AliTracker::GetBz();
300  if (mcparticle->Charge() != 0) {
301  for (Int_t itr = 0; itr < trackRefs->GetEntriesFast(); itr++) {
302  AliTrackReference *ref = (AliTrackReference *) trackRefs->At(itr);
303  if (ref->GetTrack() != itrack) continue;
304  if (ref->DetectorId() == AliTrackReference::kTPC && mapFirstTr[itrack] == NULL) {
305  mapFirstTr[itrack] = (AliTrackReference *) ref->Clone();
306  }
307  Double_t xyz[3] = {ref->X(), ref->Y(), ref->Z()};
308  Double_t pxyz[3] = {ref->Px(), ref->Py(), ref->Pz()};
309  if (ref->P() == 0)
310  pxyz[0] += 0.00000000001; // create dummy track reference in case of 0 moment (track disappeared)
311  new((*helixArray)[itr]) AliHelix(xyz, pxyz, mcparticle->Charge() / 3., conversion);
312  }
313  mapHelix[itrack] = helixArray;
314  }
315  }
316 
317  Int_t nTrackRefs = trefs->GetEntriesFast();
318  Int_t nTPCRef = 0;
319  AliTrackReference *refNearest = 0;
320  TVectorF fdist(nTrackRefs);
321  for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR) {
322  AliTrackReference *ref = static_cast<AliTrackReference *>(trefs->UncheckedAt(itrR));
323  Double_t lDist =
324  (ref->X() - x) * (ref->X() - x) + (ref->Y() - y) * (ref->Y() - y) + (ref->Z() - z) * (ref->Z() - z);
325  fdist[itrR] = lDist;
326  }
327  Double_t dist = 250 * 250;
328  Int_t index0 = 0, index1 = 0;
329  for (Int_t itrR = 1; itrR < nTrackRefs - 1; ++itrR) {
330  if (fdist[itrR] < dist) {
331  dist = fdist[itrR];
332  if (fdist[itrR - 1] < fdist[itrR + 1]) {
333  index0 = itrR - 1;
334  index1 = itrR;
335  } else {
336  index0 = itrR;
337  index1 = itrR + 1;
338  }
339  }
340  refNearest = static_cast<AliTrackReference *>(trefs->UncheckedAt(index0));
341  if (verbose) {
342  trefs->UncheckedAt(index0)->Print();
343  trefs->UncheckedAt(index1)->Print();
344  }
345  }
346  if (returnValue == 0) return index0;
347  if (returnValue == 1) return TMath::Sqrt(dist);
348 }
349 
Int_t pdg
map< Int_t, AliTrackReference * > mapFirstTr
double Double_t
Definition: External.C:58
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, Double_t rotation=0)
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 void InitStack(Int_t iEvent)
TTree * treeTR
const Double_t kMaxRadius
TTree * treeMC
TClonesArray * trackRefs
TString wdir
TDirectoryFile * dir