test/genkine/sim/test.C

   1 // Usage in compiled mode
   2 // gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include");
   3 // gROOT->LoadMacro("test.C+");
   4 // test()
   5
   6 #if !defined(__CINT__) || defined(__MAKECINT__)
   7
   8 // Root include files
   9 #include <Riostream.h>
  10 #include <TFile.h>
  11 #include <TTree.h>
  12 #include <TBranch.h>
  13 #include <TStopwatch.h>
  14 #include <TObject.h>
  15 #include <TParticle.h>
  16
  17 // AliRoot include files
  18 #include "AliESD.h"
  19 #include "AliRunLoader.h"
  20 #include "AliRun.h"
  21 #include "AliStack.h"
  22
  23 #endif
  24
  25 void test() {
  26
  27   TStopwatch timer;
  28   timer.Start();
  29
  30   TString name;
  31
  32   // Signal file, tree, and branch
  33   name = "AliESDs.root";
  34   TFile * fSig = TFile::Open(name.Data());
  35   TTree * tSig = (TTree*)fSig->Get("esdTree");
  36   TBranch * bSig = tSig->GetBranch("ESD");
  37
  38   AliESD * esdSig = 0; // The signal ESD object is put here
  39   bSig->SetAddress(&esdSig);
  40
  41   // Run loader
  42   name = "galice.root";
  43   AliRunLoader* rlSig = AliRunLoader::Open(name.Data());
  44
  45   // gAlice
  46   rlSig->LoadgAlice();
  47   gAlice = rlSig->GetAliRun();
  48
  49   // Now load kinematics and event header
  50   rlSig->LoadKinematics();
  51   rlSig->LoadHeader();
  52
  53   // Loop on events
  54   Long64_t nevSig = rlSig->GetNumberOfEvents();
  55
  56   cout << nevSig << " events" << endl;
  57
  58   for (Int_t iev=0; iev<nevSig; iev++) {
  59     cout << "Signal event " << iev << endl;
  60
  61     // Get ESD
  62     bSig->GetEntry(iev);
  63
  64     // Get MC
  65     rlSig->GetEvent(iev);
  66
  67     // Particle stack
  68     AliStack * stackSig = rlSig->Stack();
  69
  70     Int_t nrec = esdSig->GetNumberOfTracks();
  71
  72  //-------------------------------------------------------------------------------------------------
  73     Int_t nstack = stackSig->GetNtrack();
  74     for(Int_t istack=0; istack < nstack; istack++){
  75
  76       TParticle * part = stackSig->Particle(istack);
  77
  78       // Loop on particles: check if the D* decay products are reconstructed
  79
  80       if(!part) continue;
  81
  82       if(TMath::Abs(part->GetPdgCode())== 413){
  83         cout<<"particle "<< istack << " is D*"<<endl;
  84
  85
  86         Int_t iDaughter1 = part->GetFirstDaughter();  //id of the Daughter = D^0
  87         if( iDaughter1<0) continue;
  88         TParticle* daughter1 = stackSig->Particle(iDaughter1);
  89         if(!daughter1) continue;
  90         cout<<"first daughter:  "<<daughter1->GetPdgCode()<<endl;
  91
  92         Int_t iDaughter2 = part->GetLastDaughter();  //id of the Daughter = pi+
  93         if( iDaughter2<0) continue;
  94         TParticle* daughter2 = stackSig->Particle(iDaughter2);
  95         if(!daughter2) continue;
  96         cout<<"last daughter: "<<daughter2->GetPdgCode()<<endl;
  97
  98         Int_t iD0 = -1;
  99         Int_t iPi = -1;
 100
 101         if(TMath::Abs(daughter1->GetPdgCode())== 421){
 102           iD0=iDaughter1;
 103           iPi=iDaughter2;
 104         }
 105         else if(TMath::Abs(daughter2->GetPdgCode())== 421){
 106           iD0=iDaughter2;
 107           iPi=iDaughter1;
 108         }
 109
 110         if (iD0<0)  continue;
 111
 112         TParticle* secondmother = stackSig->Particle(iD0);
 113
 114         Int_t iD0Daughter1 = secondmother->GetFirstDaughter();
 115         TParticle*  D0Daughter1 = stackSig->Particle(iD0Daughter1);
 116         Int_t iD0Daughter2 = secondmother->GetLastDaughter();
 117         TParticle*  D0Daughter2 = stackSig->Particle(iD0Daughter2);
 118
 119         for(Int_t irec=0; irec<nrec; irec++) {//loop on the ESDTree;
 120           AliESDtrack * track = esdSig->GetTrack(irec);
 121           UInt_t label = TMath::Abs(track->GetLabel());
 122           if(label<10000000) {
 123             if(label == iPi)         cout<<label<< " We found the Pi from the D* decay"<<endl;
 124             if(label == iD0Daughter1) cout<<label<<" We found the K from the D0 decay"<<endl;
 125             if(label == iD0Daughter2) cout<<label<<" We found the Pi from the D0 decay"<<endl;
 126           }
 127         }
 128       }
 129     }
 130   }
 131
 132   // end loop on kine tree
 133
 134   fSig->Close();
 135
 136   timer.Stop();
 137   timer.Print();
 138 }