AliMCTreeTools.cxx

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
  Class to enable easy correlation of the reconstructed information and MC(true)    
  WARNING : this code is for fast visualization and prototyping NOT FOR PRODUCTION SOFTWARE
          : Goals was flexibility not speed
  Example usage:
    1.) GetParticle properties at vertex (using ESd tree) 
        - e.g to compare reconstructed track momenta and MC momenta  ( mentioned in PWGPP-205)     
        esdTree->Draw("Tracks[].fIp.Pt()/AliMCTreeTools::GetValueAt(Entry$,abs(Tracks[].fLabel),0,0,0,8,1,0)","(Tracks[].fTPCncls)>80&&Tracks[].fITSncls>3","") 
    2.) Find TPC space point porperty using closest MC information  ( mentioned in ATO-168) 
        - used for loop finder as a ground cluster/tracklet true
  -- fast as used  togeter with MakeCacheTree(TTree * tree, TString varList, TString outFile, TString outTree, TCut selection);          
   author marian.ivanov@cern.ch
*/


/*
  .L $AliPhysics_SRC/PWGPP/AliMCTreeTools.cxx+
  AliMCTreeTools::SetWDir(gSystem->pwd());
  TFile * f = TFile::Open("AliESDs.root");
  esdTree->Draw("Tracks[].fIp.Pt()/AliMCTreeTools::GetValueAt(Entry$,abs(Tracks[].fLabel),0,0,0,8,1,0)","(Tracks[].fTPCncls)>80&&Tracks[].fITSncls>3","")

*/


#include "TStopwatch.h"
#include "TTree.h" 
#include "TChain.h"
#include "TVectorF.h"
#include "AliStack.h"
#include "TDatabasePDG.h"
#include "TParticle.h"
#include "TTreeStream.h"
#include "AliRunLoader.h"
#include "AliTrackReference.h"
#include "AliExternalTrackParam.h"
#include "AliHelix.h"
#include "TCut.h"
#include "AliTreePlayer.h"
#include "THn.h"
#include "TF3.h"
#include "TStatToolkit.h"
#include <stdarg.h>
#include "AliNDLocalRegression.h"
#include "AliMCTreeTools.h"

ClassImp(AliMCTreeTools)


const Double_t kMaxRadius=300; 
TClonesArray* trackRefs = 0;
map<Int_t, TClonesArray*> mapTR;
map<Int_t, AliTrackReference*> mapFirstTr;
map<Int_t, TClonesArray*> mapHelix;
//

TString wdir="";
AliStack *stack=NULL;
AliRunLoader *rl=NULL;
TTree *treeTR=0;
TTree * treeCl=0;
TTree * treeMC=0;
Double_t bz=-5;


void AliMCTreeTools::SetWDir(TString dir){
  wdir=dir;
}

//
void AliMCTreeTools::ClearCache(){
  // clear cache needed to speedup acces to MC information
  mapTR.clear();
  mapHelix.clear(); 
  mapFirstTr.clear();
}



void  AliMCTreeTools::InitStack(Int_t iEvent){
  if (rl==NULL){
    if (wdir.Contains(".zip")){
      rl = AliRunLoader::Open(Form("%s#galice.root", wdir.Data()));
    }else {
      rl = AliRunLoader::Open(Form("%s/galice.root", wdir.Data()));
    }
  }
  if (treeTR && rl->GetEventNumber()==iEvent) return; 
  rl->GetEvent(iEvent);
  rl->LoadTrackRefs();
  rl->LoadKinematics();
  //
  treeTR=rl->TreeTR();
  stack=rl->Stack();
  AliMCTreeTools::ClearCache();
}

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){
  //
  // GetValueAt
  // return values 
  //    0 - gx    - interpolated MC point
  //    1 - gy    -
  //    2 - gz
  //    3 - r
  //    4 - phi
  //    5 - deltaPhi (pos-dir)
  //    6 - minimal distance to track ref.
  //    7 - weighted interpolated  P
  //    8 - weighted interpolated  Pt
  //    9  - P at first tpc ref
  //    10 - dEdx at local P
  //    11 - pdg code
  //    12 - P at vertex
  //    13 - P at vertex
  //    14 - lx
  //    15 - ly
  //    16 - gx of nearest reference
  //    17 - gy of nearest reference
  //    18 - gz of nearest reference
  // rotate position vector local to global
  {
    Double_t rx=  x*TMath::Cos(rotation)-y*TMath::Sin(rotation);
    Double_t ry=  x*TMath::Sin(rotation)+y*TMath::Cos(rotation);
    x=rx;
    y=ry;
  }

  Int_t index= FindNearestReference(iEvent,itrack,x,y,z,0,verbose);
  if (returnValue==11 || returnValue==12 ||  returnValue==13){
    if (stack==NULL) return -1;
    TParticle *particle = stack->Particle(itrack);
    Int_t pdgCode=(particle!=NULL) ? particle->GetPdgCode():-1;
    if (returnValue==11) return pdgCode;
    if (returnValue==12) return (particle!=NULL) ? particle->P():-1;
    if (returnValue==13) return (particle!=NULL) ? particle->Pt():-1;
  }
  TClonesArray *helixArray = mapHelix[itrack];
  TClonesArray *trArray = mapTR[itrack];
  if (helixArray==NULL) return 0;
  if (trArray==NULL) return 0;
  if (helixArray->GetEntriesFast()<2) return 0;
  AliHelix *helix0=(AliHelix*)helixArray->At(index);
  AliHelix *helix1=(AliHelix*)helixArray->At(index+1);
  AliTrackReference *ref0= (AliTrackReference *)trArray->At(index);
  AliTrackReference *ref1= (AliTrackReference *)trArray->At(index+1);
  if (ref0->GetTrack()!=itrack){
    static Int_t wCounter=0;
    wCounter++;
    return -1;
  }
  if (returnValue==16) return ref0->X();
  if (returnValue==17) return ref0->Y();
  if (returnValue==18) return ref0->Z();

  Double_t d0=TMath::Sqrt((ref0->X()-x)*(ref0->X()-x)+(ref0->Y()-y)*(ref0->Y()-y)+(ref0->Z()-z)*(ref0->Z()-z));
  Double_t d1=TMath::Sqrt((ref1->X()-x)*(ref1->X()-x)+(ref1->Y()-y)*(ref1->Y()-y)+(ref1->Z()-z)*(ref1->Z()-z));
  Double_t w0=d1/(d1+d0);
  Double_t w1=d0/(d1+d0);
  //
  if (helix0==NULL || helix1==NULL) return 0;
  if (returnValue==6) {
    Double_t value=TMath::Sqrt(TMath::Min(d0,d1));
    return value;
  }
  if (returnValue==7) {
    Double_t value =w0*ref0->P()+w1*ref1->P();
    return value;
  }
  if (returnValue==8) {
    Double_t value =w0*ref0->Pt()+w1*ref1->Pt();
    return value;
  }
  //
  if (returnValue==9){
    AliTrackReference * ref =  mapFirstTr[itrack];
    return (ref!=NULL) ? ref->P():0;
  }
   if (returnValue==10){
     TParticle *particle = stack->Particle(itrack);
     if (particle==NULL) return 0;
     TParticlePDG *mcParticle = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode());
     if (mcParticle==NULL) return 0;
     Double_t mass = mcParticle->Mass();
     Double_t wP =w0*ref0->P()+w1*ref1->P();
     Double_t dEdx= AliExternalTrackParam::BetheBlochAleph(wP/mass);
     return dEdx;
  }
  //

  Double_t xyz[3],  dxyz[3], ddxyz[3]; 
  if (interpolationType==0){
    helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
    Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4);
    if (TMath::Abs(zLoop)>0.5){
      Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
      xyz[2]-=nLoops*zLoop;
    }
  }
  if (interpolationType==1){
    helix1->Evaluate(helix1->GetPhase(x,y),xyz,dxyz,ddxyz);
    Double_t zLoop=TMath::TwoPi()/helix1->GetHelix(4);
    if (TMath::Abs(zLoop)>0.5){
      Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
      xyz[2]-=nLoops*zLoop;
    }
  }
  if (interpolationType==2){
    helix0->Evaluate(helix0->GetPhase(x,y),xyz,dxyz,ddxyz);
    Double_t zLoop=TMath::TwoPi()/helix0->GetHelix(4); 
    if (TMath::Abs(zLoop)>0.5){
      Int_t nLoops=TMath::Nint((xyz[2]-z)/zLoop);
      xyz[2]-=nLoops*zLoop;
    }
    Double_t xyz1[3],  dxyz1[3], ddxyz1[3]; 
    helix1->Evaluate(helix1->GetPhase(x,y),xyz1,dxyz1,ddxyz1);
    zLoop=TMath::TwoPi()/helix1->GetHelix(4); 
    if (TMath::Abs(zLoop)>0.5){
      Int_t nLoops=TMath::Nint((xyz1[2]-z)/zLoop);
      xyz1[2]-=nLoops*zLoop;
    }
    for (Int_t i=0; i<3;i++){
      xyz[i]*=w0;
      xyz[i]+=w1*xyz1[i];
    }
  }
  if (returnValue<3) return xyz[returnValue];
  if (returnValue==14 || returnValue==15 ){
    // global to local
    Double_t lx=  xyz[0]*TMath::Cos(rotation)+xyz[1]*TMath::Sin(rotation);
    Double_t ly=  -xyz[0]*TMath::Sin(rotation)+xyz[1]*TMath::Cos(rotation);
    if (returnValue==14) return lx;
    if (returnValue==15) return ly;

  }
  if (returnValue==3){
    return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  }
  if (returnValue==4){
    return TMath::ATan2(xyz[1],xyz[0]);
  }
  if (returnValue==5){
    return TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(dxyz[1],dxyz[0]);
  }

}

Double_t AliMCTreeTools::FindNearestReference(Int_t iEvent, Int_t itrack, Double_t x, Double_t y, Double_t z, Int_t returnValue,
                                     Int_t verbose) {
  //  
  TDatabasePDG *pdg = TDatabasePDG::Instance();
  if (itrack < 0) return 0;
  InitStack(iEvent);
  TClonesArray *trefs = mapTR[itrack];
  if (trefs == NULL) { // cache Trackrefs if not done before
    treeTR->SetBranchAddress("TrackReferences", &trackRefs);
    treeTR->GetEntry(stack->TreeKEntry(itrack));
    mapTR[itrack] = (TClonesArray *) trackRefs->Clone();
    trefs = mapTR[itrack];
    TClonesArray *helixArray = new TClonesArray("AliHelix", trackRefs->GetEntries());
    helixArray->ExpandCreateFast(trackRefs->GetEntries());
    TParticle *particle = stack->Particle(itrack);
    TParticlePDG *mcparticle = pdg->GetParticle(particle->GetPdgCode());
    if (mcparticle == NULL) return 0;
    //
    Float_t conversion = -1000 / 0.299792458 / bz; // AliTracker::GetBz();
    if (mcparticle->Charge() != 0) {
      for (Int_t itr = 0; itr < trackRefs->GetEntriesFast(); itr++) {
        AliTrackReference *ref = (AliTrackReference *) trackRefs->At(itr);
        if (ref->GetTrack() != itrack) continue;
        if (ref->DetectorId() == AliTrackReference::kTPC && mapFirstTr[itrack] == NULL) {
          mapFirstTr[itrack] = (AliTrackReference *) ref->Clone();
        }
        Double_t xyz[3] = {ref->X(), ref->Y(), ref->Z()};
        Double_t pxyz[3] = {ref->Px(), ref->Py(), ref->Pz()};
        if (ref->P() == 0)
          pxyz[0] += 0.00000000001; // create dummy track reference in case of 0 moment (track disappeared)
        new((*helixArray)[itr]) AliHelix(xyz, pxyz, mcparticle->Charge() / 3., conversion);
      }
      mapHelix[itrack] = helixArray;
    }
  }

  Int_t nTrackRefs = trefs->GetEntriesFast();
  Int_t nTPCRef = 0;
  AliTrackReference *refNearest = 0;
  TVectorF fdist(nTrackRefs);
  for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR) {
    AliTrackReference *ref = static_cast<AliTrackReference *>(trefs->UncheckedAt(itrR));
    Double_t lDist =
            (ref->X() - x) * (ref->X() - x) + (ref->Y() - y) * (ref->Y() - y) + (ref->Z() - z) * (ref->Z() - z);
    fdist[itrR] = lDist;
  }
  Double_t dist = 250 * 250;
  Int_t index0 = 0, index1 = 0;
  for (Int_t itrR = 1; itrR < nTrackRefs - 1; ++itrR) {
    if (fdist[itrR] < dist) {
      dist = fdist[itrR];
      if (fdist[itrR - 1] < fdist[itrR + 1]) {
        index0 = itrR - 1;
        index1 = itrR;
      } else {
        index0 = itrR;
        index1 = itrR + 1;
      }
    }
    refNearest = static_cast<AliTrackReference *>(trefs->UncheckedAt(index0));
    if (verbose) {
      trefs->UncheckedAt(index0)->Print();
      trefs->UncheckedAt(index1)->Print();
    }
  }
  if (returnValue == 0) return index0;
  if (returnValue == 1) return TMath::Sqrt(dist);
}