ANALYSIS/AliCentralitySelectionTask.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 //   Class AliCentralitySelectionTask
  18 //   Class to analyze determine centrality
  19 //   author: Alberica Toia
  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 <TObjString.h>
  29 #include <TString.h>
  30 #include <TCanvas.h>
  31 #include <TROOT.h>
  32 #include <TDirectory.h>
  33 #include <iostream>
  34
  35 #include "AliAnalysisManager.h"
  36 #include "AliVEvent.h"
  37 #include "AliESD.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDHeader.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDZDC.h"
  42 #include "AliESDFMD.h"
  43 #include "AliESDVZERO.h"
  44 #include "AliESDCentrality.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliAODHandler.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODVertex.h"
  49 #include "AliAODMCHeader.h"
  50 #include "AliMCEvent.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCParticle.h"
  53 #include "AliStack.h"
  54 #include "AliHeader.h"
  55 #include "AliAODMCParticle.h"
  56 #include "AliAnalysisTaskSE.h"
  57 #include "AliGenEventHeader.h"
  58 #include "AliGenHijingEventHeader.h"
  59 #include "AliPhysicsSelectionTask.h"
  60 #include "AliPhysicsSelection.h"
  61 #include "AliBackgroundSelection.h"
  62 #include "AliCentralitySelectionTask.h"
  63
  64 ClassImp(AliCentralitySelectionTask)
  65
  66
  67 //________________________________________________________________________
  68 AliCentralitySelectionTask::AliCentralitySelectionTask():
  69 AliAnalysisTaskSE(),
  70   fDebug(0),
  71   fAnalysisInput("ESD"),
  72   fIsMCInput(kFALSE),
  73   fFile(0),
  74   fFile2(0),
  75   fCentfilename(""),
  76   fCentfilename2(""),
  77   fFileList(new TList),
  78   fFileList2(new TList),
  79   fCurrentRun(-1),
  80   fCentV0M(0),
  81   fCentFMD(0),
  82   fCentTRK(0),
  83   fCentTKL(0),
  84   fCentCL0(0),
  85   fCentCL1(0),
  86   fCentV0MvsFMD(0),
  87   fCentTKLvsV0M(0),
  88   fCentZEMvsZDC(0),
  89   fHtempV0M(0),
  90   fHtempFMD(0),
  91   fHtempTRK(0),
  92   fHtempTKL(0),
  93   fHtempCL0(0),
  94   fHtempCL1(0),
  95   fHtempV0MvsFMD(0),
  96   fHtempTKLvsV0M(0),
  97   fHtempZEMvsZDC(0)
  98 {
  99   // Default constructor
 100   fFileList->SetOwner();
 101   fFileList2->SetOwner();
 102
 103   AliInfo("Centrality Selection enabled.");
 104 }
 105
 106 //________________________________________________________________________
 107 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
 108   AliAnalysisTaskSE(name),
 109   fDebug(0),
 110   fAnalysisInput("ESD"),
 111   fIsMCInput(kFALSE),
 112   fFile(0),
 113   fFile2(0),
 114   fCentfilename(""),
 115   fCentfilename2(""),
 116   fFileList(new TList),
 117   fFileList2(new TList),
 118   fCurrentRun(-1),
 119   fCentV0M(0),
 120   fCentFMD(0),
 121   fCentTRK(0),
 122   fCentTKL(0),
 123   fCentCL0(0),
 124   fCentCL1(0),
 125   fCentV0MvsFMD(0),
 126   fCentTKLvsV0M(0),
 127   fCentZEMvsZDC(0),
 128   fHtempV0M(0),
 129   fHtempFMD(0),
 130   fHtempTRK(0),
 131   fHtempTKL(0),
 132   fHtempCL0(0),
 133   fHtempCL1(0),
 134   fHtempV0MvsFMD(0),
 135   fHtempTKLvsV0M(0),
 136   fHtempZEMvsZDC(0)
 137 {
 138   // Default constructor
 139   fFileList->SetOwner();
 140   fFileList2->SetOwner();
 141
 142   AliInfo("Centrality Selection enabled.");
 143 }
 144
 145 //________________________________________________________________________
 146 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 147 {
 148   // Assignment operator
 149   if (this!=&c) {
 150     AliAnalysisTaskSE::operator=(c);
 151   }
 152   return *this;
 153 }
 154
 155 //________________________________________________________________________
 156 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 157   AliAnalysisTaskSE(ana),
 158   fDebug(ana.fDebug),
 159   fAnalysisInput(ana.fDebug),
 160   fIsMCInput(ana.fIsMCInput),
 161   fFile(ana.fFile),
 162   fFile2(ana.fFile2),
 163   fCentfilename(ana.fCentfilename),
 164   fCentfilename2(ana.fCentfilename2),
 165   fFileList(ana.fFileList),
 166   fFileList2(ana.fFileList2),
 167   fCurrentRun(ana.fCurrentRun),
 168   fCentV0M(ana.fCentV0M),
 169   fCentFMD(ana.fCentFMD),
 170   fCentTRK(ana.fCentTRK),
 171   fCentTKL(ana.fCentTKL),
 172   fCentCL0(ana.fCentCL0),
 173   fCentCL1(ana.fCentCL1),
 174   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 175   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 176   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 177   fHtempV0M(ana.fHtempV0M),
 178   fHtempFMD(ana.fHtempFMD),
 179   fHtempTRK(ana.fHtempTRK),
 180   fHtempTKL(ana.fHtempTKL),
 181   fHtempCL0(ana.fHtempCL0),
 182   fHtempCL1(ana.fHtempCL1),
 183   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 184   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 185   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
 186 {
 187   // Copy Constructor
 188 }
 189
 190 //________________________________________________________________________
 191 AliCentralitySelectionTask::~AliCentralitySelectionTask()
 192 {
 193   // Destructor
 194
 195   if (fFileList) {
 196     fFileList->Clear();
 197     delete fFileList;
 198   }
 199   fFileList = NULL;
 200
 201   if (fFileList2) {
 202     fFileList2->Clear();
 203     delete fFileList2;
 204   }
 205   fFileList2 = NULL;
 206 }
 207
 208 //________________________________________________________________________
 209 void AliCentralitySelectionTask::UserCreateOutputObjects()
 210 {
 211   // Create the output containers
 212   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 213
 214   if (fFileList->GetEntries() < 1) {
 215     AliError("No Inputfiles Added");
 216   }
 217
 218   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
 219 }
 220
 221 //________________________________________________________________________
 222 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 223 {
 224   // Execute analysis for current event:
 225   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 226
 227   Float_t  zncEnergy;          //  ZNC Energy
 228   Float_t  zpcEnergy;          //  ZPC Energy
 229   Float_t  znaEnergy;          //  ZNA Energy
 230   Float_t  zpaEnergy;          //  ZPA Energy
 231   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 232   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 233
 234   Int_t    nTracks = 0;             //  no. tracks
 235   Int_t    nTracklets = 0;          //  no. tracklets
 236   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
 237   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 238
 239   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 240   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 241   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 242   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 243
 244   AliESDCentrality *esdCent = 0;
 245
 246   if(fAnalysisInput.CompareTo("ESD")==0){
 247
 248     AliVEvent* event = InputEvent();
 249     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 250
 251     if (SetupRun(esd))
 252       return;
 253
 254     esdCent = esd->GetCentrality();
 255
 256     // ***** V0 info
 257     AliESDVZERO* esdV0 = esd->GetVZEROData();
 258     multV0A=esdV0->GetMTotV0A();
 259     multV0C=esdV0->GetMTotV0C();
 260
 261     // ***** CB info (tracklets, clusters, chips)
 262     nTracks    = event->GetNumberOfTracks();
 263
 264     const AliMultiplicity *mult = esd->GetMultiplicity();
 265
 266     nTracklets = mult->GetNumberOfTracklets();
 267
 268     for(Int_t ilay=0; ilay<6; ilay++){
 269       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 270     }
 271
 272     for(Int_t ilay=0; ilay<2; ilay++){
 273       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 274     }
 275
 276
 277     // ***** FMD info
 278     AliESDFMD *fmd = esd->GetFMDData();
 279     Float_t totalMultA = 0;
 280     Float_t totalMultC = 0;
 281     const Float_t fFMDLowCut = 0.4;
 282
 283     for(UShort_t det=1;det<=3;det++) {
 284       Int_t nRings = (det==1 ? 1 : 2);
 285       for (UShort_t ir = 0; ir < nRings; ir++) {
 286         Char_t   ring = (ir == 0 ? 'I' : 'O');
 287         UShort_t nsec = (ir == 0 ? 20  : 40);
 288         UShort_t nstr = (ir == 0 ? 512 : 256);
 289         for(UShort_t sec =0; sec < nsec;  sec++)  {
 290           for(UShort_t strip = 0; strip < nstr; strip++) {
 291
 292             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
 293             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
 294
 295             Float_t nParticles=0;
 296
 297             if(FMDmult > fFMDLowCut) {
 298               nParticles = 1.;
 299             }
 300
 301             if (det<3) totalMultA = totalMultA + nParticles;
 302             else totalMultC = totalMultC + nParticles;
 303
 304           }
 305         }
 306       }
 307     }
 308     multFMDA = totalMultA;
 309     multFMDC = totalMultC;
 310
 311     // ***** ZDC info
 312     AliESDZDC *esdZDC = esd->GetESDZDC();
 313     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 314     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 315     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 316     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 317     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
 318     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
 319
 320   }
 321   else if(fAnalysisInput.CompareTo("AOD")==0){
 322     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 323     // to be implemented
 324     printf("  AOD analysis not yet implemented!!!\n\n");
 325     return;
 326   }
 327
 328   // ***** Centrality Selection
 329   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
 330   ///  else     printf("  Centrality by V0 not available!!!\n\n");
 331   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 332   //  else     printf("  Centrality by FMD not available!!!\n\n");
 333   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 334   //  else     printf("  Centrality by TRK not available!!!\n\n");
 335   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 336   //  else     printf("  Centrality by TKL not available!!!\n\n");
 337   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 338   //  else     printf("  Centrality by CL0 not available!!!\n\n");
 339   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
 340   ///  else     printf("  Centrality by CL1 not available!!!\n\n");
 341
 342   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 343   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");
 344   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 345   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");
 346   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
 347   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");
 348
 349   esdCent->SetCentralityV0M(fCentV0M);
 350   esdCent->SetCentralityFMD(fCentFMD);
 351   esdCent->SetCentralityTRK(fCentTRK);
 352   esdCent->SetCentralityTKL(fCentTKL);
 353   esdCent->SetCentralityCL0(fCentCL0);
 354   esdCent->SetCentralityCL1(fCentCL1);
 355   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 356   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 357   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 358 }
 359
 360 //________________________________________________________________________
 361 void AliCentralitySelectionTask::ReadCentralityHistos()
 362 {
 363 //  Read centrality histograms
 364     TDirectory *owd = gDirectory;
 365     fFile  = TFile::Open(fCentfilename);
 366     owd->cd();
 367     fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 368     fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 369     fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 370     fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 371     fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 372     fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 373 }
 374
 375 //________________________________________________________________________
 376 void AliCentralitySelectionTask::ReadCentralityHistos2()
 377 {
 378 //  Read centrality histograms
 379     TDirectory *owd = gDirectory;
 380     fFile2  = TFile::Open(fCentfilename2);
 381     owd->cd();
 382     fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 383     fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 384     fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 385 }
 386
 387 //________________________________________________________________________
 388 void AliCentralitySelectionTask::SetPercentileFile(TString filename)
 389 {
 390 // Set the percentile file name
 391   fCentfilename = filename;
 392   ReadCentralityHistos();
 393 }
 394
 395 //________________________________________________________________________
 396 void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
 397 {
 398 // Set the percentile file name
 399   fCentfilename2 = filename;
 400   ReadCentralityHistos2();
 401 }
 402
 403 //________________________________________________________________________
 404 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 405 {
 406   // Terminate analysis
 407   if (fFile && fFile->IsOpen())
 408     fFile->Close();
 409   if (fFile2 && fFile2->IsOpen())
 410     fFile2->Close();
 411 }
 412 //________________________________________________________________________
 413 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
 414 {
 415   // Setup files for run
 416
 417   if (!esd)
 418     return -1;
 419
 420   // check if something to be done
 421   if (fCurrentRun == esd->GetRunNumber())
 422     return 0;
 423   else
 424     fCurrentRun = esd->GetRunNumber();
 425
 426   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 427
 428   Int_t runNo = fCurrentRun;
 429
 430   // CHANGE HERE FOR RUN RANGES
 431   if ( runNo == 137162 ) runNo = 137161;
 432   else if ( runNo == 137365 ) runNo = 137366;
 433   // CHANGE HERE FOR RUN RANGES
 434
 435   TString runName(Form("%d", runNo));
 436   TString fileName("");
 437   Bool_t isRunKnown = kFALSE;
 438
 439   // Check if run is in fileList
 440   // if not, take the last name in the list
 441   for ( Int_t idx=0 ; idx < fFileList->GetEntries(); ++idx ) {
 442
 443     TString str((dynamic_cast<TObjString*>(fFileList->At(idx)))->GetString());
 444     if (str.Contains(runName)) {
 445       fileName += str;
 446       isRunKnown = kTRUE;
 447       break;
 448     }
 449   }
 450
 451   if (!isRunKnown) {
 452     if (fFileList->Last()) {
 453       fileName += (dynamic_cast<TObjString*>(fFileList->Last()))->GetString();
 454       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 455     }
 456   }
 457
 458   if (fileName.Contains(".root")) {
 459     AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
 460     SetPercentileFile(fileName.Data());
 461     return 0;
 462   }
 463
 464   return -1;
 465 }