]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/base/AliAnalysisTaskdNdetaMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGUD / base / AliAnalysisTaskdNdetaMC.cxx
CommitLineData
34d2adbe 1
2//----------------------------------------------------------------
3// Implementation of Class AliAnalysisTaskdNdetaMC
4//
5// Task used to analize simulations at generation level (i.e. only
6// needs galice.root and Kinematics.root).
7//
8// The tasks produces multiplicity, eta and pt histograms for
9// different event classes (listed in the enum in the header file).
10// The multiplicity histogram are further produced for different eta
11// ranges, defined in the constructor.
12//
13// Author: Michele Floris, CERN
14//
15//----------------------------------------------------------------
16
17
18#include "TChain.h"
19#include "TTree.h"
20#include "TH1F.h"
21#include "TCanvas.h"
22
23#include "AliAnalysisTask.h"
24#include "AliAnalysisManager.h"
25
26#include "AliVEvent.h"
27#include "AliESDEvent.h"
28#include "AliMCEvent.h"
29#include "AliAODEvent.h"
30
31#include "AliAnalysisTaskdNdetaMC.h"
32#include "TGraphErrors.h"
33#include "AliLog.h"
34#include "AliStack.h"
35#include "AliGenEventHeader.h"
36#include "AliGenPythiaEventHeader.h"
37
38#include <iostream>
39#include "AliPDG.h"
40#include "AliGenDPMjetEventHeader.h"
41#include "TFile.h"
42
43using namespace std;
44
45ClassImp(AliAnalysisTaskdNdetaMC)
46
4283b1a9 47Float_t AliAnalysisTaskdNdetaMC::fEtaMax = 0.5;
ef7f4ed2 48Int_t AliAnalysisTaskdNdetaMC::fPDGCodes[] = {211,2212,321,-11,-13,
49 -211,-2212,-321,11,13,
50 3122,-3122,3122,-3122,
51 0} ; // 0 ==> all others
32d9fb69 52const char * AliAnalysisTaskdNdetaMC::fPartNames[] = {"PionPos", "ProtonPos", "KaonPos", "ePos", "muPos",
53 "PionNeg", "ProtonNeg", "KaonNeg", "eNeg", "muNeg",
ef7f4ed2 54 "Lambda", "Lambda_bar", "LambdaInclusive", "Lambda_barInclusive",
32d9fb69 55 "Others"} ;
4283b1a9 56
34d2adbe 57
58AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC()
59 : AliAnalysisTaskSE(), fNchDens(0), fMyOut(0),
32d9fb69 60 fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
34d2adbe 61{
62
63 // default constructor
64 for(Int_t ihist = 0; ihist < kNHist; ihist++){
65 fHistEta[ihist] = 0;
66 fHistPt[ihist] = 0;
4283b1a9 67 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
34d2adbe 68 fHistMult[ihist][ihist2] = 0;
69 }
32d9fb69 70 for(Int_t ipart = 0; ipart < kNPart+1; ipart++){ // +1 => all others
4283b1a9 71 fHistPtID[ihist][ipart] = 0;
72 }
34d2adbe 73 }
74
4283b1a9 75
76
77
78 fEtaBins[kEta05] = 0.5;
79 fEtaBins[kEta10] = 1.0;
80 fEtaBins[kEta14] = 1.3;
34d2adbe 81
32d9fb69 82
83
34d2adbe 84 fEtaMax=0.5;
85
86}
87
88//________________________________________________________________________
89AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name)
90 : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
32d9fb69 91 fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
34d2adbe 92{
93 // constructor
94
95 for(Int_t ihist = 0; ihist < kNHist; ihist++){
96 fHistEta[ihist] = 0;
97 fHistPt[ihist] = 0;
4283b1a9 98 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
34d2adbe 99 fHistMult[ihist][ihist2] = 0;
100 }
32d9fb69 101 for(Int_t ipart = 0; ipart < kNPart+1; ipart++){ // +1 ==> all others
4283b1a9 102 fHistPtID[ihist][ipart] = 0;
103 }
34d2adbe 104 }
4283b1a9 105 fEtaBins[kEta05] = 0.5;
106 fEtaBins[kEta10] = 1.0;
107 fEtaBins[kEta14] = 1.3;
34d2adbe 108
109
110 fEtaMax=0.5;
111
112 AliPDG::AddParticlesToPdgDataBase();
113
114 DefineOutput(1, TList::Class());
115}
116
117AliAnalysisTaskdNdetaMC::~AliAnalysisTaskdNdetaMC() {
118
119 // destructor
120
121 if(fMyOut) {
122 // fMyOut owns the histos
123 delete fMyOut;
124 fMyOut = 0;
125 }
126
127}
128
129
130AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name, const char * fname)
131 : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
32d9fb69 132 fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
34d2adbe 133{
134 // This constructor open list from a dndeta file (useful to finalize after merging)
135 for(Int_t ihist = 0; ihist < kNHist; ihist++){
136 fHistEta[ihist] = 0;
137 fHistPt[ihist] = 0;
4283b1a9 138 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
34d2adbe 139 fHistMult[ihist][ihist2] = 0;
140 }
32d9fb69 141 for(Int_t ipart = 0; ipart < kNPart +1; ipart++){ // all others
4283b1a9 142 fHistPtID[ihist][ipart] = 0;
143 }
34d2adbe 144 }
4283b1a9 145 fEtaBins[kEta05] = 0.5;
146 fEtaBins[kEta10] = 1.0;
147 fEtaBins[kEta14] = 1.3;
34d2adbe 148
149
150 fEtaMax=0.5;
151
152 AliPDG::AddParticlesToPdgDataBase();
153
154 TFile * f =new TFile(fname);
155 if (!f) AliFatal(Form("Cannot open file %s!",fname));
156 fMyOut = (TList*) f->Get("coutput");
157 if (!fMyOut) AliFatal(Form("Cannot get output from file %s!",fname));
158
159 DefineOutput(1, TList::Class());
160}
161
162//________________________________________________________________________
163void AliAnalysisTaskdNdetaMC::UserCreateOutputObjects()
164{
165 // Create histograms
166 // Called once
167
168 fMyOut = new TList();
169
170 fHistEta[kHistINEL] = BookHetaHist("fHistdNdetaMCInel", "MC dN/deta distribution Inel");
171 fHistEta[kHistNSD] = BookHetaHist("fHistdNdetaMCNSD", "MC dN/deta distribution NSD");
172 fHistEta[kHistSiD] = BookHetaHist("fHistdNdetaMCSiD", "MC dN/deta distribution SiD");
173 fHistEta[kHistND] = BookHetaHist("fHistdNdetaMCND", "MC dN/deta distribution Non-Diffractive");
174
175 fHistEta[kHistHL] = BookHetaHist("fHistdNdetaMCHL", "MC dN/deta distribution, at least 1 particle |eta| < 1");
176
177
178 fHistPt[kHistINEL] = BookHptHist("fHistdNdptMCInel", "MC dN/dpt distribution Inel (|#eta| < 0.8)");
179 fHistPt[kHistNSD] = BookHptHist("fHistdNdptMCNSD", "MC dN/dpt distribution NSD (|#eta| < 0.8)");
180 fHistPt[kHistSiD] = BookHptHist("fHistdNdptMCSiD", "MC dN/dpt distribution SiD (|#eta| < 0.8)");
181 fHistPt[kHistND] = BookHptHist("fHistdNdptMCND", "MC dN/dpt distribution Non-Diffractive (|#eta| < 0.8)");
182 fHistPt[kHistHL] = BookHptHist("fHistdNdptMCHL", "MC dN/dpt distribution at least 1 particle |eta| < 1 (|#eta| < 0.8)");
183
184
185 const char * labelType[] = {"INEL", "NSD", "SiD", "ND", "HL"};
186 for(Int_t ihist = 0; ihist < kNHist; ihist++){ // type
4283b1a9 187 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++) { // eta range
34d2adbe 188 fHistMult[ihist][ihist2] = BookMultHisto(Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]),
189 Form("(dN/dN_{ch})_{|#eta_{max}|<%1.1f} (%s)",fEtaBins[ihist2],labelType[ihist]));
190 }
191 }
4283b1a9 192 for(Int_t ihist = 0; ihist < kNHist; ihist++){ // type
32d9fb69 193 for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // particle (<= for all others)
4283b1a9 194 fHistPtID[ihist][ipart] = BookHptHist(Form("fHistPtID_%s_%s",labelType[ihist],fPartNames[ipart]),
195 Form("fHistPtID (%s), %s - |y| < 0.5",labelType[ihist],fPartNames[ipart]));
196 }
197 }
32d9fb69 198 fHistSpecies = new TH1F ("fHistSpecies", "Species contributing to primaries (|#eta| < 0.8)", kNPart+1,-0.5, kNPart+0.5); // One bin for "others"
199 for(Int_t ibin = 1; ibin <= kNPart; ibin++){
200 fHistSpecies->GetXaxis()->SetBinLabel(ibin, fPartNames[ibin-1]);
201 }
202 fHistSpecies->GetXaxis()->SetBinLabel(kNPart+1, "Others");
203 fMyOut->Add(fHistSpecies);
204
205
34d2adbe 206
207 fNchDens = new TGraphErrors();
208 fNchDens -> SetName ("fNchDens");
209 fNchDens -> SetTitle ("Charged tracks density at mid-rapidity (|#eta| < fEtaMax)");
210 fNchDens->SetMarkerStyle(kFullCircle);
211 fMyOut->Add(fNchDens);
212
213 fHistNParticlesAtMidRapidity = new TH1I("fHistNParticlesAtMidRapidity","Number of particles at midrapidity", kNHist, -0.5, kNHist-0.5);
214 fMyOut->Add(fHistNParticlesAtMidRapidity);
215
216 fHistIev = new TH1I("fHistIev","Number of particles at midrapidity", kNHist, -0.5, kNHist-0.5);
217 fMyOut->Add(fHistIev);
218
219
220 fMyOut->SetOwner();
221
222 // Suppress annoying printout
223 AliLog::SetGlobalLogLevel(AliLog::kError);
224
225}
226
227TH1F* AliAnalysisTaskdNdetaMC::BookHetaHist(const char * name, const char * title) {
228
229 // Book Eta histos
230 TH1F * h = new TH1F(name, title, 200, -10., 10.);
231 h->GetXaxis()->SetTitle("#eta");
232 h->GetYaxis()->SetTitle("dN/d#eta");
233 h->SetMarkerStyle(kFullCircle);
234 h->Sumw2();
235 fMyOut->Add(h);
236 return h;
237
238}
239
240TH1F* AliAnalysisTaskdNdetaMC::BookHptHist(const char * name, const char * title) {
241
242 // Book Pt histos
32d9fb69 243
244 const Float_t templBins[] = {0.05,0.1,0.15,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2,2.4,2.6};
245 Int_t nbinsTempl=31;
246
247 TH1F * h = new TH1F(name, title, nbinsTempl, templBins);
34d2adbe 248 h->GetXaxis()->SetTitle("p_{T}");
249 h->GetYaxis()->SetTitle("dN/dp_{T}");
250 h->SetMarkerStyle(kFullCircle);
251 h->Sumw2();
252 fMyOut->Add(h);
253 return h;
254
255}
256
257TH1F* AliAnalysisTaskdNdetaMC::BookMultHisto(const char * name, const char * title) {
258
259 // Book multiplicity histos
260 Int_t maxmult = 100;
261 TH1F * h = new TH1F(name, title, maxmult, -0.5, maxmult-0.5);
262 h->GetXaxis()->SetTitle("N_{ch}");
263 h->GetYaxis()->SetTitle("dN/dN_{ch}");
264 h->SetMarkerStyle(kFullCircle);
265 h->Sumw2();
266 fMyOut->Add(h);
267 return h;
268
269}
270
271//________________________________________________________________________
272void AliAnalysisTaskdNdetaMC::UserExec(Option_t *)
273{
274 // Main loop
275 // Called for each event
276
277 // also a AliEvent...
278 // AliVEvent* mcEvent = MCEvent();
279 AliMCEvent* mcEvent = MCEvent();
87500cef 280 if (!mcEvent) {
281 Printf("ERROR: Could not retrieve MC event");
282 return;
283 }
34d2adbe 284 AliGenPythiaEventHeader * headPy = 0;
285 AliGenDPMjetEventHeader * headPho = 0;
286 AliGenEventHeader * htmp = mcEvent->GenEventHeader();
287 if(!htmp) {
288 AliError("Cannot Get MC Header!!");
289 return;
290 }
291 if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
292 headPy = (AliGenPythiaEventHeader*) htmp;
293 } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
294 headPho = (AliGenDPMjetEventHeader*) htmp;
295 } else {
296 cout << "Unknown header" << endl;
297
298 }
299
34d2adbe 300
301 // Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
302 // Check if the evend is single diffractive
303
304 Bool_t isSD = kFALSE;
305 Bool_t isND = kFALSE;
306 if(headPy) {
307 // cout << "Process: " << headPy->ProcessType() << endl;
308 if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
309 isSD = kTRUE; // is single difractive
310 }
311 if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
312 isND = kTRUE; // is non-diffractive
313 }
314
315 } else if (headPho) {
316 if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
317 isSD = kTRUE;
318 }
319 if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
320 isND = kTRUE;
321 }
ef7f4ed2 322 }
34d2adbe 323
324 // HL definition: is there at least one particle in |eta|<1?
325 Bool_t isThereOneCentralPart = kFALSE;
326 for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
327 AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
328 if (!track) {
329 Printf("ERROR: Could not receive track %d", iTrack);
330 continue;
331 }
332 Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
333 if (isPrimary && track->Charge() != 0){
334 if (track->Eta() > -1 && track->Eta() < 1) {
335 isThereOneCentralPart = kTRUE;
336 break;
337 }
338 }
339 }
340
ef7f4ed2 341
34d2adbe 342
343
4283b1a9 344 static const Float_t ymax =0.5;
34d2adbe 345 // Track loop
4283b1a9 346 Int_t multiplicity[kNHist][kNEtaHist] = {{0}};
ef7f4ed2 347
348 // first loop to determine multiplicity
349 for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
350 AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
351 if (!track) {
352 Printf("ERROR: Could not receive track %d", iTrack);
353 continue;
354 }
355 Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
356 if (isPrimary && track->Charge() != 0){
357 for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
358 if(track->Eta() > -fEtaBins[ihist] && track->Eta() < fEtaBins[ihist]) {
359 multiplicity[kHistINEL][ihist]++;
360 if(!isSD) multiplicity[kHistNSD][ihist]++;
361 if(isSD) multiplicity[kHistSiD][ihist]++;
362 if(isND) multiplicity[kHistND] [ihist]++;
363 if(isThereOneCentralPart) multiplicity[kHistHL] [ihist]++;
364
365 }
366 }
367 }
368 }
369
370 // if(multiplicity[kHistINEL][kEta05] < 1450 || multiplicity[kHistINEL][kEta05] > 1750) return; // FIXME: ONLY CENTRAL EVENTS 1600 +- 150
371 // fHistIev->Fill(kHistINEL);
372 // if (!isSD) fHistIev->Fill(kHistNSD);
373 // if (isSD) fHistIev->Fill(kHistSiD);
374 // if (isND) fHistIev->Fill(kHistND);
375 // if (isThereOneCentralPart) fHistIev->Fill(kHistHL);
376 // if(!(Int_t(fHistIev->GetBinContent(fHistIev->FindBin(kHistINEL)))%500))
377
378 cout << "Event " << Int_t(fHistIev->GetBinContent(fHistIev->FindBin(kHistINEL))) << endl;
379
380
34d2adbe 381 for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
382 AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
383 if (!track) {
384 Printf("ERROR: Could not receive track %d", iTrack);
385 continue;
386 }
387 Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
388 if (isPrimary && track->Charge() != 0){
389 Bool_t isEtaLess08 = (track->Eta() > -0.8 && track->Eta() < 0.8);
390
391 fHistEta[kHistINEL]->Fill(track->Eta());
392 if (isEtaLess08) fHistPt[kHistINEL] ->Fill(track->Pt());
393 if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistINEL);
4283b1a9 394
32d9fb69 395 Int_t particleID = kNPart;
4283b1a9 396 for(Int_t ipart = 0; ipart < kNPart; ipart++){
32d9fb69 397 if(track->PdgCode() == fPDGCodes[ipart]) particleID = ipart; // Found one otf the expected particles, will be used to fille species histos
ef7f4ed2 398
32d9fb69 399 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) {
ef7f4ed2 400 cout << "Found " << iTrack <<" " << fPartNames[ipart] << endl;
32d9fb69 401 fHistPtID[kHistINEL][ipart]->Fill(track->Pt());
402 }
403 else if(track->Y() > -ymax && track->Y() < ymax) {
404 // cout << "Filling others " << track->Pt()<< endl;
405 fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
406 // fHistPtID[kHistINEL][kNPart]->Print();
407 }
4283b1a9 408 }
34d2adbe 409
32d9fb69 410 if(isEtaLess08) {
411 fHistSpecies->Fill(particleID); //
412 if (particleID == kNPart) {
413 AliInfo(Form("Found other particle: [%d]", track->PdgCode()));
414 }
415 }
4283b1a9 416
417
34d2adbe 418 if(!isSD) {
419 fHistEta[kHistNSD]->Fill(track->Eta());
420 if (isEtaLess08) fHistPt [kHistNSD]->Fill(track->Pt());
421 if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistNSD);
4283b1a9 422 for(Int_t ipart = 0; ipart < kNPart; ipart++){
423 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) fHistPtID[kHistNSD][ipart]->Fill(track->Pt());
32d9fb69 424 // else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
4283b1a9 425 }
34d2adbe 426 }
427 if(isSD) {
428 fHistEta[kHistSiD]->Fill(track->Eta());
429 if (isEtaLess08) fHistPt [kHistSiD]->Fill(track->Pt());
430 if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistSiD);
4283b1a9 431 for(Int_t ipart = 0; ipart < kNPart; ipart++){
432 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) fHistPtID[kHistSiD][ipart]->Fill(track->Pt());
32d9fb69 433 // else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
4283b1a9 434 }
34d2adbe 435 }
436 if (isND) {
437 fHistEta[kHistND]->Fill(track->Eta());
438 if (isEtaLess08) fHistPt [kHistND]->Fill(track->Pt());
439 if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistND);
ef7f4ed2 440 for(Int_t ipart = 0; ipart < kNPart; ipart++){
4283b1a9 441 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) fHistPtID[kHistND][ipart]->Fill(track->Pt());
32d9fb69 442 // else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
4283b1a9 443 }
444
34d2adbe 445 }
446 if (isThereOneCentralPart) {
447 fHistEta[kHistHL]->Fill(track->Eta());
448 if (isEtaLess08) fHistPt [kHistHL]->Fill(track->Pt());
449 if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistHL);
4283b1a9 450 for(Int_t ipart = 0; ipart < kNPart; ipart++){
451 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) fHistPtID[kHistHL][ipart]->Fill(track->Pt());
32d9fb69 452 // else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
4283b1a9 453 }
454
34d2adbe 455 }
456
457 // fill array of multiplicity for different classes of events
ef7f4ed2 458 // and in different eta ranges
459 // MOVED IN A SEPARATE LOOP
460 // for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
461 // if(track->Eta() > -fEtaBins[ihist] && track->Eta() < fEtaBins[ihist]) {
462 // multiplicity[kHistINEL][ihist]++;
463 // if(!isSD) multiplicity[kHistNSD][ihist]++;
464 // if(isSD) multiplicity[kHistSiD][ihist]++;
465 // if(isND) multiplicity[kHistND] [ihist]++;
466 // if(isThereOneCentralPart) multiplicity[kHistHL] [ihist]++;
34d2adbe 467
ef7f4ed2 468 // }
469 //}
470 } // if (primary and charged)
471 // extra loop to check for lambdas
472 for(Int_t ipart = kLambda; ipart < kLambdaInclusive; ipart++){
473 // if(track->PdgCode() == fPDGCodes[ipart]) particleID = ipart; // Found one otf the expected particles, will be used to fille species histos
474 if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) {
475 // cout << "Found (LAMBDA)" << " " << iTrack << " " << isPrimary << " " << fPartNames[ipart] << " " << fPartNames[ipart+2] << endl;
476 if(isPrimary) {
477 // cout << " - prim" << endl;
478 fHistPtID[kHistINEL][ipart]->Fill(track->Pt());
479 } else if(mcEvent->Stack()->IsSecondaryFromWeakDecay(iTrack)) {
480 fHistPtID[kHistINEL][ipart+2]->Fill(track->Pt()); // inclusive lambdas are offset by 2 elements
481 }
34d2adbe 482 }
483 }
ef7f4ed2 484
34d2adbe 485 } //track loop
486
487
488
489 // Fill multiplicity histos
4283b1a9 490 for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
34d2adbe 491 fHistMult[kHistINEL][ihist] ->Fill(multiplicity[kHistINEL][ihist]);
492 if(!isSD) fHistMult[kHistNSD][ihist]->Fill(multiplicity[kHistNSD][ihist]);
493 if(isSD) fHistMult[kHistSiD][ihist]->Fill(multiplicity[kHistSiD][ihist]);
494 if(isND) fHistMult[kHistND] [ihist]->Fill(multiplicity[kHistND] [ihist]);
495 if(isThereOneCentralPart) fHistMult[kHistHL] [ihist]->Fill(multiplicity[kHistHL] [ihist]);
496 }
497
498 // Post output data.
499 PostData(1, fMyOut);
500}
501
502//________________________________________________________________________
503void AliAnalysisTaskdNdetaMC::Terminate(Option_t *)
504{
505 // Draw result to the screen
506 // Called once at the end of the query
507
508 // fHistPt = dynamic_cast<TH1F*> (GetOutputData(1));
509 fMyOut = dynamic_cast<TList*> (GetOutputData(1));
510
511 Finalize();
512
513}
514
515
516
517void AliAnalysisTaskdNdetaMC::Finalize() {
518
519 // Scale histos and computes dNdeta
520
521 fHistEta[kHistINEL] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCInel");
522 fHistEta[kHistNSD] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCNSD" );
523 fHistEta[kHistSiD] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCSiD" );
524 fHistEta[kHistND] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCND" );
525 fHistEta[kHistHL] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCHL" );
526 fHistPt[kHistINEL] = (TH1F*) fMyOut->FindObject("fHistdNdptMCInel");
527 fHistPt[kHistNSD] = (TH1F*) fMyOut->FindObject("fHistdNdptMCNSD" );
528 fHistPt[kHistSiD] = (TH1F*) fMyOut->FindObject("fHistdNdptMCSiD" );
529 fHistPt[kHistND] = (TH1F*) fMyOut->FindObject("fHistdNdptMCND" );
530 fHistPt[kHistHL] = (TH1F*) fMyOut->FindObject("fHistdNdptMCHL" );
531 fNchDens = (TGraphErrors*) fMyOut->FindObject("fNchDens");
532
533 const char * labelType[] = {"INEL", "NSD", "SiD", "ND", "HL"};
534 for(Int_t ihist = 0; ihist < kNHist; ihist++){
4283b1a9 535 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
34d2adbe 536 fHistMult[ihist][ihist2] = (TH1F*) fMyOut->FindObject(Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]));
537 if (!fHistMult[ihist][ihist2]) cout << "Cannot get histo " << Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]) << endl;
538 }
32d9fb69 539 for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // <= , all others
4283b1a9 540 fHistPtID[ihist][ipart] = (TH1F*) fMyOut->FindObject(Form("fHistPtID_%s_%s",labelType[ihist],fPartNames[ipart]));
541 if(!fHistPtID[ihist][ipart]) AliWarning(Form("Cannot get histo fHistPtID_%s_%s",labelType[ihist],fPartNames[ipart]));
542 }
34d2adbe 543 }
544
545 fHistIev = (TH1I*) fMyOut->FindObject("fHistIev" );
546 fHistNParticlesAtMidRapidity = (TH1I*) fMyOut->FindObject("fHistNParticlesAtMidRapidity" );
547
548
549// fHistIev->Draw();
550
551 cout << "Eta Max: " << fEtaMax << endl;
552
553
554 for(Int_t ihist = 0; ihist < kNHist; ihist++){
555
556 Int_t iev = (Int_t) fHistIev->GetBinContent(fHistIev->FindBin(ihist));
557 Int_t npm = (Int_t) fHistNParticlesAtMidRapidity->GetBinContent(fHistNParticlesAtMidRapidity->FindBin(ihist));
558
559 // Renormalize dNdeta to the number of events
560 // cout << fHistEta[ihist] << " " << iev << endl;
561
562 // compute density at midrapidity (Delta_y = 1);
563 Int_t firstbin = fHistEta[ihist]->FindBin(-fEtaMax);
564 Int_t lastbin = fHistEta[ihist]->FindBin(fEtaMax);
565
566 if (!fSkipNormalization) {
567 fHistEta[ihist]->Scale(1./iev, "width");
568 fHistPt[ihist] ->Scale(1./iev, "width");
4283b1a9 569 for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
34d2adbe 570 fHistMult[ihist][ihist2]->Scale(1./iev);
571 }
32d9fb69 572 for(Int_t ipart = 0; ipart <= kNPart; ipart++){
4283b1a9 573 fHistPtID[ihist][ipart]->Scale(1./iev, "width");
574 }
575
34d2adbe 576 }
577
578
579 Double_t meaneta = 0;
580 Double_t meanerr = 0;
581 Double_t sumweight = 0;
582 for(Int_t ibin = firstbin; ibin <=lastbin ; ibin++){
583 Double_t x = fHistEta[ihist]->GetBinContent(ibin);
584 Double_t xerr = fHistEta[ihist]->GetBinError(ibin);
585
586
587 Double_t xerr2 = xerr*xerr;
588 if(xerr2){
589 // cout << "xe2 " << xerr2 << endl;
590 Double_t weight = 1. / xerr2;
591 sumweight += weight;
592 meaneta += weight * x;
593 }
594
595 }
596
597 if(sumweight){
598 meaneta /= sumweight;
599 meanerr = TMath::Sqrt(1./ sumweight);
600 }
601 else {
602 meaneta = 0;
603 meanerr = 0;
604 }
605
606
607 Double_t range = 2*fEtaMax;
608
609// meaneta /= fHistEta[ihist]->GetBinWidth(firstbin);
610// meanerr /= fHistEta[ihist]->GetBinWidth(firstbin);
611 cout << "Histo: " << fHistEta[ihist]->GetName() << endl;
612 cout << " Evts: " << iev << endl;
613
614 cout << " Density at midrapidity: " << meaneta << "+-" << meanerr << endl;
615
616 // Direct computation
617 Double_t errdir = TMath::Sqrt(npm) / iev;
618 Double_t etadir = Double_t(npm) / iev;
619
620 cout << " Density at midrapidity: (dir) " << etadir/range << "+-" << errdir/range << endl;
621
622 fNchDens->SetPoint (ihist, ihist+1, etadir);
623 fNchDens->SetPointError(ihist, 0, errdir);
624
625 cout << " Density at midrapidity: (TH1, Eta<0.5) " << fHistMult[ihist][0]->GetMean() << "+-" << fHistMult[ihist][0]->GetMeanError() << endl;
626
627 }
628
629 // Draw fHistEtaAll
630 TCanvas *c1 = new TCanvas("AliAnalysisTaskdNdetaMC","dNdetaMC",10,10,800,400);
631 // c1->cd(1)->SetLogy();
632 c1->Divide(2,1);
633 c1->cd(1);
634 fHistEta[0]->DrawCopy("");
635 fHistEta[1]->SetLineColor(kRed);
636 fHistEta[1]->SetMarkerColor(kRed);
637 fHistEta[1]->SetMarkerStyle(kOpenSquare);
638 fHistEta[1]->DrawCopy("same");
639 c1->cd(2);
640 fNchDens->Draw("AP");
641
642}