]>
Commit | Line | Data |
---|---|---|
af2b48fb | 1 | void ReadAOD_MCBranch(char *filein="AliAODs.root"){ |
2 | ||
3 | //============================================================================= | |
4 | // Macro to read AOD+MC Branch: | |
5 | // in this example the macro reads generated and reconstructed JPsi | |
6 | // Plots refer to all MC generated particles, to MC particles which have | |
7 | // been reconstructed and reconstructed tracks | |
8 | // | |
9 | // Each aod event has 3 MC particles (J/Psi, mu1, mu2) | |
10 | // To access MC muons: mctrackall->IsPhysicalPrimary() | |
11 | // To access JPsi muons: mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary() | |
12 | ||
13 | // AliStack::IsPhysicalPrimary() | |
14 | // Test if a particle is a physical primary according to the following definition: | |
15 | // Particles produced in the collision including products of strong and | |
16 | // electromagnetic decay and excluding feed-down from weak decays of strange | |
17 | // particles. | |
18 | //============================================================================= | |
19 | ||
20 | gStyle->SetFrameFillColor(0); | |
21 | gStyle->SetCanvasColor(0); | |
22 | gStyle->SetPalette(1); | |
23 | ||
24 | // Open AOD file | |
25 | TFile *file = new TFile(filein); | |
26 | TTree *aodTree = file->Get("aodTree"); | |
27 | ||
28 | // Read AOD event | |
29 | AliAODEvent *event = new AliAODEvent(); | |
30 | event->ReadFromTree(aodTree); | |
31 | ||
32 | // Read MC info | |
33 | TClonesArray *mcarray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); | |
34 | if(!mcarray) { | |
35 | printf("No MC info in AOD!\n"); | |
36 | return; | |
37 | } | |
38 | ||
39 | // Book Histos | |
40 | TH2D *hpt_mc_rec = new TH2D("hpt_mc_rec","Pt correlations MC-REC",50,0.,10.,50,0.,10.); | |
41 | TH1D *hpt_rec = new TH1D("hpt_rec","Pt REC tracks",20,0.,10.); | |
42 | TH1D *hpt_mc = new TH1D("hpt_mc","PT MC particles corresponding to rec tracks",20,0.,10.); | |
43 | TH1D *hpt_allmc = new TH1D("hpt_allmc","PT all MC particles",20,0.,10.); | |
44 | TH1D *hpt_primarymc = new TH1D("hpt_primarymc","PT MC primaries (J/Psi and muons)",20,0.,10.); | |
45 | TH1D *hpt_physprimarymc = new TH1D("hpt_physprimarymc","PT Muons MC",20,0.,10.); | |
46 | TH1D *hpt_notphysprimarymc = new TH1D("hpt_notphysprimarymc","Pt JPsi MC",20,0.,10.); | |
47 | TH1D *hpt_diff = new TH1D("hpt_diff","Pt: Rec - MC",100,-5.,5.); | |
48 | ||
49 | TH1D *hCharge_primarymc = new TH1D("hCharge_primarymc","Charge MC primary particle",10,-5.,5.); | |
50 | TH1D *hCharge_physprimmc = new TH1D("hCharge_physprimmc","Charge Muons MC",10,-5.,5.); | |
51 | TH1D *hCharge_notphysprimmc = new TH1D("hCharge_notphysprimmc","Charge JPsi MC",10,-5.,5.); | |
52 | TH1D *hCharge_rec = new TH1D("hCharge_rec","hCharge_rec",10,-5.,5.); | |
53 | TH1D *hCharge_mc = new TH1D("hCharge_mc","hCharge_mc",10,-5.,5.); | |
54 | ||
55 | TH1D *hMass_mc = new TH1D("hMass_mc","Mass MC particles",100,0.,5..); | |
56 | TH1D *hMass_rec = new TH1D("hMass_rec","Mass rec tracks",100,0.,5..); | |
57 | TH1D *hMass_allmc = new TH1D("hMass_allmc","Mass all MC particles",100,0.,5..); | |
58 | TH1D *hMass_primarymc = new TH1D("hMass_primarymc","Mass primary particles MC",100,0.,5..); | |
59 | TH1D *hMass_physprimarymc = new TH1D("hMass_physprimarymc","Mass muons MC",100,0.,5..); | |
60 | TH1D *hMass_notphysprimarymc = new TH1D("hMass_notphysprimarymc","Mass JPsi MC",100,0.,5..); | |
61 | TH1D *hMassJPsi_mc = new TH1D("hMassJPsi_mc","Mass JPsi MC",100,0.,5..); | |
62 | ||
63 | // Loop on AOD events | |
64 | Int_t nev=0; | |
65 | ||
66 | while(aodTree->GetEvent(nev++)){ | |
67 | ||
68 | // Loop on all MC tracks | |
69 | for(int ii=0;ii<mcarray->GetEntries();ii++){ | |
70 | AliAODMCParticle *mctrackall = (AliAODMCParticle*) mcarray->At(ii); | |
71 | // all particles | |
72 | hpt_allmc->Fill(mctrackall->Pt()); | |
73 | hMass_allmc->Fill(mctrackall->M()); | |
74 | ||
75 | // Muons | |
76 | if(mctrackall->IsPhysicalPrimary()) { | |
77 | hCharge_physprimmc->Fill(mctrackall->Charge()); | |
78 | hpt_physprimarymc->Fill(mctrackall->Pt()); | |
79 | hMass_physprimarymc->Fill(mctrackall->M()); | |
80 | } | |
81 | ||
82 | if(mctrackall->IsPrimary()) { | |
83 | hpt_primarymc->Fill(mctrackall->Pt()); | |
84 | hMass_primarymc->Fill(mctrackall->M()); | |
85 | hCharge_primarymc->Fill(mctrackall->Charge()); | |
86 | } | |
87 | ||
88 | // JPsi | |
89 | if(mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()) { | |
90 | hCharge_notphysprimmc->Fill(mctrackall->Charge()); | |
91 | hMass_notphysprimarymc->Fill(mctrackall->M()); | |
92 | hpt_notphysprimarymc->Fill(mctrackall->Pt()); | |
93 | } | |
94 | } | |
95 | ||
96 | // Loop on reconstructed tracks | |
97 | for(int i=0;i<event->GetNumberOfTracks();i++){ | |
98 | AliAODTrack *track = dynamic_cast<AliAODTrack*>(event->GetTrack(i)); | |
99 | if(!track) { | |
100 | printf("Error: no track! \n"); | |
101 | continue; | |
102 | } | |
103 | hpt_rec->Fill(track->Pt()); | |
104 | hMass_rec->Fill(track->M()); | |
105 | hCharge_rec->Fill(track->Charge()); | |
106 | ||
107 | // Read MC particle corresponding to reconstructed tracks | |
108 | AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(track->GetLabel()); | |
109 | if(!mctrack){ | |
110 | printf("Error: no MC track! \n"); | |
111 | continue; | |
112 | } | |
113 | //printf("nev=%d i=%d mc=%f rec=%f\n",nev-1,i,mctrack->Pt(),track->Pt()); | |
114 | hpt_mc->Fill(mctrack->Pt()); | |
115 | hpt_mc_rec->Fill(track->Pt(),mctrack->Pt()); | |
116 | hpt_diff->Fill(track->Pt()-mctrack->Pt()); | |
117 | hMass_mc->Fill(mctrack->M()); | |
118 | hCharge_mc->Fill(mctrack->Charge()); | |
119 | } | |
120 | } | |
121 | ||
122 | // Plots | |
123 | TCanvas *c = new TCanvas("c","PT spectra",20,20,600,600); | |
124 | c->Divide(3,3); | |
125 | c->cd(1); | |
126 | hpt_primarymc->Draw(); | |
127 | hpt_primarymc->GetXaxis()->SetTitle("Pt MC primary particles"); | |
128 | c->cd(2); | |
129 | hpt_mc->GetXaxis()->SetTitle("Pt MC particles"); | |
130 | hpt_mc->Draw(); | |
131 | hpt_mc->SetLineColor(4); | |
132 | c->cd(3); | |
133 | hpt_rec->GetXaxis()->SetTitle("Rec Pt track"); | |
134 | hpt_rec->Draw(); | |
135 | hpt_rec->SetLineColor(2); | |
136 | c->cd(4); | |
137 | hpt_physprimarymc->Draw(); | |
138 | hpt_physprimarymc->GetXaxis()->SetTitle("Pt MC muon track"); | |
139 | c->cd(5); | |
140 | hpt_notphysprimarymc->Draw(); | |
141 | hpt_notphysprimarymc->GetXaxis()->SetTitle("Pt MC JPsi track"); | |
142 | c->cd(7); | |
143 | hpt_mc_rec->Draw("colz"); | |
144 | hpt_mc_rec->GetXaxis()->SetTitle("Pt rec track"); | |
145 | hpt_mc_rec->GetYaxis()->SetTitle("Pt MC track"); | |
146 | c->cd(8); | |
147 | hpt_diff->Draw(); | |
148 | hpt_diff->GetXaxis()->SetTitle("Rec-MC Pt"); | |
149 | ||
150 | TCanvas *c1 = new TCanvas("c1","Mass distributions",40,40,600,600); | |
151 | c1->Divide(2,2); | |
152 | c1->cd(1); | |
153 | hMass_primarymc->Draw(); | |
154 | hMass_primarymc->GetXaxis()->SetTitle("Mass MC primary particles"); | |
155 | c1->cd(2); | |
156 | hMass_physprimarymc->Draw(); | |
157 | hMass_physprimarymc->GetXaxis()->SetTitle("Mass MC muon particles"); | |
158 | c1->cd(3); | |
159 | hMass_notphysprimarymc->Draw(); | |
160 | hMass_notphysprimarymc->GetXaxis()->SetTitle("Mass MC J/Psi"); | |
161 | c1->cd(4); | |
162 | hMass_rec->Draw(); | |
163 | hMass_rec->GetXaxis()->SetTitle("Mass Rec tracks"); | |
164 | ||
165 | TCanvas *c2= new TCanvas("c2","Charge distributions",60,60,600,600); | |
166 | c2->Divide(2,2); | |
167 | c2->cd(1); | |
168 | hCharge_primarymc->Draw(); | |
169 | hCharge_primarymc->GetXaxis()->SetTitle("Charge MC primary particles"); | |
170 | c2->cd(2); | |
171 | hCharge_physprimmc->Draw(); | |
172 | hCharge_physprimmc->GetXaxis()->SetTitle("Charge muon MC");; | |
173 | c2->cd(3); | |
174 | hCharge_notphysprimmc->Draw(); | |
175 | hCharge_notphysprimmc->GetXaxis()->SetTitle("Charge J/Psi MC");; | |
176 | c2->cd(4); | |
177 | hCharge_rec->Draw(); | |
178 | hCharge_rec->GetXaxis()->SetTitle("Charge rec tracks"); | |
179 | } |