3928b038 |
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 |
2eeabe87 |
18 | #include "AliESDEvent.h" |
3928b038 |
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"); |
3928b038 |
36 | |
af885e0f |
37 | AliESDEvent * esdSig = new AliESDEvent; // The signal ESD object is put here |
38 | esdSig->ReadFromTree(tSig); |
3928b038 |
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 |
af885e0f |
61 | tSig->GetEntry(iev); |
3928b038 |
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 | } |