AliITSv11Hybrid is replaced by AliITSv11
[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 "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 }