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 <TString.h>
  29 #include <TCanvas.h>
  30 #include <iostream>
  31
  32 #include "AliAnalysisManager.h"
  33 #include "AliVEvent.h"
  34 #include "AliESD.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDHeader.h"
  37 #include "AliESDInputHandler.h"
  38 #include "AliESDZDC.h"
  39 #include "AliESDFMD.h"
  40 #include "AliESDVZERO.h"
  41 #include "AliESDCentrality.h"
  42 #include "AliMultiplicity.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODMCHeader.h"
  47 #include "AliMCEvent.h"
  48 #include "AliMCEventHandler.h"
  49 #include "AliMCParticle.h"
  50 #include "AliStack.h"
  51 #include "AliHeader.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliAnalysisTaskSE.h"
  54 #include "AliGenEventHeader.h"
  55 #include "AliGenHijingEventHeader.h"
  56 #include "AliPhysicsSelectionTask.h"
  57 #include "AliPhysicsSelection.h"
  58 #include "AliBackgroundSelection.h"
  59 #include "AliCentralitySelectionTask.h"
  60
  61 ClassImp(AliCentralitySelectionTask)
  62
  63
  64 //________________________________________________________________________
  65 AliCentralitySelectionTask::AliCentralitySelectionTask():
  66 AliAnalysisTaskSE(),
  67   fDebug(0),
  68   fAnalysisInput("ESD"),
  69   fIsMCInput(kFALSE),
  70   fFile(0),
  71   fFile2(0),
  72   fCentfilename(0),
  73   fCentfilename2(0),
  74   fCentV0M(0),
  75   fCentFMD(0),
  76   fCentTRK(0),
  77   fCentTKL(0),
  78   fCentCL0(0),
  79   fCentV0MvsFMD(0),
  80   fCentTKLvsV0M(0),
  81   fCentZEMvsZDC(0),
  82   fHtempV0M(0),
  83   fHtempFMD(0),
  84   fHtempTRK(0),
  85   fHtempTKL(0),
  86   fHtempCL0(0),
  87   fHtempV0MvsFMD(0),
  88   fHtempTKLvsV0M(0),
  89   fHtempZEMvsZDC(0)
  90 {
  91   // Default constructor
  92 }
  93
  94 //________________________________________________________________________
  95 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
  96   AliAnalysisTaskSE(name),
  97   fDebug(0),
  98   fAnalysisInput("ESD"),
  99   fIsMCInput(kFALSE),
 100   fFile(0),
 101   fFile2(0),
 102   fCentfilename(0),
 103   fCentfilename2(0),
 104   fCentV0M(0),
 105   fCentFMD(0),
 106   fCentTRK(0),
 107   fCentTKL(0),
 108   fCentCL0(0),
 109   fCentV0MvsFMD(0),
 110   fCentTKLvsV0M(0),
 111   fCentZEMvsZDC(0),
 112   fHtempV0M(0),
 113   fHtempFMD(0),
 114   fHtempTRK(0),
 115   fHtempTKL(0),
 116   fHtempCL0(0),
 117   fHtempV0MvsFMD(0),
 118   fHtempTKLvsV0M(0),
 119   fHtempZEMvsZDC(0)
 120 {
 121   // Default constructor
 122 }
 123
 124 //________________________________________________________________________
 125 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 126 {
 127   // Assignment operator
 128   if (this!=&c) {
 129     AliAnalysisTaskSE::operator=(c);
 130   }
 131   return *this;
 132 }
 133
 134 //________________________________________________________________________
 135 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 136   AliAnalysisTaskSE(ana),
 137   fDebug(ana.fDebug),
 138   fAnalysisInput(ana.fDebug),
 139   fIsMCInput(ana.fIsMCInput),
 140   fFile(ana.fFile),
 141   fFile2(ana.fFile2),
 142   fCentfilename(ana.fCentfilename),
 143   fCentfilename2(ana.fCentfilename2),
 144   fCentV0M(ana.fCentV0M),
 145   fCentFMD(ana.fCentFMD),
 146   fCentTRK(ana.fCentTRK),
 147   fCentTKL(ana.fCentTKL),
 148   fCentCL0(ana.fCentCL0),
 149   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 150   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 151   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 152   fHtempV0M(ana.fHtempV0M),
 153   fHtempFMD(ana.fHtempFMD),
 154   fHtempTRK(ana.fHtempTRK),
 155   fHtempTKL(ana.fHtempTKL),
 156   fHtempCL0(ana.fHtempCL0),
 157   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 158   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 159   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
 160 {
 161   // Copy Constructor
 162 }
 163
 164 //________________________________________________________________________
 165  AliCentralitySelectionTask::~AliCentralitySelectionTask()
 166  {
 167    // Destructor
 168  }
 169
 170 //________________________________________________________________________
 171 void AliCentralitySelectionTask::UserCreateOutputObjects()
 172 {
 173   // Create the output containers
 174   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 175 }
 176
 177 //________________________________________________________________________
 178 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 179 {
 180   // Execute analysis for current event:
 181   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 182
 183   Float_t  zncEnergy;               //  ZNC Energy
 184   Float_t  zpcEnergy;               //  ZPC Energy
 185   Float_t  znaEnergy;               //  ZNA Energy
 186   Float_t  zpaEnergy;               //  ZPA Energy
 187   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 188   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 189
 190   Int_t    nTracks    = 0;          //  no. tracks
 191   Int_t    nTracklets = 0;          //  no. tracklets
 192   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
 193   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 194
 195   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 196   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 197   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 198   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 199
 200   AliESDCentrality *esdCent = 0;
 201
 202   if(fAnalysisInput.CompareTo("ESD")==0){
 203
 204     AliVEvent* event = InputEvent();
 205     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 206
 207     esdCent = esd->GetCentrality();
 208
 209     // ***** V0 info
 210     AliESDVZERO* esdV0 = esd->GetVZEROData();
 211     multV0A=esdV0->GetMTotV0A();
 212     multV0C=esdV0->GetMTotV0C();
 213
 214     // ***** CB info (tracklets, clusters, chips)
 215     nTracks    = event->GetNumberOfTracks();
 216
 217     const AliMultiplicity *mult = esd->GetMultiplicity();
 218
 219     nTracklets = mult->GetNumberOfTracklets();
 220
 221     for(Int_t ilay=0; ilay<6; ilay++){
 222       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 223     }
 224
 225     for(Int_t ilay=0; ilay<2; ilay++){
 226       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 227     }
 228
 229
 230     // ***** FMD info
 231     AliESDFMD *fmd = esd->GetFMDData();
 232     Float_t totalMultA = 0;
 233     Float_t totalMultC = 0;
 234     const Float_t fFMDLowCut = 0.4;
 235
 236     for(UShort_t det=1;det<=3;det++) {
 237       Int_t nRings = (det==1 ? 1 : 2);
 238       for (UShort_t ir = 0; ir < nRings; ir++) {
 239         Char_t   ring = (ir == 0 ? 'I' : 'O');
 240         UShort_t nsec = (ir == 0 ? 20  : 40);
 241         UShort_t nstr = (ir == 0 ? 512 : 256);
 242         for(UShort_t sec =0; sec < nsec;  sec++)  {
 243           for(UShort_t strip = 0; strip < nstr; strip++) {
 244
 245             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
 246             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
 247
 248             Float_t nParticles=0;
 249
 250             if(FMDmult > fFMDLowCut) {
 251               nParticles = 1.;
 252             }
 253
 254             if (det<3) totalMultA = totalMultA + nParticles;
 255             else totalMultC = totalMultC + nParticles;
 256
 257           }
 258         }
 259       }
 260     }
 261     multFMDA = totalMultA;
 262     multFMDC = totalMultC;
 263
 264     // ***** ZDC info
 265     AliESDZDC *esdZDC = esd->GetESDZDC();
 266     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 267     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 268     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 269     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 270     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
 271     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
 272
 273   }
 274   else if(fAnalysisInput.CompareTo("AOD")==0){
 275     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 276     // to be implemented
 277     printf("  AOD analysis not yet implemented!!!\n\n");
 278     return;
 279   }
 280
 281   // ***** Centrality Selection
 282   fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
 283   fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 284   fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 285   fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 286   fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 287
 288   fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 289   fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 290   fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
 291
 292   esdCent->SetCentralityV0M(fCentV0M);
 293   esdCent->SetCentralityFMD(fCentFMD);
 294   esdCent->SetCentralityTRK(fCentTRK);
 295   esdCent->SetCentralityTKL(fCentTKL);
 296   esdCent->SetCentralityCL0(fCentCL0);
 297   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 298   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 299   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 300 }
 301
 302 //________________________________________________________________________
 303 void AliCentralitySelectionTask::ReadCentralityHistos()
 304 {
 305   fFile  = new TFile(fCentfilename);
 306   fHtempV0M  = (TH1D*) (fFile->Get("hmultV0_percentile"));
 307   fHtempFMD  = (TH1D*) (fFile->Get("hmultFMD_percentile"));
 308   fHtempTRK  = (TH1D*) (fFile->Get("hNtracks_percentile"));
 309   fHtempTKL  = (TH1D*) (fFile->Get("hNtracklets_percentile"));
 310   fHtempCL0  = (TH1D*) (fFile->Get("hNclusters0_percentile"));
 311 }
 312
 313 void AliCentralitySelectionTask::ReadCentralityHistos2()
 314 {
 315   fFile2  = new TFile(fCentfilename2);
 316   fHtempV0MvsFMD  = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 317   fHtempTKLvsV0M  = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 318   fHtempZEMvsZDC  = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 319 }
 320
 321 void AliCentralitySelectionTask::SetPercentileFile(TString filename)
 322 {
 323   fCentfilename = filename;
 324   ReadCentralityHistos();
 325 }
 326
 327 void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
 328 {
 329   fCentfilename2 = filename;
 330   ReadCentralityHistos2();
 331 }
 332
 333 //________________________________________________________________________
 334 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 335 {
 336   // Terminate analysis
 337   fFile->Close();
 338   fFile2->Close();
 339 }
 340 //________________________________________________________________________
 341