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