HLT/trigger/RunD0offline.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 //sec. vertex
  38 double v2[3]={0,0,0};
  39
  40 // sigle track cuts
  41 const Double_t kPtCut = 0.5;  // GeV/c
  42 const Double_t kd0Cut = 50.; // micron
  43 const Double_t kd0CutHigh = 200.; // micron
  44
  45
  46 //cuts for combined tracks
  47 const Double_t cuts[7] = {0.005,     // cuts[0] = lowest V0 cut  (cm)
  48                           0.015,     // cuts[1] = highest V0 cut (cm)
  49                           0.05,      // cuts[2] = inv. mass cut (diferense) (Gev/c)
  50                           0.95,      // cuts[3] = max cosine value for pointing angle
  51                           -5000,     // cuts[4] = d0d0
  52                           0.8,       // cuts[5] = costhetastar
  53                           0.5};      // cuts[6] = ptchild
  54 //cut for distance of closest aprach
  55 double cutDCA=0.01;
  56
  57 // this function applies single track cuts
  58 Bool_t TrkCuts(const AliITStrackV2& trk);
  59
  60 // this function creates TObjArrays with positive and negative tracks
  61 void   SelectTracks(TTree& itsTree,
  62                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
  63                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
  64
  65 //void GetPrimaryVertex(int i,Char_t* path="./");
  66
  67 Int_t iTrkP,iTrkN,itsEntries;
  68 Char_t trksName[100];
  69 Int_t nTrksP=0,nTrksN=0;
  70 Int_t nD0=0;
  71 int ev=0;
  72 double mom[6];
  73
  74 void RunD0offline(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
  75
  76   const Char_t *name="AliD0offline";
  77   cerr<<'\n'<<name<<"...\n";
  78   gBenchmark->Start(name);
  79
  80   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
  81
  82   // Open file with ITS tracks
  83   Char_t fnameTrack[1024];
  84   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  85   TFile* itstrks = TFile::Open(fnameTrack);
  86
  87    // tracks from ITS
  88   sprintf(trksName,"TreeT_ITS_%d",ev);
  89   TTree *itsTree=(TTree*)itstrks->Get(trksName);
  90   if(!itsTree) continue;
  91   itsEntries = (Int_t)itsTree->GetEntries();
  92   printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
  93
  94   //Getting primary Vertex
  95   GetPrimaryVertex(0,path);
  96
  97   // call function which applies sigle track selection and
  98   // separetes positives and negatives
  99   TObjArray trksP(itsEntries/2);
 100   Int_t *itsEntryP = new Int_t[itsEntries];
 101   TObjArray trksN(itsEntries/2);
 102   Int_t *itsEntryN = new Int_t[itsEntries];
 103   SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
 104
 105   cout<<"#pos: "<<nTrksP<<endl;
 106   cout<<"#neg: "<<nTrksN<<endl;
 107
 108   //the offline stuff
 109   // define the cuts for vertexing
 110   Double_t vtxcuts[]={33., // max. allowed chi2
 111                       0.0, // min. allowed negative daughter's impact param
 112                       0.0, // min. allowed positive daughter's impact param
 113                       1.0, // max. allowed DCA between the daughter tracks
 114                      -1.0, // min. allowed cosine of V0's pointing angle
 115                       0.0, // min. radius of the fiducial volume
 116                       2.9};// max. radius of the fiducial volume
 117
 118   // create the AliV0vertexer object
 119   AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
 120
 121   AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
 122
 123   double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
 124
 125   for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 126     postrack = (AliITStrackV2*)trksP.At(iTrkP);
 127     for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 128       negtrack = (AliITStrackV2*)trksN.At(iTrkN);
 129       D0.SetTracks(postrack,negtrack);
 130       //
 131       // ----------- DCA MINIMIZATION ------------------
 132       //
 133       // find the DCA and propagate the tracks to the DCA
 134       double dca = vertexer2->PropagateToDCA(negtrack,postrack);
 135
 136       if(dca<cutDCA){
 137         // define the AliV0vertex object
 138         AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
 139         // get position of the vertex
 140         vertex2->GetXYZ(v2[0],v2[1],v2[2]);
 141         delete vertex2;
 142         if(D0.pTchild()){
 143           if(D0.d0d0()){
 144             if(D0.FindV0offline(v2)){
 145
 146               // momenta of the tracks at the vertex
 147               ptP = 1./TMath::Abs(postrack->Get1Pt());
 148               alphaP = postrack->GetAlpha();
 149               phiP = alphaP+TMath::ASin(postrack->GetSnp());
 150               mom[0] = ptP*TMath::Cos(phiP);
 151               mom[1] = ptP*TMath::Sin(phiP);
 152               mom[2] = ptP*postrack->GetTgl();
 153
 154               ptN = 1./TMath::Abs(negtrack->Get1Pt());
 155               alphaN = negtrack->GetAlpha();
 156               phiN = alphaN+TMath::ASin(negtrack->GetSnp());
 157               mom[3] = ptN*TMath::Cos(phiN);
 158               mom[4] = ptN*TMath::Sin(phiN);
 159               mom[5] = ptN*negtrack->GetTgl();
 160
 161               D0.SetMomenta(mom);
 162
 163               if(D0.FindInvMass()){
 164                 if(D0.CosThetaStar()){
 165                   if(D0.PointingAngle()){
 166                     nD0++;
 167                   }
 168                 }
 169               }
 170             }
 171           }
 172         }
 173       }
 174     }
 175   }
 176   cout<<"#D0: "<<nD0<<endl;
 177   gBenchmark->Stop(name);
 178   gBenchmark->Show(name);
 179 }
 180 //___________________________________________________________________________
 181 void   SelectTracks(TTree& itsTree,
 182                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 183                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
 184 //
 185 // this function creates two TObjArrays with positive and negative tracks
 186 //
 187   nTrksP=0,nTrksN=0;
 188
 189
 190   Int_t entr = (Int_t)itsTree.GetEntries();
 191
 192   // trasfer tracks from tree to arrays
 193   for(Int_t i=0; i<entr; i++) {
 194
 195     AliITStrackV2 *itstrack = new AliITStrackV2;
 196     itsTree.SetBranchAddress("tracks",&itstrack);
 197
 198     itsTree.GetEvent(i);
 199
 200     // single track selection
 201     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
 202
 203     if(itstrack->Get1Pt()>0.) { // negative track
 204       trksN.AddLast(itstrack);
 205       itsEntryN[nTrksN] = i;
 206       nTrksN++;
 207     } else {                    // positive track
 208       trksP.AddLast(itstrack);
 209       itsEntryP[nTrksP] = i;
 210       nTrksP++;
 211     }
 212
 213   }
 214
 215   return;
 216 }
 217 //____________________________________________________________________________
 218 Bool_t TrkCuts(const AliITStrackV2& trk) {
 219 //
 220 // this function tells if track passes some kinematical cuts
 221 //
 222   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
 223   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
 224   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
 225
 226   return kTRUE;
 227 }
 228 //____________________________________________________________________________
 229 void GetPrimaryVertex(int i,Char_t* path="./") {
 230
 231   int event=i;
 232
 233   Char_t falice[1024];
 234   sprintf(falice,"%s/galice.root",path);
 235   TFile * galice = new TFile(falice);
 236
 237   TDirectory * curdir;
 238
 239   Char_t vname[20];
 240   galice->cd();
 241
 242   sprintf(vname,"Vertex_%d",event);
 243   TArrayF o = 0;
 244   o.Set(3);
 245   AliHeader * header = 0;
 246
 247   TTree * treeE = (TTree*)gDirectory->Get("TE");
 248   treeE->SetBranchAddress("Header",&header);
 249   treeE->GetEntry(event);
 250   AliGenEventHeader* genHeader = header->GenEventHeader();
 251   if(genHeader){
 252     genHeader->PrimaryVertex(o);
 253     primaryvertex[0] = (Double_t)o[0];
 254     primaryvertex[1] = (Double_t)o[1];
 255     primaryvertex[2] = (Double_t)o[2];
 256   }
 257   else{
 258     printf("Can't find Header");
 259   }
 260   delete header;
 261   delete galice;
 262 }