PWGPP/EVCHAR/AliAnalysisTaskCentrality.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //                                                         //
  18 //      Class to analyze centrality measurements           //
  19 //                                                         //
  20 /////////////////////////////////////////////////////////////
  21
  22 #include <TTree.h>
  23 #include <TList.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TProfile.h>
  27 #include <TFile.h>
  28 #include <TString.h>
  29 #include <TCanvas.h>
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliVEvent.h"
  33 #include "AliESD.h"
  34 #include "AliESDEvent.h"
  35 #include "AliESDHeader.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliESDZDC.h"
  38 #include "AliESDFMD.h"
  39 #include "AliESDVZERO.h"
  40 #include "AliMultiplicity.h"
  41 #include "AliAODHandler.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODVertex.h"
  44 #include "AliAODMCHeader.h"
  45 #include "AliMCEvent.h"
  46 #include "AliMCEventHandler.h"
  47 #include "AliMCParticle.h"
  48 #include "AliStack.h"
  49 #include "AliHeader.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliAnalysisTaskSE.h"
  52 #include "AliGenEventHeader.h"
  53 #include "AliGenHijingEventHeader.h"
  54 #include "AliPhysicsSelectionTask.h"
  55 #include "AliPhysicsSelection.h"
  56 #include "AliBackgroundSelection.h"
  57 #include "AliAnalysisTaskCentrality.h"
  58
  59 ClassImp(AliAnalysisTaskCentrality)
  60
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskCentrality::AliAnalysisTaskCentrality():
  64   AliAnalysisTaskSE(),
  65   fDebug(0),
  66   fAnalysisInput("ESD"),
  67   fIsMCInput(kFALSE),
  68   fOutput(0x0),
  69   fhEzdc(0x0),
  70   fhEzem(0x0),
  71   fhNtracks(0x0),
  72   fhNtracklets(0x0),
  73   fhNclusters0(0x0),
  74   fhmultV0(0x0),
  75   fhmultFMD(0x0),
  76   fhEzemvsEzdc(0x0),
  77   fhNtracksvsEzdc(0x0),
  78   fhNtrackletsvsEzdc(0x0),
  79   fhNclusters0vsEzdc(0x0),
  80   fhmultV0vsEzdc(0x0),
  81   fhmultFMDvsEzdc(0x0),
  82   fhNtracksvsEzem(0x0),
  83   fhNtrackletsvsEzem(0x0),
  84   fhNclusters0vsEzem(0x0),
  85   fhmultV0vsEzem(0x0),
  86   fhmultFMDvsEzem(0x0),
  87   fhNtracksvsmultV0(0x0),
  88   fhNtrackletsvsmultV0(0x0),
  89   fhNclusters0vsmultV0(0x0),
  90   fhNtracksvsmultFMD(0x0),
  91   fhNtrackletsvsmultFMD(0x0),
  92   fhNclusters0vsmultFMD(0x0),
  93   fhmultV0vsmultFMD(0x0),
  94   fNev(0),
  95   fBeamEnergy(0),
  96   fNmyTracksgen(0),
  97   fTrigClass(),
  98   fxVertex(0),
  99   fyVertex(0),
 100   fzVertex(0),
 101   fVertexer3d(0),
 102   fbMC(0),
 103   fNpartTargMC(0),
 104   fNpartProjMC(0),
 105   fNNColl(0),
 106   fNNwColl(0),
 107   fNwNColl(0),
 108   fNwNwColl(0),
 109   fNTracklets(0),
 110   fNSingleClusters(0),
 111   fNClusters(),
 112   fNChips(),
 113   fbZDC(0),
 114   fNpartZDC(0),
 115   fbZDCA(0),
 116   fNpartZDCA(0),
 117   fbZDCC(0),
 118   fNpartZDCC(0),
 119   fESDFlag(0),
 120   fZNCEnergy(0),
 121   fZPCEnergy(0),
 122   fZNAEnergy(0),
 123   fZPAEnergy(0),
 124   fZEM1Energy(0),
 125   fZEM2Energy(0),
 126   fZNCtower(),
 127   fZPCtower(),
 128   fZNAtower(),
 129   fZPAtower(),
 130   fCentrZNC(),
 131   fCentrZNA(),
 132   fNTracks(0),
 133   fNPmdTracks(0),
 134   fMultV0A(0),
 135   fMultV0C(0),
 136   fMultFMDA(0),
 137   fMultFMDC(0)
 138 {
 139    // Default constructor
 140 }
 141
 142 //________________________________________________________________________
 143 AliAnalysisTaskCentrality::AliAnalysisTaskCentrality(const char *name):
 144   AliAnalysisTaskSE(name),
 145   fDebug(0),
 146   fAnalysisInput("ESD"),
 147   fIsMCInput(kFALSE),
 148   fOutput(0x0),
 149   fhEzdc(0x0),
 150   fhEzem(0x0),
 151   fhNtracks(0x0),
 152   fhNtracklets(0x0),
 153   fhNclusters0(0x0),
 154   fhmultV0(0x0),
 155   fhmultFMD(0x0),
 156   fhEzemvsEzdc(0x0),
 157   fhNtracksvsEzdc(0x0),
 158   fhNtrackletsvsEzdc(0x0),
 159   fhNclusters0vsEzdc(0x0),
 160   fhmultV0vsEzdc(0x0),
 161   fhmultFMDvsEzdc(0x0),
 162   fhNtracksvsEzem(0x0),
 163   fhNtrackletsvsEzem(0x0),
 164   fhNclusters0vsEzem(0x0),
 165   fhmultV0vsEzem(0x0),
 166   fhmultFMDvsEzem(0x0),
 167   fhNtracksvsmultV0(0x0),
 168   fhNtrackletsvsmultV0(0x0),
 169   fhNclusters0vsmultV0(0x0),
 170   fhNtracksvsmultFMD(0x0),
 171   fhNtrackletsvsmultFMD(0x0),
 172   fhNclusters0vsmultFMD(0x0),
 173   fhmultV0vsmultFMD(0x0),
 174   fNev(0),
 175   fBeamEnergy(0),
 176   fNmyTracksgen(0),
 177   fTrigClass(),
 178   fxVertex(0),
 179   fyVertex(0),
 180   fzVertex(0),
 181   fVertexer3d(0),
 182   fbMC(0),
 183   fNpartTargMC(0),
 184   fNpartProjMC(0),
 185   fNNColl(0),
 186   fNNwColl(0),
 187   fNwNColl(0),
 188   fNwNwColl(0),
 189   fNTracklets(0),
 190   fNSingleClusters(0),
 191   fNClusters(),
 192   fNChips(),
 193   fbZDC(0),
 194   fNpartZDC(0),
 195   fbZDCA(0),
 196   fNpartZDCA(0),
 197   fbZDCC(0),
 198   fNpartZDCC(0),
 199   fESDFlag(0),
 200   fZNCEnergy(0),
 201   fZPCEnergy(0),
 202   fZNAEnergy(0),
 203   fZPAEnergy(0),
 204   fZEM1Energy(0),
 205   fZEM2Energy(0),
 206   fZNCtower(),
 207   fZPCtower(),
 208   fZNAtower(),
 209   fZPAtower(),
 210   fCentrZNC(),
 211   fCentrZNA(),
 212   fNTracks(0),
 213   fNPmdTracks(0),
 214   fMultV0A(0),
 215   fMultV0C(0),
 216   fMultFMDA(0),
 217   fMultFMDC(0)
 218 {
 219   // Default constructor
 220
 221   // Output slot #1 writes into a TList container
 222   DefineOutput(1, TList::Class());
 223
 224 }
 225
 226 //________________________________________________________________________
 227 AliAnalysisTaskCentrality& AliAnalysisTaskCentrality::operator=(const AliAnalysisTaskCentrality& c)
 228 {
 229   //
 230   // Assignment operator
 231   //
 232   if (this!=&c) {
 233     AliAnalysisTaskSE::operator=(c);
 234   }
 235   return *this;
 236 }
 237
 238 //________________________________________________________________________
 239 AliAnalysisTaskCentrality::AliAnalysisTaskCentrality(const AliAnalysisTaskCentrality& ana):
 240   AliAnalysisTaskSE(ana),
 241   fDebug(ana.fDebug),
 242   fAnalysisInput(ana.fDebug),
 243   fIsMCInput(ana.fIsMCInput),
 244   fOutput(ana.fOutput),
 245   fhEzdc(ana.fhEzdc),
 246   fhEzem(ana.fhEzem),
 247   fhNtracks(ana.fhNtracks),
 248   fhNtracklets(ana.fhNtracklets),
 249   fhNclusters0(ana.fhNclusters0),
 250   fhmultV0(ana.fhmultV0),
 251   fhmultFMD(ana.fhmultFMD),
 252   fhEzemvsEzdc(ana.fhEzemvsEzdc),
 253   fhNtracksvsEzdc(ana.fhNtracksvsEzdc),
 254   fhNtrackletsvsEzdc(ana.fhNtrackletsvsEzdc),
 255   fhNclusters0vsEzdc(ana.fhNclusters0vsEzdc),
 256   fhmultV0vsEzdc(ana.fhmultV0vsEzdc),
 257   fhmultFMDvsEzdc(ana.fhmultFMDvsEzdc),
 258   fhNtracksvsEzem(ana.fhNtracksvsEzem),
 259   fhNtrackletsvsEzem(ana.fhNtrackletsvsEzem),
 260   fhNclusters0vsEzem(ana.fhNclusters0vsEzem),
 261   fhmultV0vsEzem(ana.fhmultV0vsEzem),
 262   fhmultFMDvsEzem(ana.fhmultFMDvsEzem),
 263   fhNtracksvsmultV0(ana.fhNtracksvsmultV0),
 264   fhNtrackletsvsmultV0(ana.fhNtrackletsvsmultV0),
 265   fhNclusters0vsmultV0(ana.fhNclusters0vsmultV0),
 266   fhNtracksvsmultFMD(ana.fhNtracksvsmultFMD),
 267   fhNtrackletsvsmultFMD(ana.fhNtrackletsvsmultFMD),
 268   fhNclusters0vsmultFMD(ana.fhNclusters0vsmultFMD),
 269   fhmultV0vsmultFMD(ana.fhmultV0vsmultFMD),
 270   fNev(ana.fNev),
 271   fBeamEnergy(ana.fBeamEnergy),
 272   fNmyTracksgen(ana.fNmyTracksgen),
 273   fxVertex(ana.fxVertex),
 274   fyVertex(ana.fyVertex),
 275   fzVertex(ana.fzVertex),
 276   fVertexer3d(ana.fVertexer3d),
 277   fbMC(ana.fbMC),
 278   fNpartTargMC(ana.fNpartTargMC),
 279   fNpartProjMC(ana.fNpartProjMC),
 280   fNNColl(ana.fNNColl),
 281   fNNwColl(ana.fNNwColl),
 282   fNwNColl(ana.fNwNColl),
 283   fNwNwColl(ana.fNwNwColl),
 284   fNTracklets(ana.fNTracklets),
 285   fNSingleClusters(ana.fNSingleClusters),
 286   fbZDC(ana.fbZDC),
 287   fNpartZDC(ana.fNpartZDC),
 288   fbZDCA(ana.fbZDCA),
 289   fNpartZDCA(ana.fNpartZDCA),
 290   fbZDCC(ana.fbZDCC),
 291   fNpartZDCC(ana.fNpartZDCC),
 292   fESDFlag(ana.fESDFlag),
 293   fZNCEnergy(ana.fZNCEnergy),
 294   fZPCEnergy(ana.fZPCEnergy),
 295   fZNAEnergy(ana.fZNAEnergy),
 296   fZPAEnergy(ana.fZPAEnergy),
 297   fZEM1Energy(ana.fZEM1Energy),
 298   fZEM2Energy(ana.fZEM2Energy),
 299   fZNCtower(),
 300   fZPCtower(),
 301   fZNAtower(),
 302   fZPAtower(),
 303   fCentrZNC(),
 304   fCentrZNA(),
 305   fNTracks(ana.fNTracks),
 306   fNPmdTracks(ana.fNPmdTracks),
 307   fMultV0A(ana.fMultV0A),
 308   fMultV0C(ana.fMultV0C),
 309   fMultFMDA(ana.fMultFMDA),
 310   fMultFMDC(ana.fMultFMDC)
 311 {
 312   //
 313   // Copy Constructor
 314   //
 315
 316   for (Int_t i=0; i<5; ++i) {
 317     fZNCtower[i] = ana.fZNCtower[i];
 318     fZPCtower[i] = ana.fZPCtower[i];
 319     fZNAtower[i] = ana.fZNAtower[i];
 320     fZPAtower[i] = ana.fZPAtower[i];
 321   }
 322   for (Int_t i=0; i<2; ++i) {
 323     fCentrZNC[i] = ana.fCentrZNC[i];
 324     fCentrZNA[i] = ana.fCentrZNA[i];
 325   }
 326 }
 327
 328 //________________________________________________________________________
 329  AliAnalysisTaskCentrality::~AliAnalysisTaskCentrality()
 330  {
 331    // Destructor
 332    if(fOutput){
 333      delete fOutput; fOutput=0;
 334    }
 335  }
 336
 337 //________________________________________________________________________
 338 void AliAnalysisTaskCentrality::UserCreateOutputObjects()
 339 {
 340
 341   // Create the output containers
 342   if(fDebug>1) printf("AnalysisTaskZDCpp::UserCreateOutputObjects() \n");
 343
 344   // Several histograms are more conveniently managed in a TList
 345   fOutput = new TList();
 346   fOutput->SetOwner();
 347   fOutput->SetName("OutputHistos");
 348
 349   fhEzdc         = new TH1F("hEzdc","hEzdc",500,0,150);
 350   fhEzem         = new TH1F("hEzem","hEzem",500,0,5);
 351   fhNtracks      = new TH1F("hNtracks","hNtracks",500,0,17000);
 352   fhNtracklets   = new TH1F("hNtracklets","hNtracklets",500,0,10000);
 353   fhNclusters0   = new TH1F("hNclusters0","hNclusters0",500,0,15000);
 354   fhmultV0       = new TH1F("hmultV0","hmultV0",500,0,30000);
 355   fhmultFMD      = new TH1F("hmultFMD","hmultFMD",500,0,24000);
 356
 357   fhEzemvsEzdc         = new TProfile("hEzemvsEzdc","hEzemvsEzdc",500,0,5,"");
 358   fhNtracksvsEzdc      = new TProfile("hNtracksvsEzdc","hNtracksvsEzdc",500,0,17000,"");
 359   fhNtrackletsvsEzdc   = new TProfile("hNtrackletsvsEzdc","hNtrackletsvsEzdc",500,0,10000,"");
 360   fhNclusters0vsEzdc   = new TProfile("hNclusters0vsEzdc","hNclusters0vsEzdc",500,0,15000,"");
 361   fhmultV0vsEzdc       = new TProfile("hmultV0vsEzdc","hmultV0vsEzdc",500,0,30000,"");
 362   fhmultFMDvsEzdc      = new TProfile("hmultFMDvsEzdc","hmultFMDvsEzdc",500,0,24000,"");
 363   fhNtracksvsEzem      = new TProfile("hNtracksvsEzem","hNtracksvsEzem",500,0,17000,"");
 364   fhNtrackletsvsEzem   = new TProfile("hNtrackletsvsEzem","hNtrackletsvsEzem",500,0,10000,"");
 365   fhNclusters0vsEzem   = new TProfile("hNclusters0vsEzem","hNclusters0vsEzem",500,0,15000,"");
 366   fhmultV0vsEzem       = new TProfile("hmultV0vsEzem","hmultV0vsEzem",500,0,30000,"");
 367   fhmultFMDvsEzem      = new TProfile("hmultFMDvsEzem","hmultFMDvsEzem",500,0,24000,"");
 368   fhNtracksvsmultV0    = new TProfile("hNtracksvsmultV0","hNtracksvsmultV0",500,0,17000,"");
 369   fhNtrackletsvsmultV0 = new TProfile("hNtrackletsvsmultV0","hNtrackletsvsmultV0",500,0,10000,"");
 370   fhNclusters0vsmultV0 = new TProfile("hNclusters0vsmultV0","hNclusters0vsmultV0",500,0,15000,"");
 371   fhNtracksvsmultFMD   = new TProfile("hNtracksvsmultFMD","hNtracksvsmultFMD",500,0,17000,"");
 372   fhNtrackletsvsmultFMD= new TProfile("hNtrackletsvsmultFMD","hNtrackletsvsmultFMD",500,0,10000,"");
 373   fhNclusters0vsmultFMD= new TProfile("hNclusters0vsmultFMD","hNclusters0vsmultFMD",500,0,15000,"");
 374   fhmultV0vsmultFMD    = new TProfile("hmultV0vsmultFMD","hmultV0vsmultFMD",500,0,30000,"");
 375
 376   fhEzdc         ->GetXaxis()->SetTitle("E_{ZDC}[TeV]");
 377   fhEzem         ->GetXaxis()->SetTitle("E_{ZEM}[TeV]");
 378   fhNtracks      ->GetXaxis()->SetTitle("N_{tracks}");
 379   fhNtracklets   ->GetXaxis()->SetTitle("N_{tracklets}");
 380   fhNclusters0   ->GetXaxis()->SetTitle("N_{clusters0}");
 381   fhmultV0       ->GetXaxis()->SetTitle("V0 mult");
 382   fhmultFMD      ->GetXaxis()->SetTitle("FMD mult");
 383
 384   fhEzemvsEzdc         ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 385   fhNtracksvsEzdc      ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 386   fhNtrackletsvsEzdc   ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 387   fhNclusters0vsEzdc   ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 388   fhmultV0vsEzdc       ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 389   fhmultFMDvsEzdc      ->GetYaxis()->SetTitle("E_{ZDC}[TeV]");
 390   fhNtracksvsEzem      ->GetYaxis()->SetTitle("E_{ZEM}[TeV]");
 391   fhNtrackletsvsEzem   ->GetYaxis()->SetTitle("E_{ZEM}[TeV]");
 392   fhNclusters0vsEzem   ->GetYaxis()->SetTitle("E_{ZEM}[TeV]");
 393   fhmultV0vsEzem       ->GetYaxis()->SetTitle("E_{ZEM}[TeV]");
 394   fhmultFMDvsEzem      ->GetYaxis()->SetTitle("E_{ZEM}[TeV]");
 395   fhNtracksvsmultV0    ->GetYaxis()->SetTitle("V0 mult");
 396   fhNtrackletsvsmultV0 ->GetYaxis()->SetTitle("V0 mult");
 397   fhNclusters0vsmultV0 ->GetYaxis()->SetTitle("V0 mult");
 398   fhNtracksvsmultFMD   ->GetYaxis()->SetTitle("FMD mult");
 399   fhNtrackletsvsmultFMD->GetYaxis()->SetTitle("FMD mult");
 400   fhNclusters0vsmultFMD->GetYaxis()->SetTitle("FMD mult");
 401   fhmultV0vsmultFMD    ->GetYaxis()->SetTitle("FMD mult");
 402
 403   fhEzemvsEzdc         ->GetXaxis()->SetTitle("E_{ZEM}[TeV]");
 404   fhNtracksvsEzdc      ->GetXaxis()->SetTitle("N_{tracks}");
 405   fhNtrackletsvsEzdc   ->GetXaxis()->SetTitle("N_{tracklets}");
 406   fhNclusters0vsEzdc   ->GetXaxis()->SetTitle("N_{clusters0}");
 407   fhmultV0vsEzdc       ->GetXaxis()->SetTitle("V0 mult");
 408   fhmultFMDvsEzdc      ->GetXaxis()->SetTitle("FMD mult");
 409   fhNtracksvsEzem      ->GetXaxis()->SetTitle("N_{tracks}");
 410   fhNtrackletsvsEzem   ->GetXaxis()->SetTitle("N_{tracklets}");
 411   fhNclusters0vsEzem   ->GetXaxis()->SetTitle("N_{clusters0}");
 412   fhmultV0vsEzem       ->GetXaxis()->SetTitle("V0 mult");
 413   fhmultFMDvsEzem      ->GetXaxis()->SetTitle("FMD mult");
 414   fhNtracksvsmultV0    ->GetXaxis()->SetTitle("N_{tracks}");
 415   fhNtrackletsvsmultV0 ->GetXaxis()->SetTitle("N_{tracklets}");
 416   fhNclusters0vsmultV0 ->GetXaxis()->SetTitle("N_{clusters0}");
 417   fhNtracksvsmultFMD   ->GetXaxis()->SetTitle("N_{tracks}");
 418   fhNtrackletsvsmultFMD->GetXaxis()->SetTitle("N_{tracklets}");
 419   fhNclusters0vsmultFMD->GetXaxis()->SetTitle("N_{clusters}");
 420   fhmultV0vsmultFMD    ->GetXaxis()->SetTitle("V0 mult");
 421
 422   fOutput->Add(fhEzdc);
 423   fOutput->Add(fhEzem);
 424   fOutput->Add(fhNtracks);
 425   fOutput->Add(fhNtracklets);
 426   fOutput->Add(fhNclusters0);
 427   fOutput->Add(fhmultV0);
 428   fOutput->Add(fhmultFMD);
 429
 430   fOutput->Add(fhEzemvsEzdc);
 431   fOutput->Add(fhNtracksvsEzdc);
 432   fOutput->Add(fhNtrackletsvsEzdc);
 433   fOutput->Add(fhNclusters0vsEzdc);
 434   fOutput->Add(fhmultV0vsEzdc);
 435   fOutput->Add(fhmultFMDvsEzdc);
 436   fOutput->Add(fhNtracksvsEzem);
 437   fOutput->Add(fhNtrackletsvsEzem);
 438   fOutput->Add(fhNclusters0vsEzem);
 439   fOutput->Add(fhmultV0vsEzem);
 440   fOutput->Add(fhmultFMDvsEzem);
 441   fOutput->Add(fhNtracksvsmultV0);
 442   fOutput->Add(fhNtrackletsvsmultV0);
 443   fOutput->Add(fhNclusters0vsmultV0);
 444   fOutput->Add(fhNtracksvsmultFMD);
 445   fOutput->Add(fhNtrackletsvsmultFMD);
 446   fOutput->Add(fhNclusters0vsmultFMD);
 447   fOutput->Add(fhmultV0vsmultFMD);
 448
 449   PostData(1, fOutput);
 450 }
 451
 452 //________________________________________________________________________
 453 void AliAnalysisTaskCentrality::UserExec(Option_t */*option*/)
 454 {
 455   // Execute analysis for current event:
 456   if(fDebug>1) printf(" **** AliAnalysisTaskCentrality::UserExec() \n");
 457
 458   if(fAnalysisInput.CompareTo("ESD")==0){
 459
 460     AliVEvent* event = InputEvent();
 461     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 462     if (!esd) return;
 463       fNev++;
 464
 465       fNTracks    = event->GetNumberOfTracks();
 466       fNPmdTracks = esd->GetNumberOfPmdTracks();
 467
 468       AliESDVZERO* esdV0 = esd->GetVZEROData();
 469       fMultV0A=esdV0->GetMTotV0A();
 470       fMultV0C=esdV0->GetMTotV0C();
 471
 472       if(fIsMCInput){
 473
 474         AliMCEvent* mcEvent = MCEvent();
 475         if (!mcEvent) {
 476           printf("   Could not retrieve MC event!!!\n");
 477           return;
 478         }
 479
 480         fNmyTracksgen = 0;
 481         AliStack *stack = 0x0; // needed for MC studies
 482         stack = MCEvent()->Stack();
 483         for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) {
 484           //get properties of mc particle
 485           AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
 486           // Primaries only
 487           if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
 488           //charged tracks only
 489           if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
 490           //same cuts as on ESDtracks
 491 //        if(TMath::Abs(mcP->Eta())>0.9)continue;
 492 //        if(mcP->Pt()<0.2)continue;
 493 //        if(mcP->Pt()>200)continue;
 494
 495           fNmyTracksgen ++;
 496         }
 497
 498         AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 499         if(!genHeader){
 500           printf("  Event generator header not available!!!\n");
 501           return;
 502         }
 503
 504         if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
 505           fbMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
 506           Int_t specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
 507           Int_t specProtonProj  = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
 508           Int_t specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn();
 509           Int_t specProtonTarg  = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp();
 510           fNpartTargMC = (Int_t)(208.-(specNeutronTarg+specProtonTarg));
 511           fNpartProjMC = (Int_t)(208.-(specNeutronProj+specProtonProj));
 512           fNNColl   = ((AliGenHijingEventHeader*) genHeader)->NN();
 513           fNNwColl  = ((AliGenHijingEventHeader*) genHeader)->NNw();
 514           fNwNColl  = ((AliGenHijingEventHeader*) genHeader)->NwN();
 515           fNwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw();
 516         }
 517
 518       }
 519
 520       fBeamEnergy = esd->GetBeamEnergy();
 521
 522       // ***** Trigger selection
 523       TString triggerClass = esd->GetFiredTriggerClasses();
 524       sprintf(fTrigClass,"%s",triggerClass.Data());
 525
 526       const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
 527       fxVertex = vertex->GetX();
 528       fyVertex = vertex->GetY();
 529       fzVertex = vertex->GetZ();
 530       if(vertex->IsFromVertexer3D()) fVertexer3d = kTRUE;
 531       else fVertexer3d = kFALSE;
 532       Double_t vertex3[3];
 533       vertex->GetXYZ(vertex3);
 534
 535       const AliMultiplicity *mult = esd->GetMultiplicity();
 536       fNTracklets = mult->GetNumberOfTracklets();
 537
 538       for(Int_t ilay=0; ilay<6; ilay++){
 539         fNClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 540       }
 541       fNSingleClusters = mult->GetNumberOfSingleClusters();
 542
 543       for(Int_t ilay=0; ilay<2; ilay++){
 544         fNChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 545       }
 546
 547
 548       AliESDFMD *fmd = esd->GetFMDData();
 549       Float_t totalMultA = 0;
 550       Float_t totalMultC = 0;
 551       const Float_t fFMDLowCut = 0.4;
 552
 553       for(UShort_t det=1;det<=3;det++) {
 554         Int_t nRings = (det==1 ? 1 : 2);
 555         for (UShort_t ir = 0; ir < nRings; ir++) {
 556           Char_t   ring = (ir == 0 ? 'I' : 'O');
 557           UShort_t nsec = (ir == 0 ? 20  : 40);
 558           UShort_t nstr = (ir == 0 ? 512 : 256);
 559           for(UShort_t sec =0; sec < nsec;  sec++)  {
 560             for(UShort_t strip = 0; strip < nstr; strip++) {
 561
 562               Float_t fmdmult = fmd->Multiplicity(det,ring,sec,strip);
 563               if(fmdmult == 0 || fmdmult == AliESDFMD::kInvalidMult) continue;
 564
 565               Float_t nParticles=0;
 566
 567                 if(fmdmult > fFMDLowCut) {
 568                   nParticles = 1.;
 569                 }
 570
 571               if (det<3) totalMultA = totalMultA + nParticles;
 572               else totalMultC = totalMultC + nParticles;
 573
 574             }
 575           }
 576         }
 577       }
 578       fMultFMDA = totalMultA;
 579       fMultFMDC = totalMultC;
 580
 581       AliESDZDC *esdZDC = esd->GetESDZDC();
 582       fESDFlag =  esdZDC->GetESDQuality();
 583       fZNCEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 584       fZPCEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 585       fZNAEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 586       fZPAEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 587       fZEM1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
 588       fZEM2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
 589
 590       fbZDC = esdZDC->GetImpactParameter();
 591       fNpartZDC = esdZDC->GetZDCParticipants();
 592       fbZDCA = esdZDC->GetImpactParamSideA();
 593       fNpartZDCA = esdZDC->GetZDCPartSideA();
 594       fbZDCC = esdZDC->GetImpactParamSideC();
 595       fNpartZDCC = esdZDC->GetZDCPartSideC();
 596
 597       const Double_t * towZNC = esdZDC->GetZN1TowerEnergy();
 598       const Double_t * towZPC = esdZDC->GetZP1TowerEnergy();
 599       const Double_t * towZNA = esdZDC->GetZN2TowerEnergy();
 600       const Double_t * towZPA = esdZDC->GetZP2TowerEnergy();
 601
 602       for(Int_t it=0; it<5; it++){
 603          fZNCtower[it] = (Float_t) (towZNC[it]);
 604          fZPCtower[it] = (Float_t) (towZPC[it]);
 605          fZNAtower[it] = (Float_t) (towZNA[it]);
 606          fZPAtower[it] = (Float_t) (towZPA[it]);
 607       }
 608
 609       Double_t xyZNC[2]={-99.,-99.}, xyZNA[2]={-99.,-99.};
 610       esdZDC->GetZNCentroidInPbPb(fBeamEnergy, xyZNC, xyZNA);
 611       for(Int_t it=0; it<2; it++){
 612          fCentrZNC[it] = xyZNC[it];
 613          fCentrZNA[it] = xyZNA[it];
 614       }
 615
 616       // filling histos
 617       fhEzdc         ->Fill((fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 618       fhEzem         ->Fill(fZEM1Energy+fZEM2Energy);
 619       fhNtracks      ->Fill(fNTracks);
 620       fhNtracklets   ->Fill(fNTracklets);
 621       fhNclusters0   ->Fill(fNClusters[0]);
 622       fhmultV0       ->Fill(fMultV0A+fMultV0C);
 623       fhmultFMD      ->Fill(fMultFMDA+fMultFMDC);
 624       fhEzemvsEzdc         ->Fill(fZEM1Energy+fZEM2Energy, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 625       fhNtracksvsEzdc      ->Fill(fNTracks, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 626       fhNtrackletsvsEzdc   ->Fill(fNTracklets,  (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 627       fhNclusters0vsEzdc   ->Fill(fNClusters[0],  (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 628       fhmultV0vsEzdc       ->Fill(fMultV0A+fMultV0C,  (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 629       fhmultFMDvsEzdc      ->Fill(fMultFMDA+fMultFMDC,  (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
 630       fhNtracksvsEzem      ->Fill(fNTracks, fZEM1Energy+fZEM2Energy);
 631       fhNtrackletsvsEzem   ->Fill(fNTracklets, fZEM1Energy+fZEM2Energy);
 632       fhNclusters0vsEzem   ->Fill(fNClusters[0], fZEM1Energy+fZEM2Energy);
 633       fhmultV0vsEzem       ->Fill(fMultV0A+fMultV0C, fZEM1Energy+fZEM2Energy);
 634       fhmultFMDvsEzem      ->Fill(fMultFMDA+fMultFMDC, fZEM1Energy+fZEM2Energy);
 635       fhNtracksvsmultV0    ->Fill(fNTracks,fMultV0A+fMultV0C);
 636       fhNtrackletsvsmultV0 ->Fill(fNTracklets,fMultV0A+fMultV0C);
 637       fhNclusters0vsmultV0 ->Fill(fNClusters[0],fMultV0A+fMultV0C);
 638       fhNtracksvsmultFMD   ->Fill(fNTracks,fMultFMDA+fMultFMDC);
 639       fhNtrackletsvsmultFMD->Fill(fNTracklets,fMultFMDA+fMultFMDC);
 640       fhNclusters0vsmultFMD->Fill(fNClusters[0],fMultFMDA+fMultFMDC);
 641       fhmultV0vsmultFMD    ->Fill(fMultV0A+fMultV0C,fMultFMDA+fMultFMDC);
 642   }
 643   else if(fAnalysisInput.CompareTo("AOD")==0){
 644     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 645     // to be implemented
 646     printf("  AOD analysis not yet implemented!!!\n\n");
 647     return;
 648   }
 649   PostData(1, fOutput);
 650 }
 651
 652 //________________________________________________________________________
 653 void AliAnalysisTaskCentrality::Terminate(Option_t */*option*/)
 654 {
 655   // Terminate analysis
 656 }