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