HLT/trigger/RunD0Trigger.C

   1 //#define __COMPILE__
   2 //#ifdef __COMPILE__
   3 #if !defined(__CINT__) || defined(__MAKECINT__)
   4 //-- --- standard headers-------------
   5 #include <iostream.h>
   6 //--------Root headers ---------------
   7 #include <TSystem.h>
   8 #include <TFile.h>
   9 #include <TString.h>
  10 #include <TStopwatch.h>
  11 #include <TObject.h>
  12 #include <TVector3.h>
  13 #include <TTree.h>
  14 #include <TParticle.h>
  15 #include <TArray.h>
  16 //----- AliRoot headers ---------------
  17 #include "alles.h"
  18 #include "AliRun.h"
  19 #include "AliKalmanTrack.h"
  20 #include "AliITStrackV2.h"
  21 #include "AliHeader.h"
  22 #include "AliGenEventHeader.h"
  23 #include "AliV0vertex.h"
  24 #include "AliV0vertexer.h"
  25 #include "AliITSVertex.h"
  26 #include "AliITSVertexer.h"
  27 #include "AliITSVertexerTracks.h"
  28 #include "AliD0Trigger.h"
  29 #endif
  30 //-------------------------------------
  31 // field (T)
  32 const Double_t kBz = 0.4;
  33
  34 // primary vertex
  35 Double_t primaryvertex[3] = {0.,0.,0,};
  36
  37 // sigle track cuts
  38 const Double_t kPtCut = 0.5;  // GeV/c
  39 const Double_t kd0Cut = 50.; // micron
  40 const Double_t kd0CutHigh = 200.; // micron
  41
  42
  43 //cuts for combined tracks
  44 // cuts[0] = lowest V0 cut
  45 // cuts[1] = highest V0 cut
  46 // cuts[2] = inv. mass cut (diiferense)
  47 // cuts[3] = max value for pointing angle
  48 const Double_t cuts[4] = {0.005,
  49                           0.02,
  50                           0.5,
  51                           0.95};
  52
  53 // this function applies single track cuts
  54 Bool_t TrkCuts(const AliITStrackV2& trk);
  55
  56 // this function creates TObjArrays with positive and negative tracks
  57 void   SelectTracks(TTree& itsTree,
  58                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
  59                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
  60
  61 //void GetPrimaryVertex(int i,Char_t* path="./");
  62
  63 Int_t iTrkP,iTrkN,itsEntries;
  64 Char_t trksName[100];
  65 Int_t nTrksP=0,nTrksN=0;
  66 Int_t nD0=0;
  67 int ev=0;
  68
  69 void   RunD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
  70
  71   const Char_t *name="AliD0Trigger";
  72   cerr<<'\n'<<name<<"...\n";
  73   gBenchmark->Start(name);
  74
  75   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
  76
  77   // Open file with ITS tracks
  78   Char_t fnameTrack[1024];
  79   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  80   TFile* itstrks = TFile::Open(fnameTrack);
  81
  82    // tracks from ITS
  83   sprintf(trksName,"TreeT_ITS_%d",ev);
  84   TTree *itsTree=(TTree*)itstrks->Get(trksName);
  85   if(!itsTree) continue;
  86   itsEntries = (Int_t)itsTree->GetEntries();
  87   printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
  88
  89   //Getting primary Vertex
  90   GetPrimaryVertex(0,path);
  91
  92   // call function which applies sigle track selection and
  93   // separetes positives and negatives
  94   TObjArray trksP(itsEntries/2);
  95   Int_t *itsEntryP = new Int_t[itsEntries];
  96   TObjArray trksN(itsEntries/2);
  97   Int_t *itsEntryN = new Int_t[itsEntries];
  98   SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
  99
 100   cout<<"#pos: "<<nTrksP<<endl;
 101   cout<<"#neg: "<<nTrksN<<endl;
 102
 103   AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
 104
 105   for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 106     postrack = (AliITStrackV2*)trksP.At(iTrkP);
 107     for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 108       negtrack = (AliITStrackV2*)trksN.At(iTrkN);
 109       D0.SetTracks(postrack,negtrack);
 110       if(D0.FindV0()){
 111         D0.FindMomentaAtVertex();
 112         if(D0.FindInvMass()){
 113           if(D0.PointingAngle()){
 114             nD0++;
 115           }
 116         }
 117       }
 118     }
 119   }
 120   cout<<"#D0: "<<nD0<<endl;
 121   gBenchmark->Stop(name);
 122   gBenchmark->Show(name);
 123 }
 124 //___________________________________________________________________________
 125 void   SelectTracks(TTree& itsTree,
 126                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 127                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
 128 //
 129 // this function creates two TObjArrays with positive and negative tracks
 130 //
 131   nTrksP=0,nTrksN=0;
 132
 133
 134   Int_t entr = (Int_t)itsTree.GetEntries();
 135
 136   // trasfer tracks from tree to arrays
 137   for(Int_t i=0; i<entr; i++) {
 138
 139     AliITStrackV2 *itstrack = new AliITStrackV2;
 140     itsTree.SetBranchAddress("tracks",&itstrack);
 141
 142     itsTree.GetEvent(i);
 143
 144     // single track selection
 145     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
 146
 147     if(itstrack->Get1Pt()>0.) { // negative track
 148       trksN.AddLast(itstrack);
 149       itsEntryN[nTrksN] = i;
 150       nTrksN++;
 151     } else {                    // positive track
 152       trksP.AddLast(itstrack);
 153       itsEntryP[nTrksP] = i;
 154       nTrksP++;
 155     }
 156
 157   }
 158
 159   return;
 160 }
 161 //____________________________________________________________________________
 162 Bool_t TrkCuts(const AliITStrackV2& trk) {
 163 //
 164 // this function tells if track passes some kinematical cuts
 165 //
 166   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
 167   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
 168   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
 169
 170   return kTRUE;
 171 }
 172 //____________________________________________________________________________
 173 void GetPrimaryVertex(int i,Char_t* path="./") {
 174
 175   int event=i;
 176
 177   Char_t falice[1024];
 178   sprintf(falice,"%s/galice.root",path);
 179   TFile * galice = new TFile(falice);
 180
 181   TDirectory * curdir;
 182
 183   Char_t vname[20];
 184   galice->cd();
 185
 186   sprintf(vname,"Vertex_%d",event);
 187   TArrayF o = 0;
 188   o.Set(3);
 189   AliHeader * header = 0;
 190
 191   TTree * treeE = (TTree*)gDirectory->Get("TE");
 192   treeE->SetBranchAddress("Header",&header);
 193   treeE->GetEntry(event);
 194   AliGenEventHeader* genHeader = header->GenEventHeader();
 195   if(genHeader){
 196     genHeader->PrimaryVertex(o);
 197     primaryvertex[0] = (Double_t)o[0];
 198     primaryvertex[1] = (Double_t)o[1];
 199     primaryvertex[2] = (Double_t)o[2];
 200   }
 201   else{
 202     printf("Can't find Header");
 203   }
 204   delete header;
 205   delete galice;
 206 }