PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.cxx

   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
  43 using namespace std;
  44
  45 ClassImp(AliAnalysisTaskdNdetaMC)
  46
  47 Float_t      AliAnalysisTaskdNdetaMC::fEtaMax = 0.5;
  48 Int_t        AliAnalysisTaskdNdetaMC::fPDGCodes[]  = {211,2212,321,-11,-13,-211,-2212,-321,11,13,0} ; // 0 ==> all others
  49 const char * AliAnalysisTaskdNdetaMC::fPartNames[] = {"PionPos", "ProtonPos", "KaonPos", "ePos", "muPos",
  50                                                       "PionNeg", "ProtonNeg", "KaonNeg", "eNeg", "muNeg",
  51                                                       "Others"} ;
  52
  53
  54 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC()
  55   : AliAnalysisTaskSE(), fNchDens(0), fMyOut(0),
  56     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
  57 {
  58
  59   // default constructor
  60   for(Int_t ihist = 0; ihist < kNHist; ihist++){
  61     fHistEta[ihist]  = 0;
  62     fHistPt[ihist]  = 0;
  63     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
  64       fHistMult[ihist][ihist2] = 0;
  65     }
  66     for(Int_t ipart = 0; ipart < kNPart+1; ipart++){ //  +1 => all others
  67       fHistPtID[ihist][ipart] = 0;
  68     }
  69   }
  70
  71
  72
  73
  74   fEtaBins[kEta05] = 0.5;
  75   fEtaBins[kEta10] = 1.0;
  76   fEtaBins[kEta14] = 1.3;
  77
  78
  79
  80   fEtaMax=0.5;
  81
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name)
  86   : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
  87     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
  88 {
  89   // constructor
  90
  91   for(Int_t ihist = 0; ihist < kNHist; ihist++){
  92     fHistEta[ihist]  = 0;
  93     fHistPt[ihist]  = 0;
  94     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
  95       fHistMult[ihist][ihist2] = 0;
  96     }
  97     for(Int_t ipart = 0; ipart < kNPart+1; ipart++){  // +1 ==> all others
  98       fHistPtID[ihist][ipart] = 0;
  99     }
 100   }
 101   fEtaBins[kEta05] = 0.5;
 102   fEtaBins[kEta10] = 1.0;
 103   fEtaBins[kEta14] = 1.3;
 104
 105
 106   fEtaMax=0.5;
 107
 108   AliPDG::AddParticlesToPdgDataBase();
 109
 110   DefineOutput(1, TList::Class());
 111 }
 112
 113 AliAnalysisTaskdNdetaMC::~AliAnalysisTaskdNdetaMC() {
 114
 115   // destructor
 116
 117   if(fMyOut) {
 118     // fMyOut owns the histos
 119     delete fMyOut;
 120     fMyOut = 0;
 121   }
 122
 123 }
 124
 125
 126 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name, const char * fname)
 127   : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
 128     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
 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;
 134     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 135       fHistMult[ihist][ihist2] = 0;
 136     }
 137     for(Int_t ipart = 0; ipart < kNPart +1; ipart++){ // all others
 138       fHistPtID[ihist][ipart] = 0;
 139     }
 140   }
 141   fEtaBins[kEta05] = 0.5;
 142   fEtaBins[kEta10] = 1.0;
 143   fEtaBins[kEta14] = 1.3;
 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 //________________________________________________________________________
 159 void 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
 183     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++) { // eta range
 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   }
 188   for(Int_t ihist = 0; ihist < kNHist; ihist++){ // type
 189     for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // particle (<= for all others)
 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   }
 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
 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
 223 TH1F* 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
 236 TH1F* AliAnalysisTaskdNdetaMC::BookHptHist(const char * name, const char * title) {
 237
 238   // Book Pt histos
 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);
 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
 253 TH1F* 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 //________________________________________________________________________
 268 void 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();
 276   AliGenPythiaEventHeader * headPy  = 0;
 277   AliGenDPMjetEventHeader * headPho = 0;
 278   AliGenEventHeader * htmp = mcEvent->GenEventHeader();
 279   if(!htmp) {
 280     AliError("Cannot Get MC Header!!");
 281     return;
 282   }
 283   if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
 284     headPy =  (AliGenPythiaEventHeader*) htmp;
 285   } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
 286     headPho = (AliGenDPMjetEventHeader*) htmp;
 287   } else {
 288     cout << "Unknown header" << endl;
 289
 290   }
 291
 292   if (!mcEvent) {
 293      Printf("ERROR: Could not retrieve MC event");
 294      return;
 295   }
 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
 337
 338
 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
 348   static const Float_t ymax =0.5;
 349   // Track loop
 350   Int_t multiplicity[kNHist][kNEtaHist] = {{0}};
 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);
 364
 365       Int_t particleID = kNPart;
 366       for(Int_t ipart = 0; ipart < kNPart; ipart++){
 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         }
 376       }
 377
 378       if(isEtaLess08) {
 379         fHistSpecies->Fill(particleID); //
 380         if (particleID == kNPart) {
 381           AliInfo(Form("Found other particle: [%d]", track->PdgCode()));
 382         }
 383       }
 384
 385
 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);
 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());
 392           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 393         }
 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);
 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());
 401           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 402         }
 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);
 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());
 410           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 411         }
 412
 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);
 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());
 420           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 421         }
 422
 423       }
 424
 425       // fill array of multiplicity for different classes of events
 426       // and in different eta ranges
 427       for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
 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
 443   for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
 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 //________________________________________________________________________
 456 void 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
 470 void 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++){
 488     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 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     }
 492     for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // <= , all others
 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     }
 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");
 522       for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 523         fHistMult[ihist][ihist2]->Scale(1./iev);
 524       }
 525       for(Int_t ipart = 0; ipart <= kNPart; ipart++){
 526         fHistPtID[ihist][ipart]->Scale(1./iev, "width");
 527       }
 528
 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 }