AliPhysics  ff07904 (ff07904)
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  AliMCTreeTools::SetWDir(gSystem->pwd());
33  TFile * f = TFile::Open("AliESDs.root");
34  esdTree->Draw("Tracks[].fIp.Pt()/AliMCTreeTools::GetValueAt(Entry$,abs(Tracks[].fLabel),0,0,0,8,1,0)","(Tracks[].fTPCncls)>80&&Tracks[].fITSncls>3","")
35  //
36  AliMCTreeTools::MakeChains();
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 std::map<Int_t, TClonesArray*> mapTR;
67 std::map<Int_t, AliTrackReference*> mapFirstTr;
68 std::map<Int_t, TClonesArray*> mapHelix;
69 //
70 
72 AliStack *stack=NULL;
73 AliRunLoader *rl=NULL;
78 TChain *AliMCTreeTools::fKineChain=0;
79 TChain *AliMCTreeTools::fTRChain=0;
80 
81 void AliMCTreeTools::SetWDir(TString dir){
82  wdir=dir;
83 }
84 
85 //
87  // clear cache needed to speedup acces to MC information
88  mapTR.clear();
89  mapHelix.clear();
90  mapFirstTr.clear();
91 }
92 
93 
94 
96  if (rl==NULL){
97  if (wdir.Contains(".zip")){
98  rl = AliRunLoader::Open(Form("%s#galice.root", wdir.Data()));
99  }else {
100  rl = AliRunLoader::Open(Form("%s/galice.root", wdir.Data()));
101  }
102  }
103  if (treeTR && rl->GetEventNumber()==iEvent) return;
104  rl->GetEvent(iEvent);
105  rl->LoadTrackRefs();
106  rl->LoadKinematics();
107  //
108  treeTR=rl->TreeTR();
109  stack=rl->Stack();
111 }
112 
123 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){
124  //
125  // GetValueAt
126  // return values
127  // 0 - gx - interpolated MC point
128  // 1 - gy -
129  // 2 - gz
130  // 3 - r
131  // 4 - phi
132  // 5 - deltaPhi (pos-dir)
133  // 6 - minimal distance to track ref.
134  // 7 - weighted interpolated P
135  // 8 - weighted interpolated Pt
136  // 9 - P at first tpc ref
137  // 10 - dEdx at local P
138  // 11 - pdg code
139  // 12 - P at vertex
140  // 13 - P at vertex
141  // 14 - lx
142  // 15 - ly
143  // 16 - gx of nearest reference
144  // 17 - gy of nearest reference
145  // 18 - gz of nearest reference
146  // rotate position vector local to global
147  {
148  Double_t rx= x*TMath::Cos(rotation)-y*TMath::Sin(rotation);
149  Double_t ry= x*TMath::Sin(rotation)+y*TMath::Cos(rotation);
150  x=rx;
151  y=ry;
152  }
153 
154  Int_t index= FindNearestReference(iEvent,itrack,x,y,z,0,verbose);
155  if (returnValue==11 || returnValue==12 || returnValue==13){
156  if (stack==NULL) return -1;
157  TParticle *particle = stack->Particle(itrack);
158  Int_t pdgCode=(particle!=NULL) ? particle->GetPdgCode():-1;
159  if (returnValue==11) return pdgCode;
160  if (returnValue==12) return (particle!=NULL) ? particle->P():-1;
161  if (returnValue==13) return (particle!=NULL) ? particle->Pt():-1;
162  }
163  TClonesArray *helixArray = mapHelix[itrack];
164  TClonesArray *trArray = mapTR[itrack];
165  if (helixArray==NULL) return 0;
166  if (trArray==NULL) return 0;
167  if (helixArray->GetEntriesFast()<2) return 0;
168  AliHelix *helix0=(AliHelix*)helixArray->At(index);
169  AliHelix *helix1=(AliHelix*)helixArray->At(index+1);
170  AliTrackReference *ref0= (AliTrackReference *)trArray->At(index);
171  AliTrackReference *ref1= (AliTrackReference *)trArray->At(index+1);
172  if (ref0->GetTrack()!=itrack){
173  static Int_t wCounter=0;
174  wCounter++;
175  return -1;
176  }
177  if (returnValue==16) return ref0->X();
178  if (returnValue==17) return ref0->Y();
179  if (returnValue==18) return ref0->Z();
180 
181  Double_t d0=TMath::Sqrt((ref0->X()-x)*(ref0->X()-x)+(ref0->Y()-y)*(ref0->Y()-y)+(ref0->Z()-z)*(ref0->Z()-z));
182  Double_t d1=TMath::Sqrt((ref1->X()-x)*(ref1->X()-x)+(ref1->Y()-y)*(ref1->Y()-y)+(ref1->Z()-z)*(ref1->Z()-z));
183  Double_t w0=d1/(d1+d0);
184  Double_t w1=d0/(d1+d0);
185  //
186  if (helix0==NULL || helix1==NULL) return 0;
187  if (returnValue==6) {
188  Double_t value=TMath::Sqrt(TMath::Min(d0,d1));
189  return value;
190  }
191  if (returnValue==7) {
192  Double_t value =w0*ref0->P()+w1*ref1->P();
193  return value;
194  }
195  if (returnValue==8) {
196  Double_t value =w0*ref0->Pt()+w1*ref1->Pt();
197  return value;
198  }
199  //
200  if (returnValue==9){
201  AliTrackReference * ref = mapFirstTr[itrack];
202  return (ref!=NULL) ? ref->P():0;
203  }
204  if (returnValue==10){
205  TParticle *particle = stack->Particle(itrack);
206  if (particle==NULL) return 0;
207  TParticlePDG *mcParticle = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode());
208  if (mcParticle==NULL) return 0;
209  Double_t mass = mcParticle->Mass();
210  Double_t wP =w0*ref0->P()+w1*ref1->P();
211  Double_t dEdx= AliExternalTrackParam::BetheBlochAleph(wP/mass);
212  return dEdx;
213  }
214  //
215 
216  Double_t xyz[3], dxyz[3], ddxyz[3];
217  if (interpolationType==0){
218  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
219  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
220  if (TMath::Abs(zLoop)>0.5){
221  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
222  xyz[2]-=nLoops*zLoop;
223  }
224  }
225  if (interpolationType==1){
226  helix1->Evaluate(helix1->GetPhase(x,y),xyz,dxyz,ddxyz);
227  Double_t zLoop=TMath::TwoPi()/helix1->GetHelix(4);
228  if (TMath::Abs(zLoop)>0.5){
229  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
230  xyz[2]-=nLoops*zLoop;
231  }
232  }
233  if (interpolationType==2){
234  helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
235  Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
236  if (TMath::Abs(zLoop)>0.5){
237  Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
238  xyz[2]-=nLoops*zLoop;
239  }
240  Double_t xyz1[3], dxyz1[3], ddxyz1[3];
241  helix1->Evaluate(helix1->GetPhase(x,y),xyz1,dxyz1,ddxyz1);
242  zLoop=TMath::TwoPi()/helix1->GetHelix(4);
243  if (TMath::Abs(zLoop)>0.5){
244  Int_t nLoops=TMath::Nint((xyz1[2]-z)/zLoop);
245  xyz1[2]-=nLoops*zLoop;
246  }
247  for (Int_t i=0; i<3;i++){
248  xyz[i]*=w0;
249  xyz[i]+=w1*xyz1[i];
250  }
251  }
252  if (returnValue<3) return xyz[returnValue];
253  if (returnValue==14 || returnValue==15 ){
254  // global to local
255  Double_t lx= xyz[0]*TMath::Cos(rotation)+xyz[1]*TMath::Sin(rotation);
256  Double_t ly= -xyz[0]*TMath::Sin(rotation)+xyz[1]*TMath::Cos(rotation);
257  if (returnValue==14) return lx;
258  if (returnValue==15) return ly;
259 
260  }
261  if (returnValue==3){
262  return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
263  }
264  if (returnValue==4){
265  return TMath::ATan2(xyz[1],xyz[0]);
266  }
267  if (returnValue==5){
268  return TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(dxyz[1],dxyz[0]);
269  }
270 
271 }
272 
283  Int_t verbose) {
284  //
285  TDatabasePDG *pdg = TDatabasePDG::Instance();
286  if (itrack < 0) return 0;
287  InitStack(iEvent);
288  TClonesArray *trefs = mapTR[itrack];
289  if (trefs == NULL) { // cache Trackrefs if not done before
290  treeTR->SetBranchAddress("TrackReferences", &trackRefs);
291  treeTR->GetEntry(stack->TreeKEntry(itrack));
292  mapTR[itrack] = (TClonesArray *) trackRefs->Clone();
293  trefs = mapTR[itrack];
294  TClonesArray *helixArray = new TClonesArray("AliHelix", trackRefs->GetEntries());
295  helixArray->ExpandCreateFast(trackRefs->GetEntries());
296  TParticle *particle = stack->Particle(itrack);
297  TParticlePDG *mcparticle = pdg->GetParticle(particle->GetPdgCode());
298  if (mcparticle == NULL) return 0;
299  //
300  Float_t conversion = -1000 / 0.299792458 / bz; // AliTracker::GetBz();
301  if (mcparticle->Charge() != 0) {
302  for (Int_t itr = 0; itr < trackRefs->GetEntriesFast(); itr++) {
303  AliTrackReference *ref = (AliTrackReference *) trackRefs->At(itr);
304  if (ref->GetTrack() != itrack) continue;
305  if (ref->DetectorId() == AliTrackReference::kTPC && mapFirstTr[itrack] == NULL) {
306  mapFirstTr[itrack] = (AliTrackReference *) ref->Clone();
307  }
308  Double_t xyz[3] = {ref->X(), ref->Y(), ref->Z()};
309  Double_t pxyz[3] = {ref->Px(), ref->Py(), ref->Pz()};
310  if (ref->P() == 0)
311  pxyz[0] += 0.00000000001; // create dummy track reference in case of 0 moment (track disappeared)
312  new((*helixArray)[itr]) AliHelix(xyz, pxyz, mcparticle->Charge() / 3., conversion);
313  }
314  mapHelix[itrack] = helixArray;
315  }
316  }
317 
318  Int_t nTrackRefs = trefs->GetEntriesFast();
319  Int_t nTPCRef = 0;
320  AliTrackReference *refNearest = 0;
321  TVectorF fdist(nTrackRefs);
322  for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR) {
323  AliTrackReference *ref = static_cast<AliTrackReference *>(trefs->UncheckedAt(itrR));
324  Double_t lDist =
325  (ref->X() - x) * (ref->X() - x) + (ref->Y() - y) * (ref->Y() - y) + (ref->Z() - z) * (ref->Z() - z);
326  fdist[itrR] = lDist;
327  }
328  Double_t dist = 250 * 250;
329  Int_t index0 = 0, index1 = 0;
330  for (Int_t itrR = 1; itrR < nTrackRefs - 1; ++itrR) {
331  if (fdist[itrR] < dist) {
332  dist = fdist[itrR];
333  if (fdist[itrR - 1] < fdist[itrR + 1]) {
334  index0 = itrR - 1;
335  index1 = itrR;
336  } else {
337  index0 = itrR;
338  index1 = itrR + 1;
339  }
340  }
341  refNearest = static_cast<AliTrackReference *>(trefs->UncheckedAt(index0));
342  if (verbose) {
343  trefs->UncheckedAt(index0)->Print();
344  trefs->UncheckedAt(index1)->Print();
345  }
346  }
347  if (returnValue == 0) return index0;
348  if (returnValue == 1) return TMath::Sqrt(dist);
349 }
350 
354 TChain * AliMCTreeTools::MakeKineChain(const char *inputList){
355  TChain * fKineChain = new TChain("KineChain","KineChain");
356  TChain * fTRChain = new TChain("TRChain","TRChain");
357  TObjArray *kineArray = NULL;
358  if (inputList) {
359  kineArray = gSystem->GetFromPipe("cat kinematics.list").Tokenize("\n");
360  }else{
361  kineArray=new TObjArray();
362  kineArray->AddLast(new TObjString("Kinematics.root"));
363  }
364  Int_t nFiles=kineArray->GetEntries();
365  for (Int_t iFile=0; iFile<nFiles; iFile++){
366  TFile * fkine = TFile::Open(kineArray->At(iFile)->GetName());
367  if (fkine==NULL) continue;
368  TList * kineList = fkine->GetListOfKeys();
369  for (Int_t iKey=0; iKey<kineList->GetEntries();iKey++){
370  fKineChain->AddFile(kineArray->At(iFile)->GetName(), TChain::kBigNumber, TString::Format("%s/TreeK",kineList->At(iKey)->GetName()).Data());
371  }
372  TString refName=kineArray->At(iFile)->GetName();
373  refName.ReplaceAll("Kinematics.root", "TrackRefs.root");
374  TFile * ftr = TFile::Open(refName.Data());
375  TList * trList = ftr->GetListOfKeys();
376  for (Int_t iKey=0; iKey<trList->GetEntries();iKey++){
377  fTRChain->AddFile(refName.Data(), TChain::kBigNumber, TString::Format("%s/TreeTR",trList->At(iKey)->GetName()).Data());
378  }
379  }
380  fKineChain->AddFriend(fTRChain,"TR");
381  fTRChain->AddFriend(fKineChain,"K.");
382  return fKineChain;
383 }
Int_t pdg
std::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)
static TChain * fKineChain
Definition: AliMCTreeTools.h:9
Double_t mass
AliRunLoader * rl
TSystem * gSystem
Double_t bz
std::map< Int_t, TClonesArray * > mapTR
AliStack * stack
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)
static TChain * fTRChain
int Int_t
Definition: External.C:63
std::map< Int_t, TClonesArray * > mapHelix
float Float_t
Definition: External.C:68
static void ClearCache()
TTree * treeCl
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
static void InitStack(Int_t iEvent)
static TChain * MakeKineChain(const char *inputList="kinematics.list")
TTree * treeTR
const Double_t kMaxRadius
TTree * treeMC
TClonesArray * trackRefs
TString wdir
TDirectoryFile * dir