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 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;        // 0.5 GeV/c
  42 const Double_t kd0Cut = 50.;        // 50  micron
  43 const Double_t kd0CutHigh = 400.;   // 200 micron
  44
  45
  46 //cuts for combined tracks
  47 const Double_t cuts[7] = {0.005,     // 0.005 cuts[0] = lowest V0 cut  (cm)
  48                           0.800,     // 0.015 cuts[1] = highest V0 cut (cm)
  49                           0.012,     // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
  50                           0.8,       // 0.8   cuts[3] = min. cosine value for pointing angle
  51                           -5000,     // -5000 cuts[4] = d0d0
  52                           0,       // 0.8   cuts[5] = costhetastar
  53                           0.5};      // 0.5   cuts[6] = ptchild
  54 //cut for distance of closest aprach
  55 double cutDCA=0.05;   //0.05
  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 //void PtD0(Char_t* path="./");
  68
  69 Int_t iTrkP,iTrkN,itsEntries;
  70 Char_t trksName[100];
  71 Int_t nTrksP=0,nTrksN=0;
  72 Int_t nD0=0;
  73 int ev=0;
  74 double mom[6];
  75 int event[10000];
  76 int index=0;
  77
  78 void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
  79   Char_t falice[1024];
  80   sprintf(falice,"%s/galice.root",path);
  81   TFile *f = new TFile(falice);
  82   gAlice=(AliRun*)f->Get("gAlice");
  83   int nEvent=gAlice->GetEventsPerRun();
  84   //int nEvent=20;
  85   delete gAlice;
  86
  87   const Char_t *name="AliD0offline";
  88   cerr<<'\n'<<name<<"...\n";
  89   cout<<"#Events: "<<nEvent<<endl;
  90   gBenchmark->Start(name);
  91
  92   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
  93
  94   // Open file with ITS tracks
  95   Char_t fnameTrack[1024];
  96   //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  97   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  98   TFile* itstrks = TFile::Open(fnameTrack);
  99
 100   //the offline stuff
 101   // define the cuts for vertexing
 102   Double_t vtxcuts[]={33., // max. allowed chi2
 103                       0.0, // min. allowed negative daughter's impact param
 104                       0.0, // min. allowed positive daughter's impact param
 105                       1.0, // max. allowed DCA between the daughter tracks
 106                       -1.0, // min. allowed cosine of V0's pointing angle
 107                       0.0, // min. radius of the fiducial volume
 108                       2.9};// max. radius of the fiducial volume
 109
 110   TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
 111   TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
 112   TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
 113   TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
 114
 115   for(ev=0;ev<nEvent;ev++) {
 116
 117     cout<<"\n Event: "<<ev<<endl;
 118
 119     // tracks from ITS
 120     sprintf(trksName,"TreeT_ITS_%d",ev);
 121     TTree *itsTree=(TTree*)itstrks->Get(trksName);
 122     if(!itsTree) continue;
 123     itsEntries = (Int_t)itsTree->GetEntries();
 124     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
 125
 126     //Getting primary Vertex
 127     GetPrimaryVertex(ev,path);
 128
 129     // call function which applies sigle track selection and
 130     // separetes positives and negatives
 131     TObjArray trksP(itsEntries/2);
 132     Int_t *itsEntryP = new Int_t[itsEntries];
 133     TObjArray trksN(itsEntries/2);
 134     Int_t *itsEntryN = new Int_t[itsEntries];
 135     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
 136
 137     cout<<"#pos: "<<nTrksP<<endl;
 138     cout<<"#neg: "<<nTrksN<<endl;
 139
 140     // create the AliV0vertexer object
 141     AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
 142
 143     AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
 144
 145     double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
 146
 147     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 148       postrack = (AliITStrackV2*)trksP.At(iTrkP);
 149       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 150         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
 151         D0->SetTracks(postrack,negtrack);
 152         //
 153         // ----------- DCA MINIMIZATION ------------------
 154         //
 155         // find the DCA and propagate the tracks to the DCA
 156         double dca = vertexer2->PropagateToDCA(negtrack,postrack);
 157
 158         if(dca<cutDCA){
 159           // define the AliV0vertex object
 160           AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
 161           // get position of the vertex
 162           vertex2->GetXYZ(v2[0],v2[1],v2[2]);
 163           delete vertex2;
 164           //if(D0->FindV0offline(v2) && D0->d0d0()){
 165           if(D0->d0d0()){
 166             D0->SetV0(v2);
 167             D0->FindMomentaOffline();
 168             if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
 169               nD0++;
 170               event[index]=ev; index++;
 171               if(h){
 172                 h1->Fill(D0->Pt());
 173                 h3->Fill(D0->Eta());
 174               }
 175             }
 176           }
 177         }
 178       }
 179     }
 180   }
 181   cout<<"\n#D0: "<<nD0<<endl;
 182   for(int i=0;i<nD0;i++)
 183     cout<<event[i]<<endl;
 184   gBenchmark->Stop(name);
 185   gBenchmark->Show(name);
 186
 187   if(h){
 188     if(!PtD0){
 189       TCanvas *c = new TCanvas("c","c",700,1000);
 190       c->Divide(1,2);
 191       c->cd(1);
 192       h1->Draw();
 193       pt = new TPaveText(7.3,0.4,11,1.8,"br");
 194       pt->SetFillColor(0);
 195       pt->SetBorderSize(1);
 196       pt->AddText("Cuts:");
 197       Char_t st[1024];
 198       sprintf(st,"First Pt:         %g",kPtCut);
 199       pt->AddText(st);
 200       sprintf(st,"d0 low:           %g",kd0Cut);
 201       pt->AddText(st);
 202       sprintf(st,"d0 high:          %g",kd0CutHigh);
 203       pt->AddText(st);
 204       sprintf(st,"V0 low:           %g",cuts[0]);
 205       pt->AddText(st);
 206       sprintf(st,"V0 high:          %g",cuts[1]);
 207       pt->AddText(st);
 208       sprintf(st,"InvMass Diff:     %g",cuts[2]);
 209       pt->AddText(st);
 210       sprintf(st,"cosPointingAngle: %g",cuts[3]);
 211       pt->AddText(st);
 212       sprintf(st,"d0d0:             %g",cuts[4]);
 213       pt->AddText(st);
 214       sprintf(st,"cosThetaStar:     %g",cuts[5]);
 215       pt->AddText(st);
 216       sprintf(st,"PtChild:          %g",cuts[6]);
 217       pt->AddText(st);
 218       sprintf(st,"DCA:              %g",cutDCA);
 219       pt->AddText(st);
 220       pt->Draw();
 221       c_1->Modified();
 222       c->cd();
 223       c->cd(2);
 224       h3->Draw();
 225     }
 226     if(PtD0){
 227       PtD0(path,h2,h4);
 228       TCanvas *c = new TCanvas("c","c",1000,700);
 229       c->Divide(2,2);
 230       c->cd(1);
 231       h1->Draw();
 232       pt = new TPaveText(7.3,0.4,11,1.8,"br");
 233       pt->SetFillColor(0);
 234       pt->SetBorderSize(1);
 235       pt->AddText("Cuts:");
 236       Char_t st[1024];
 237       sprintf(st,"First Pt:         %g",kPtCut);
 238       pt->AddText(st);
 239       sprintf(st,"d0 low:           %g",kd0Cut);
 240       pt->AddText(st);
 241       sprintf(st,"d0 high:          %g",kd0CutHigh);
 242       pt->AddText(st);
 243       sprintf(st,"V0 low:           %g",cuts[0]);
 244       pt->AddText(st);
 245       sprintf(st,"V0 high:          %g",cuts[1]);
 246       pt->AddText(st);
 247       sprintf(st,"InvMass Diff:     %g",cuts[2]);
 248       pt->AddText(st);
 249       sprintf(st,"cosPointingAngle: %g",cuts[3]);
 250       pt->AddText(st);
 251       sprintf(st,"d0d0:             %g",cuts[4]);
 252       pt->AddText(st);
 253       sprintf(st,"cosThetaStar:     %g",cuts[5]);
 254       pt->AddText(st);
 255       sprintf(st,"PtChild:          %g",cuts[6]);
 256       pt->AddText(st);
 257       sprintf(st,"DCA:              %g",cutDCA);
 258       pt->AddText(st);
 259       pt->Draw();
 260       c_1->Modified();
 261       c->cd();
 262       c->cd(2);
 263       h3->Draw();
 264       c->cd(3);
 265       h2->Draw();
 266       c->cd(4);
 267       h4->Draw();
 268     }
 269     //TFile *outfile = new TFile("results.root","recreate");
 270     //h1->Write();
 271     //outfile->Close();
 272
 273   }
 274 }
 275 //___________________________________________________________________________
 276 void   SelectTracks(TTree& itsTree,
 277                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 278                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
 279 //
 280 // this function creates two TObjArrays with positive and negative tracks
 281 //
 282   nTrksP=0,nTrksN=0;
 283
 284
 285   Int_t entr = (Int_t)itsTree.GetEntries();
 286
 287   // trasfer tracks from tree to arrays
 288   for(Int_t i=0; i<entr; i++) {
 289
 290     AliITStrackV2 *itstrack = new AliITStrackV2;
 291     itsTree.SetBranchAddress("tracks",&itstrack);
 292
 293     itsTree.GetEvent(i);
 294
 295     // single track selection
 296     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
 297
 298     if(itstrack->Get1Pt()>0.) { // negative track
 299       trksN.AddLast(itstrack);
 300       itsEntryN[nTrksN] = i;
 301       nTrksN++;
 302     } else {                    // positive track
 303       trksP.AddLast(itstrack);
 304       itsEntryP[nTrksP] = i;
 305       nTrksP++;
 306     }
 307
 308   }
 309
 310   return;
 311 }
 312 //____________________________________________________________________________
 313 Bool_t TrkCuts(const AliITStrackV2& trk) {
 314 //
 315 // this function tells if track passes some kinematical cuts
 316 //
 317   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
 318   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
 319   //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
 320
 321   return kTRUE;
 322 }
 323 //____________________________________________________________________________
 324 void GetPrimaryVertex(int i,Char_t* path="./") {
 325
 326   int event=i;
 327
 328   Char_t falice[1024];
 329   sprintf(falice,"%s/galice.root",path);
 330   TFile * galice = new TFile(falice);
 331
 332   TDirectory * curdir;
 333
 334   Char_t vname[20];
 335   galice->cd();
 336
 337   sprintf(vname,"Vertex_%d",event);
 338   TArrayF o = 0;
 339   o.Set(3);
 340   AliHeader * header = 0;
 341
 342   TTree * treeE = (TTree*)gDirectory->Get("TE");
 343   treeE->SetBranchAddress("Header",&header);
 344   treeE->GetEntry(event);
 345   AliGenEventHeader* genHeader = header->GenEventHeader();
 346   if(genHeader){
 347     genHeader->PrimaryVertex(o);
 348     primaryvertex[0] = (Double_t)o[0];
 349     primaryvertex[1] = (Double_t)o[1];
 350     primaryvertex[2] = (Double_t)o[2];
 351   }
 352   else{
 353     printf("Can't find Header");
 354   }
 355   delete header;
 356   delete galice;
 357 }
 358 //____________________________________________________________________________
 359 PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
 360
 361   Char_t falice[1024];
 362   sprintf(falice,"%s/galice.root",path);
 363   TFile *f = new TFile(falice);
 364   gAlice=(AliRun*)f->Get("gAlice");
 365
 366   TParticle *p;
 367   int nd0=0;
 368   int nkminus =0;
 369   int npipluss = 0;
 370   int nEvent=gAlice->GetEventsPerRun();
 371
 372   for (Int_t i = 0; i <nEvent; i++) {
 373     cout<<"Event: "<<i<<endl;
 374     gAlice->GetEvent(i);
 375     Int_t nPart = gAlice->GetNtrack();
 376     for (Int_t iPart = 0; iPart < nPart; iPart++) {
 377       //cout<<"Particlenr.: "<<iPart<<endl;
 378       p = (TParticle*)gAlice->Particle(iPart);
 379       if (p->GetPdgCode()==421){
 380         if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
 381         h4->Fill(p.Eta());
 382         nd0++;
 383       }
 384       if (p->GetPdgCode()==-321){
 385         nkminus++;
 386       }
 387       if (p->GetPdgCode()==211){
 388         npipluss++;
 389       }
 390     }
 391   }
 392   delete gAlice;
 393 }