PWGUD/base/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,
  49                                                       -211,-2212,-321,11,13,
  50                                                       3122,-3122,3122,-3122,
  51                                                       0} ; // 0 ==> all others
  52 const char * AliAnalysisTaskdNdetaMC::fPartNames[] = {"PionPos", "ProtonPos", "KaonPos", "ePos", "muPos",
  53                                                       "PionNeg", "ProtonNeg", "KaonNeg", "eNeg", "muNeg",
  54                                                       "Lambda", "Lambda_bar", "LambdaInclusive", "Lambda_barInclusive",
  55                                                       "Others"} ;
  56
  57
  58 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC()
  59   : AliAnalysisTaskSE(), fNchDens(0), fMyOut(0),
  60     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
  61 {
  62
  63   // default constructor
  64   for(Int_t ihist = 0; ihist < kNHist; ihist++){
  65     fHistEta[ihist]  = 0;
  66     fHistPt[ihist]  = 0;
  67     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
  68       fHistMult[ihist][ihist2] = 0;
  69     }
  70     for(Int_t ipart = 0; ipart < kNPart+1; ipart++){ //  +1 => all others
  71       fHistPtID[ihist][ipart] = 0;
  72     }
  73   }
  74
  75
  76
  77
  78   fEtaBins[kEta05] = 0.5;
  79   fEtaBins[kEta10] = 1.0;
  80   fEtaBins[kEta14] = 1.3;
  81
  82
  83
  84   fEtaMax=0.5;
  85
  86 }
  87
  88 //________________________________________________________________________
  89 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name)
  90   : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
  91     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
  92 {
  93   // constructor
  94
  95   for(Int_t ihist = 0; ihist < kNHist; ihist++){
  96     fHistEta[ihist]  = 0;
  97     fHistPt[ihist]  = 0;
  98     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
  99       fHistMult[ihist][ihist2] = 0;
 100     }
 101     for(Int_t ipart = 0; ipart < kNPart+1; ipart++){  // +1 ==> all others
 102       fHistPtID[ihist][ipart] = 0;
 103     }
 104   }
 105   fEtaBins[kEta05] = 0.5;
 106   fEtaBins[kEta10] = 1.0;
 107   fEtaBins[kEta14] = 1.3;
 108
 109
 110   fEtaMax=0.5;
 111
 112   AliPDG::AddParticlesToPdgDataBase();
 113
 114   DefineOutput(1, TList::Class());
 115 }
 116
 117 AliAnalysisTaskdNdetaMC::~AliAnalysisTaskdNdetaMC() {
 118
 119   // destructor
 120
 121   if(fMyOut) {
 122     // fMyOut owns the histos
 123     delete fMyOut;
 124     fMyOut = 0;
 125   }
 126
 127 }
 128
 129
 130 AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name, const char * fname)
 131   : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
 132     fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0), fHistSpecies(0)
 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;
 138     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 139       fHistMult[ihist][ihist2] = 0;
 140     }
 141     for(Int_t ipart = 0; ipart < kNPart +1; ipart++){ // all others
 142       fHistPtID[ihist][ipart] = 0;
 143     }
 144   }
 145   fEtaBins[kEta05] = 0.5;
 146   fEtaBins[kEta10] = 1.0;
 147   fEtaBins[kEta14] = 1.3;
 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 //________________________________________________________________________
 163 void 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
 187     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++) { // eta range
 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   }
 192   for(Int_t ihist = 0; ihist < kNHist; ihist++){ // type
 193     for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // particle (<= for all others)
 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   }
 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
 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
 227 TH1F* 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
 240 TH1F* AliAnalysisTaskdNdetaMC::BookHptHist(const char * name, const char * title) {
 241
 242   // Book Pt histos
 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);
 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
 257 TH1F* 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 //________________________________________________________________________
 272 void 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();
 280   if (!mcEvent) {
 281      Printf("ERROR: Could not retrieve MC event");
 282      return;
 283   }
 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
 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     }
 322   }
 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
 341
 342
 343
 344   static const Float_t ymax =0.5;
 345   // Track loop
 346   Int_t multiplicity[kNHist][kNEtaHist] = {{0}};
 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
 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);
 394
 395       Int_t particleID = kNPart;
 396       for(Int_t ipart = 0; ipart < kNPart; ipart++){
 397         if(track->PdgCode() == fPDGCodes[ipart]) particleID = ipart; // Found one otf the expected particles, will be used to fille species histos
 398
 399         if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart]) {
 400           cout << "Found " << iTrack <<" "  << fPartNames[ipart] << endl;
 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         }
 408       }
 409
 410       if(isEtaLess08) {
 411         fHistSpecies->Fill(particleID); //
 412         if (particleID == kNPart) {
 413           AliInfo(Form("Found other particle: [%d]", track->PdgCode()));
 414         }
 415       }
 416
 417
 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);
 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());
 424           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 425         }
 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);
 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());
 433           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 434         }
 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);
 440         for(Int_t ipart = 0; ipart < kNPart; ipart++){
 441           if(track->Y() > -ymax && track->Y() < ymax && track->PdgCode() == fPDGCodes[ipart])  fHistPtID[kHistND][ipart]->Fill(track->Pt());
 442           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 443         }
 444
 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);
 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());
 452           //      else fHistPtID[kHistINEL][kNPart]->Fill(track->Pt()); // all others
 453         }
 454
 455       }
 456
 457       // fill array of multiplicity for different classes of events
 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]++;
 467
 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         }
 482       }
 483     }
 484
 485   } //track loop
 486
 487
 488
 489   // Fill multiplicity histos
 490   for(Int_t ihist = 0; ihist < kNEtaHist; ihist++){
 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 //________________________________________________________________________
 503 void 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
 517 void 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++){
 535     for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 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     }
 539     for(Int_t ipart = 0; ipart <= kNPart; ipart++){ // <= , all others
 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     }
 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");
 569       for(Int_t ihist2 = 0; ihist2 < kNEtaHist; ihist2++){
 570         fHistMult[ihist][ihist2]->Scale(1./iev);
 571       }
 572       for(Int_t ipart = 0; ipart <= kNPart; ipart++){
 573         fHistPtID[ihist][ipart]->Scale(1./iev, "width");
 574       }
 575
 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 }