HLT/trigger/RunD0offline.C

   1 //#define __COMPILE__
   2 //#ifdef __COMPILE__
   3 #if !defined(__CINT__) || defined(__MAKECINT__)
   4 //-- --- standard headers-------------
   5 #include <iostream.h>
   6 #include <fstream.h>
   7 //--------Root headers ---------------
   8 #include <TSystem.h>
   9 #include <TFile.h>
  10 #include <TString.h>
  11 #include <TStopwatch.h>
  12 #include <TObject.h>
  13 #include <TVector3.h>
  14 #include <TTree.h>
  15 #include <TObjArray>
  16 #include <TParticle.h>
  17 #include <TArray.h>
  18 //----- AliRoot headers ---------------
  19 #include "alles.h"
  20 #include "AliRun.h"
  21 #include "AliKalmanTrack.h"
  22 #include "AliITStrackV2.h"
  23 #include "AliHeader.h"
  24 #include "AliGenEventHeader.h"
  25 #include "AliV0vertex.h"
  26 #include "AliV0vertexer.h"
  27 #include "AliITSVertex.h"
  28 #include "AliITSVertexer.h"
  29 #include "AliITSVertexerTracks.h"
  30 #include "AliD0Trigger.h"
  31 #endif
  32 //-------------------------------------
  33 typedef struct {
  34   Int_t lab;
  35   Int_t pdg;
  36   Int_t mumlab;
  37   Int_t mumpdg;
  38   Float_t Vx,Vy,Vz;
  39   Float_t Px,Py,Pz;
  40 } RECTRACK;
  41
  42 // field (T)
  43 const Double_t kBz = 0.4;
  44
  45 // primary vertex
  46 double primaryvertex[3]={0.,0.,0.};
  47
  48 //sec. vertex
  49 double v2[3]={0,0,0};
  50
  51 // sigle track cuts
  52 const Double_t kPtCut = 0.5;        // 0.5 GeV/c
  53 const Double_t kd0Cut = 50.;        // 50  micron
  54 const Double_t kd0CutHigh = 400.;   // 200 micron
  55
  56
  57 //cuts for combined tracks
  58 const Double_t cuts[7] = {0.005,     // 0.005 cuts[0] = lowest V0 cut  (cm)
  59                           0.800,     // 0.015 cuts[1] = highest V0 cut (cm)
  60                           0.012,     // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
  61                           0.8,       // 0.8   cuts[3] = min. cosine value for pointing angle
  62                           -5000,     // -5000 cuts[4] = d0d0
  63                           0,         // 0.8   cuts[5] = costhetastar
  64                           0.5};      // 0.5   cuts[6] = ptchild
  65 //cut for distance of closest aprach
  66 double cutDCA=0.05;   //0.05
  67
  68 // this function applies single track cuts
  69 Bool_t TrkCuts(const AliITStrackV2& trk);
  70
  71 // this function creates TObjArrays with positive and negative tracks
  72 void   SelectTracks(TTree& itsTree,
  73                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
  74                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
  75
  76 //void GetPrimaryVertex(int i,Char_t* path="./");
  77
  78 //void PtD0(Char_t* path="./");
  79
  80 Int_t iTrkP,iTrkN,itsEntries;
  81 Char_t trksName[100],refsName[100];
  82 Int_t nTrksP=0,nTrksN=0;
  83 Int_t nD0=0;
  84 int ev=0;
  85 double mom[6];
  86 int event[10000];
  87 int index=0;
  88 Bool_t isSignal;
  89 Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
  90 Int_t pdg[2],mum[2],mumlab[2];
  91 RECTRACK rectrk;
  92
  93 void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
  94
  95   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
  96
  97   // Open file with ITS tracks
  98   Char_t fnameTrack[1024];
  99   //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
 100   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
 101   TFile* itstrks = TFile::Open(fnameTrack);
 102
 103   Char_t refFile[1024];
 104   //sprintf(fnameTrack,"%s/ITStracksRefFile.root",path);
 105   sprintf(refFile,"%s/ITStracksRefFile.root",path);
 106   TFile* itsrefs = TFile::Open(refFile);
 107
 108
 109   //the offline stuff
 110   // define the cuts for vertexing
 111   Double_t vtxcuts[]={33., // max. allowed chi2
 112                       0.0, // min. allowed negative daughter's impact param
 113                       0.0, // min. allowed positive daughter's impact param
 114                       1.0, // max. allowed DCA between the daughter tracks
 115                       -1.0, // min. allowed cosine of V0's pointing angle
 116                       0.0, // min. radius of the fiducial volume
 117                       2.9};// max. radius of the fiducial volume
 118
 119   TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
 120   TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
 121   TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
 122   TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
 123
 124   Char_t falice[1024];
 125   sprintf(falice,"%s/galice.root",path);
 126   TFile *f = new TFile(falice);
 127   gAlice=(AliRun*)f->Get("gAlice");
 128   int nEvent=gAlice->GetEventsPerRun();
 129   //int nEvent=20;
 130   cout<<"#Events: "<<nEvent<<endl;
 131   delete gAlice;
 132
 133   const Char_t *name="AliD0offline";
 134   cerr<<'\n'<<name<<"...\n";
 135   gBenchmark->Start(name);
 136
 137   for(ev=0;ev<nEvent;ev++) {
 138
 139     cout<<"\n Event: "<<ev<<endl;
 140
 141     // tracks from ITS
 142     sprintf(trksName,"TreeT_ITS_%d",ev);
 143     TTree *itsTree=(TTree*)itstrks->Get(trksName);
 144     if(!itsTree) continue;
 145     itsEntries = (Int_t)itsTree->GetEntries();
 146     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
 147
 148     // tree from reference file
 149
 150     sprintf(refsName,"Tree_Ref_%d",ev);
 151     TTree *refTree=(TTree*)itsrefs->Get(refsName);
 152     refTree->SetBranchAddress("rectracks",&rectrk);
 153
 154     //Getting primary Vertex
 155     GetPrimaryVertex(ev,path);
 156
 157     // count the total number of events
 158     nTotEv++;
 159
 160     // call function which applies sigle track selection and
 161     // separetes positives and negatives
 162     TObjArray trksP(itsEntries/2);
 163     Int_t *itsEntryP = new Int_t[itsEntries];
 164     TObjArray trksN(itsEntries/2);
 165     Int_t *itsEntryN = new Int_t[itsEntries];
 166     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
 167
 168     cout<<"#pos: "<<nTrksP<<endl;
 169     cout<<"#neg: "<<nTrksN<<endl;
 170
 171     // create the AliV0vertexer object
 172     AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
 173
 174     AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
 175
 176     double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
 177
 178     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
 179       postrack = (AliITStrackV2*)trksP.At(iTrkP);
 180
 181       // get info on tracks PDG and mothers PDG from reference file
 182       refTree->GetEvent(itsEntryP[iTrkP]);
 183       pdg[0] = rectrk.pdg;
 184       mum[0] = rectrk.mumpdg;
 185       mumlab[0] = rectrk.mumlab;
 186
 187       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
 188         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
 189
 190         // get info on tracks PDG and mothers PDG from reference file
 191         refTree->GetEvent(itsEntryN[iTrkN]);
 192         pdg[1] = rectrk.pdg;
 193         mum[1] = rectrk.mumpdg;
 194         mumlab[1] = rectrk.mumlab;
 195
 196         D0->SetTracks(postrack,negtrack);
 197
 198         // ----------- DCA MINIMIZATION ------------------
 199         //
 200         // find the DCA and propagate the tracks to the DCA
 201         double dca = vertexer2->PropagateToDCA(negtrack,postrack);
 202
 203         if(dca<cutDCA){
 204           // define the AliV0vertex object
 205           AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
 206           // get position of the vertex
 207           vertex2->GetXYZ(v2[0],v2[1],v2[2]);
 208           delete vertex2;
 209           //if(D0->FindV0offline(v2) && D0->d0d0()){
 210           if(D0->d0d0()){
 211             D0->SetV0(v2);
 212             D0->FindMomentaOffline();
 213             if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
 214               nD0++;
 215               event[index]=ev; index++;
 216               if(h){
 217                 h1->Fill(D0->Pt());
 218                 h3->Fill(D0->Eta());
 219               }
 220
 221               if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) {
 222                 nD0recSgn++;
 223               }
 224               else {
 225                 if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
 226                    TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
 227                   nD0recBkgS++;
 228                 } else {
 229                   nD0recBkg++;
 230                 }
 231               }
 232             }
 233           }
 234         }
 235       }
 236     }
 237
 238     //delete D0;
 239     //delete itstrks;
 240     //delete itsTree;
 241     //delete trksP;
 242     //delete itsEntryP;
 243     //delete trksN;
 244     //delete itsEntryN;
 245   }
 246
 247   cout<<"\nMy #D0: "<<nD0<<"\n"<<endl;
 248   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
 249   printf("\n+++\n+++ Total number of D0 candidates:\n
 250                  +++    Sgn:   %d\n
 251                  +++    BkgS:  %d\n
 252                  +++    Bkg:   %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
 253
 254   gBenchmark->Stop(name);
 255   gBenchmark->Show(name);
 256
 257   if(h){
 258     if(!PtD0){
 259       TCanvas *c = new TCanvas("c","c",700,1000);
 260       c->Divide(1,2);
 261       c->cd(1);
 262       h1->Draw();
 263       pt = new TPaveText(7.3,0.4,11,1.8,"br");
 264       pt->SetFillColor(0);
 265       pt->SetBorderSize(1);
 266       pt->AddText("Cuts:");
 267       Char_t st[1024];
 268       sprintf(st,"First Pt:         %g",kPtCut);
 269       pt->AddText(st);
 270       sprintf(st,"d0 low:           %g",kd0Cut);
 271       pt->AddText(st);
 272       sprintf(st,"d0 high:          %g",kd0CutHigh);
 273       pt->AddText(st);
 274       sprintf(st,"V0 low:           %g",cuts[0]);
 275       pt->AddText(st);
 276       sprintf(st,"V0 high:          %g",cuts[1]);
 277       pt->AddText(st);
 278       sprintf(st,"InvMass Diff:     %g",cuts[2]);
 279       pt->AddText(st);
 280       sprintf(st,"cosPointingAngle: %g",cuts[3]);
 281       pt->AddText(st);
 282       sprintf(st,"d0d0:             %g",cuts[4]);
 283       pt->AddText(st);
 284       sprintf(st,"cosThetaStar:     %g",cuts[5]);
 285       pt->AddText(st);
 286       sprintf(st,"PtChild:          %g",cuts[6]);
 287       pt->AddText(st);
 288       sprintf(st,"DCA:              %g",cutDCA);
 289       pt->AddText(st);
 290       pt->Draw();
 291       c_1->Modified();
 292       c->cd();
 293       c->cd(2);
 294       h3->Draw();
 295     }
 296     if(PtD0){
 297       PtD0(path,h2,h4);
 298       TCanvas *c = new TCanvas("c","c",1000,700);
 299       c->Divide(2,2);
 300       c->cd(1);
 301       h1->Draw();
 302       pt = new TPaveText(7.3,0.4,11,1.8,"br");
 303       pt->SetFillColor(0);
 304       pt->SetBorderSize(1);
 305       pt->AddText("Cuts:");
 306       Char_t st[1024];
 307       sprintf(st,"First Pt:         %g",kPtCut);
 308       pt->AddText(st);
 309       sprintf(st,"d0 low:           %g",kd0Cut);
 310       pt->AddText(st);
 311       sprintf(st,"d0 high:          %g",kd0CutHigh);
 312       pt->AddText(st);
 313       sprintf(st,"V0 low:           %g",cuts[0]);
 314       pt->AddText(st);
 315       sprintf(st,"V0 high:          %g",cuts[1]);
 316       pt->AddText(st);
 317       sprintf(st,"InvMass Diff:     %g",cuts[2]);
 318       pt->AddText(st);
 319       sprintf(st,"cosPointAng: %g",cuts[3]);
 320       pt->AddText(st);
 321       sprintf(st,"d0d0:             %g",cuts[4]);
 322       pt->AddText(st);
 323       sprintf(st,"cosTheta*:     %g",cuts[5]);
 324       pt->AddText(st);
 325       sprintf(st,"PtChild:          %g",cuts[6]);
 326       pt->AddText(st);
 327       sprintf(st,"DCA:              %g",cutDCA);
 328       pt->AddText(st);
 329       pt->Draw();
 330       c_1->Modified();
 331       c->cd();
 332       c->cd(2);
 333       h3->Draw();
 334       c->cd(3);
 335       h2->Draw();
 336       c->cd(4);
 337       h4->Draw();
 338     }
 339   }
 340   if(h){
 341     if(!PtD0){
 342       Char_t outName[1024];
 343       sprintf(outName,"%s/ReconstructedD0.root",path);
 344       TFile* outroot = new TFile(outName,"recreate");
 345       h1->Write();
 346       h3->Write();
 347       outroot->Close();
 348       delete outroot;
 349     }
 350   }
 351   if(PtD0){
 352     Char_t outName[1024];
 353     sprintf(outName,"%s/ReconstructedD0.root",path);
 354     TFile* outroot = new TFile(outName,"recreate");
 355     h1->Write();
 356     h2->Write();
 357     h3->Write();
 358     h4->Write();
 359     outroot->Close();
 360     delete outroot;
 361   }
 362   Char_t foutName[1024];
 363   sprintf(foutName,"%s/Cuts",path);
 364   ofstream fout(foutName);
 365   Char_t st2[1024];
 366   sprintf(st2,"First Pt:       %g",kPtCut);
 367   fout<<st2<<endl;
 368   sprintf(st2,"d0 low:         %g",kd0Cut);
 369   fout<<st2<<endl;
 370   sprintf(st2,"d0 high:        %g",kd0CutHigh);
 371   fout<<st2<<endl;
 372   sprintf(st2,"V0 low:         %g",cuts[0]);
 373   fout<<st2<<endl;
 374   sprintf(st2,"V0 high:        %g",cuts[1]);
 375   fout<<st2<<endl;
 376   sprintf(st2,"InvMass Diff:   %g",cuts[2]);
 377   fout<<st2<<endl;
 378   sprintf(st2,"cosPointAng:    %g",cuts[3]);
 379   fout<<st2<<endl;
 380   sprintf(st2,"d0d0:           %g",cuts[4]);
 381   fout<<st2<<endl;
 382   sprintf(st2,"cosTheta*:      %g",cuts[5]);
 383   fout<<st2<<endl;
 384   sprintf(st2,"PtChild:        %g",cuts[6]);
 385   fout<<st2<<endl;
 386   sprintf(st2,"DCA:            %g",cutDCA);
 387   fout<<st2<<endl;
 388   fout.close();
 389
 390   Char_t fName[1024];
 391   sprintf(fName,"%s/Events",path);
 392   ofstream fevent(fName);
 393   for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
 394   fevent.close();
 395 }
 396 //___________________________________________________________________________
 397 void   SelectTracks(TTree& itsTree,
 398                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
 399                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
 400 //
 401 // this function creates two TObjArrays with positive and negative tracks
 402 //
 403   nTrksP=0,nTrksN=0;
 404
 405
 406   Int_t entr = (Int_t)itsTree.GetEntries();
 407
 408   // trasfer tracks from tree to arrays
 409   for(Int_t i=0; i<entr; i++) {
 410
 411     AliITStrackV2 *itstrack = new AliITStrackV2;
 412     itsTree.SetBranchAddress("tracks",&itstrack);
 413
 414     itsTree.GetEvent(i);
 415
 416     // single track selection
 417     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
 418
 419     if(itstrack->Get1Pt()>0.) { // negative track
 420       trksN.AddLast(itstrack);
 421       itsEntryN[nTrksN] = i;
 422       nTrksN++;
 423     } else {                    // positive track
 424       trksP.AddLast(itstrack);
 425       itsEntryP[nTrksP] = i;
 426       nTrksP++;
 427     }
 428
 429   }
 430
 431   return;
 432 }
 433 //____________________________________________________________________________
 434 Bool_t TrkCuts(const AliITStrackV2& trk) {
 435 //
 436 // this function tells if track passes some kinematical cuts
 437 //
 438   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
 439   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
 440   //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
 441
 442   return kTRUE;
 443 }
 444 //____________________________________________________________________________
 445 void GetPrimaryVertex(int i,Char_t* path="./") {
 446
 447   int event=i;
 448
 449   Char_t falice[1024];
 450   sprintf(falice,"%s/galice.root",path);
 451   TFile * galice = new TFile(falice);
 452
 453   TDirectory * curdir;
 454
 455   Char_t vname[20];
 456   galice->cd();
 457
 458   sprintf(vname,"Vertex_%d",event);
 459   TArrayF o = 0;
 460   o.Set(3);
 461   AliHeader * header = 0;
 462
 463   TTree * treeE = (TTree*)gDirectory->Get("TE");
 464   treeE->SetBranchAddress("Header",&header);
 465   treeE->GetEntry(event);
 466   AliGenEventHeader* genHeader = header->GenEventHeader();
 467   if(genHeader){
 468     genHeader->PrimaryVertex(o);
 469     primaryvertex[0] = (Double_t)o[0];
 470     primaryvertex[1] = (Double_t)o[1];
 471     primaryvertex[2] = (Double_t)o[2];
 472   }
 473   else{
 474     printf("Can't find Header");
 475   }
 476   delete header;
 477   delete galice;
 478 }
 479 //____________________________________________________________________________
 480 PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
 481
 482   Char_t falice[1024];
 483   sprintf(falice,"%s/galice.root",path);
 484   TFile *f = new TFile(falice);
 485   gAlice=(AliRun*)f->Get("gAlice");
 486
 487   TParticle *p;
 488   int nd0=0;
 489   int nkminus =0;
 490   int npipluss = 0;
 491   int nEvent=gAlice->GetEventsPerRun();
 492
 493   for (Int_t i = 0; i <nEvent; i++) {
 494     cout<<"Event: "<<i<<endl;
 495     gAlice->GetEvent(i);
 496     Int_t nPart = gAlice->GetNtrack();
 497     for (Int_t iPart = 0; iPart < nPart; iPart++) {
 498       //cout<<"Particlenr.: "<<iPart<<endl;
 499       p = (TParticle*)gAlice->Particle(iPart);
 500       if (p->GetPdgCode()==421){
 501         if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
 502         h4->Fill(p.Eta());
 503         nd0++;
 504       }
 505       if (p->GetPdgCode()==-321){
 506         nkminus++;
 507       }
 508       if (p->GetPdgCode()==211){
 509         npipluss++;
 510       }
 511     }
 512   }
 513   delete gAlice;
 514 }