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