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 "AliESDEvent.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
  37   AliESDEvent * esdSig = new AliESDEvent; // The signal ESD object is put here
  38   esdSig->ReadFromTree(tSig);
  39
  40   // Run loader
  41   name = "galice.root";
  42   AliRunLoader* rlSig = AliRunLoader::Open(name.Data());
  43
  44   // gAlice
  45   rlSig->LoadgAlice();
  46   gAlice = rlSig->GetAliRun();
  47
  48   // Now load kinematics and event header
  49   rlSig->LoadKinematics();
  50   rlSig->LoadHeader();
  51
  52   // Loop on events
  53   Long64_t nevSig = rlSig->GetNumberOfEvents();
  54
  55   cout << nevSig << " events" << endl;
  56
  57   for (Int_t iev=0; iev<nevSig; iev++) {
  58     cout << "Signal event " << iev << endl;
  59
  60     // Get ESD
  61     tSig->GetEntry(iev);
  62
  63     // Get MC
  64     rlSig->GetEvent(iev);
  65
  66     // Particle stack
  67     AliStack * stackSig = rlSig->Stack();
  68
  69     Int_t nrec = esdSig->GetNumberOfTracks();
  70
  71  //-------------------------------------------------------------------------------------------------
  72     Int_t nstack = stackSig->GetNtrack();
  73     for(Int_t istack=0; istack < nstack; istack++){
  74
  75       TParticle * part = stackSig->Particle(istack);
  76
  77       // Loop on particles: check if the D* decay products are reconstructed
  78
  79       if(!part) continue;
  80
  81       if(TMath::Abs(part->GetPdgCode())== 413){
  82         cout<<"particle "<< istack << " is D*"<<endl;
  83
  84
  85         Int_t iDaughter1 = part->GetFirstDaughter();  //id of the Daughter = D^0
  86         if( iDaughter1<0) continue;
  87         TParticle* daughter1 = stackSig->Particle(iDaughter1);
  88         if(!daughter1) continue;
  89         cout<<"first daughter:  "<<daughter1->GetPdgCode()<<endl;
  90
  91         Int_t iDaughter2 = part->GetLastDaughter();  //id of the Daughter = pi+
  92         if( iDaughter2<0) continue;
  93         TParticle* daughter2 = stackSig->Particle(iDaughter2);
  94         if(!daughter2) continue;
  95         cout<<"last daughter: "<<daughter2->GetPdgCode()<<endl;
  96
  97         Int_t iD0 = -1;
  98         Int_t iPi = -1;
  99
 100         if(TMath::Abs(daughter1->GetPdgCode())== 421){
 101           iD0=iDaughter1;
 102           iPi=iDaughter2;
 103         }
 104         else if(TMath::Abs(daughter2->GetPdgCode())== 421){
 105           iD0=iDaughter2;
 106           iPi=iDaughter1;
 107         }
 108
 109         if (iD0<0)  continue;
 110
 111         TParticle* secondmother = stackSig->Particle(iD0);
 112
 113         Int_t iD0Daughter1 = secondmother->GetFirstDaughter();
 114         TParticle*  D0Daughter1 = stackSig->Particle(iD0Daughter1);
 115         Int_t iD0Daughter2 = secondmother->GetLastDaughter();
 116         TParticle*  D0Daughter2 = stackSig->Particle(iD0Daughter2);
 117
 118         for(Int_t irec=0; irec<nrec; irec++) {//loop on the ESDTree;
 119           AliESDtrack * track = esdSig->GetTrack(irec);
 120           UInt_t label = TMath::Abs(track->GetLabel());
 121           if(label<10000000) {
 122             if(label == iPi)         cout<<label<< " We found the Pi from the D* decay"<<endl;
 123             if(label == iD0Daughter1) cout<<label<<" We found the K from the D0 decay"<<endl;
 124             if(label == iD0Daughter2) cout<<label<<" We found the Pi from the D0 decay"<<endl;
 125           }
 126         }
 127       }
 128     }
 129   }
 130
 131   // end loop on kine tree
 132
 133   fSig->Close();
 134
 135   timer.Stop();
 136   timer.Print();
 137 }