Structure modules as built
[u/mrichter/AliRoot.git] / 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 }