PWG/EMCAL/AliAnalysisTaskEmcal.cxx

   1 // $Id: AliAnalysisTaskEmcal.cxx 56756 2012-05-30 05:03:02Z loizides $
   2 //
   3 // Emcal base analysis task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include "AliAnalysisTaskEmcal.h"
   8
   9 #include <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TList.h>
  12 #include <TObject.h>
  13 #include <TH1F.h>
  14
  15 #include "AliAODEvent.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliESDEvent.h"
  20 #include "AliEmcalParticle.h"
  21 #include "AliEventplane.h"
  22 #include "AliInputEventHandler.h"
  23 #include "AliLog.h"
  24 #include "AliMCParticle.h"
  25 #include "AliVCluster.h"
  26 #include "AliVEventHandler.h"
  27 #include "AliVParticle.h"
  28
  29 ClassImp(AliAnalysisTaskEmcal)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
  33   AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
  34   fAnaType(kTPC),
  35   fForceBeamType(kNA),
  36   fGeneralHistograms(kFALSE),
  37   fInitialized(kFALSE),
  38   fCreateHisto(kTRUE),
  39   fTracksName(),
  40   fCaloName(),
  41   fMinCent(-999),
  42   fMaxCent(-999),
  43   fMinVz(-999),
  44   fMaxVz(-999),
  45   fOffTrigger(AliVEvent::kAny),
  46   fTrigClass(),
  47   fNbins(500),
  48   fMinBinPt(0),
  49   fMaxBinPt(250),
  50   fClusPtCut(0.15),
  51   fTrackPtCut(0.15),
  52   fTrackMinEta(-0.9),
  53   fTrackMaxEta(0.9),
  54   fTrackMinPhi(-10),
  55   fTrackMaxPhi(10),
  56   fClusTimeCutLow(-10),
  57   fClusTimeCutUp(10),
  58   fGeom(0),
  59   fTracks(0),
  60   fCaloClusters(0),
  61   fCent(0),
  62   fCentBin(-1),
  63   fEPV0(-1.0),
  64   fEPV0A(-1.0),
  65   fEPV0C(-1.0),
  66   fNVertCont(0),
  67   fBeamType(kNA),
  68   fOutput(0)
  69 {
  70   // Default constructor.
  71
  72   fVertex[0] = 0;
  73   fVertex[1] = 0;
  74   fVertex[2] = 0;
  75 }
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
  79   AliAnalysisTaskSE(name),
  80   fAnaType(kTPC),
  81   fForceBeamType(kNA),
  82   fGeneralHistograms(kFALSE),
  83   fInitialized(kFALSE),
  84   fCreateHisto(histo),
  85   fTracksName(),
  86   fCaloName(),
  87   fMinCent(-999),
  88   fMaxCent(-999),
  89   fMinVz(-999),
  90   fMaxVz(-999),
  91   fOffTrigger(AliVEvent::kAny),
  92   fTrigClass(),
  93   fNbins(500),
  94   fMinBinPt(0),
  95   fMaxBinPt(250),
  96   fClusPtCut(0.15),
  97   fTrackPtCut(0.15),
  98   fTrackMinEta(-0.9),
  99   fTrackMaxEta(0.9),
 100   fTrackMinPhi(-10),
 101   fTrackMaxPhi(10),
 102   fClusTimeCutLow(-10),
 103   fClusTimeCutUp(10),
 104   fGeom(0),
 105   fTracks(0),
 106   fCaloClusters(0),
 107   fCent(0),
 108   fCentBin(-1),
 109   fEPV0(-1.0),
 110   fEPV0A(-1.0),
 111   fEPV0C(-1.0),
 112   fNVertCont(0),
 113   fBeamType(kNA),
 114   fOutput(0)
 115 {
 116   // Standard constructor.
 117
 118   fVertex[0] = 0;
 119   fVertex[1] = 0;
 120   fVertex[2] = 0;
 121
 122   if (fCreateHisto) {
 123     DefineOutput(1, TList::Class());
 124   }
 125 }
 126
 127 //________________________________________________________________________
 128 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
 129 {
 130   // Destructor
 131 }
 132
 133 //________________________________________________________________________
 134 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
 135 {
 136   // Create user output.
 137   if (!fCreateHisto)
 138     return;
 139
 140   OpenFile(1);
 141   fOutput = new TList();
 142   fOutput->SetOwner();
 143
 144   if (!fGeneralHistograms)
 145     return;
 146
 147   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
 148   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 149   fHistCentrality->GetYaxis()->SetTitle("counts");
 150   fOutput->Add(fHistCentrality);
 151
 152   fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
 153   fHistZVertex->GetXaxis()->SetTitle("z");
 154   fHistZVertex->GetYaxis()->SetTitle("counts");
 155   fOutput->Add(fHistZVertex);
 156
 157   PostData(1, fOutput);
 158 }
 159
 160 //________________________________________________________________________
 161 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
 162 {
 163   fHistCentrality->Fill(fCent);
 164   fHistZVertex->Fill(fVertex[2]);
 165
 166   return kTRUE;
 167 }
 168
 169 //________________________________________________________________________
 170 void AliAnalysisTaskEmcal::UserExec(Option_t *)
 171 {
 172   // Main loop, called for each event.
 173
 174   if (!fInitialized)
 175     ExecOnce();
 176
 177   if (!fInitialized)
 178     return;
 179
 180   if (!RetrieveEventObjects())
 181     return;
 182
 183   if (!IsEventSelected())
 184     return;
 185
 186   if (fGeneralHistograms && fCreateHisto) {
 187     if (!FillGeneralHistograms())
 188       return;
 189   }
 190
 191   if (!Run())
 192     return;
 193
 194   if (fCreateHisto) {
 195     if (!FillHistograms())
 196       return;
 197   }
 198
 199   if (fCreateHisto && fOutput) {
 200     // information for this iteration of the UserExec in the container
 201     PostData(1, fOutput);
 202   }
 203 }
 204
 205 //________________________________________________________________________
 206 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Bool_t acceptMC) const
 207 {
 208   // Return true if cluster is accepted.
 209
 210   if (!clus)
 211     return kFALSE;
 212
 213   if (!clus->IsEMCAL())
 214     return kFALSE;
 215
 216   if (!acceptMC && clus->Chi2() == 100)
 217     return kFALSE;
 218
 219   if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow)
 220     return kFALSE;
 221
 222   TLorentzVector nPart;
 223   clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 224
 225   if (nPart.Et() < fClusPtCut)
 226     return kFALSE;
 227
 228   return kTRUE;
 229 }
 230
 231 //________________________________________________________________________
 232 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVTrack *track, Bool_t acceptMC) const
 233 {
 234   // Return true if track is accepted.
 235
 236   if (!track)
 237     return kFALSE;
 238
 239   if (!acceptMC && track->GetLabel() == 100)
 240     return kFALSE;
 241
 242   if (track->Pt() < fTrackPtCut)
 243     return kFALSE;
 244
 245   if (track->Eta() < fTrackMinEta || track->Eta() > fTrackMaxEta ||
 246       track->Phi() < fTrackMinPhi || track->Phi() > fTrackMaxPhi)
 247     return kFALSE;
 248
 249   return kTRUE;
 250 }
 251
 252 //________________________________________________________________________
 253 Bool_t AliAnalysisTaskEmcal::AcceptEmcalPart(AliEmcalParticle *part, Bool_t acceptMC) const
 254 {
 255   // Return true if EMCal particle is accepted.
 256
 257   if (!part)
 258     return kFALSE;
 259
 260   if (part->IsTrack()) {
 261     if (part->Pt() < fTrackPtCut)
 262       return kFALSE;
 263
 264     if (part->Eta() < fTrackMinEta || part->Eta() > fTrackMaxEta ||
 265         part->Phi() < fTrackMinPhi || part->Phi() > fTrackMaxPhi)
 266       return kFALSE;
 267   }
 268
 269   if (part->IsCluster()) {
 270     if (!part->IsEMCAL())
 271       return kFALSE;
 272
 273     if (part->Pt() < fClusPtCut)
 274       return kFALSE;
 275   }
 276
 277   if (!acceptMC && part->IsMC())
 278     return kFALSE;
 279
 280   return kTRUE;
 281 }
 282
 283 //________________________________________________________________________
 284 void AliAnalysisTaskEmcal::ExecOnce()
 285 {
 286   // Init the analysis.
 287
 288   if (!InputEvent()) {
 289     AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
 290     return;
 291   }
 292
 293   fGeom = AliEMCALGeometry::GetInstance();
 294   if (!fGeom)
 295     AliWarning(Form("%s: Can not create geometry", GetName()));
 296
 297   if (!fCaloName.IsNull() && !fCaloClusters) {
 298     fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 299     if (!fCaloClusters) {
 300       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCaloName.Data()));
 301       return;
 302     } else {
 303       TClass *cl = fCaloClusters->GetClass();
 304       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
 305         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCaloName.Data()));
 306         fCaloClusters = 0;
 307         return;
 308       }
 309     }
 310   }
 311
 312   if (!fTracksName.IsNull() && !fTracks) {
 313     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 314     if (!fTracks) {
 315       AliError(Form("%s: Could not retrieve tracks %s!", GetName(), fTracksName.Data()));
 316       return;
 317     } else {
 318       TClass *cl = fTracks->GetClass();
 319       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 320         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracksName.Data()));
 321         fTracks = 0;
 322         return;
 323       }
 324     }
 325   }
 326
 327   fInitialized = kTRUE;
 328 }
 329
 330 //_____________________________________________________
 331 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
 332 {
 333   // Get beam type : pp-AA-pA
 334   // ESDs have it directly, AODs get it from hardcoded run number ranges
 335
 336   if (fForceBeamType != kNA)
 337     return fForceBeamType;
 338
 339   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
 340   if (esd) {
 341     const AliESDRun *run = esd->GetESDRun();
 342     TString beamType = run->GetBeamType();
 343     if (beamType == "p-p")
 344       return kpp;
 345     else if (beamType == "A-A")
 346       return kAA;
 347     else if (beamType == "p-A")
 348       return kpA;
 349     else
 350       return kNA;
 351   } else {
 352     Int_t runNumber = InputEvent()->GetRunNumber();
 353     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
 354         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
 355     {
 356       return kAA;
 357     } else {
 358       return kpp;
 359     }
 360   }
 361 }
 362
 363 //________________________________________________________________________
 364 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
 365 {
 366   // Check if event is selected
 367
 368   if (fOffTrigger != AliVEvent::kAny) {
 369     UInt_t res = 0;
 370     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 371     if (eev) {
 372       res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 373     } else {
 374       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 375       if (aev) {
 376         res = aev->GetHeader()->GetOfflineTrigger();
 377       }
 378     }
 379     if ((res & fOffTrigger) == 0)
 380       return kFALSE;
 381   }
 382
 383   if (!fTrigClass.IsNull()) {
 384     TString fired;
 385     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
 386     if (eev) {
 387       fired = eev->GetFiredTriggerClasses();
 388     } else {
 389       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
 390       if (aev) {
 391         fired = aev->GetFiredTriggerClasses();
 392       }
 393     }
 394     if (!fired.Contains("-B-"))
 395       return kFALSE;
 396     TObjArray *arr = fTrigClass.Tokenize("|");
 397     if (!arr)
 398       return kFALSE;
 399     Bool_t match = 0;
 400     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
 401       TObject *obj = arr->At(i);
 402       if (!obj)
 403         continue;
 404       if (fired.Contains(obj->GetName())) {
 405         match = 1;
 406         break;
 407       }
 408     }
 409     delete arr;
 410     if (!match)
 411       return kFALSE;
 412   }
 413
 414   if ((fMinCent != -999) && (fMaxCent != -999)) {
 415     if (fCent<fMinCent)
 416       return kFALSE;
 417     if (fCent>fMaxCent)
 418       return kFALSE;
 419   }
 420
 421   if ((fMinVz != -999) && (fMaxVz != -999)) {
 422     if (fNVertCont == 0 )
 423       return kFALSE;
 424     Double_t vz = fVertex[2];
 425     if (vz<fMinVz)
 426       return kFALSE;
 427     if (vz>fMaxVz)
 428       return kFALSE;
 429   }
 430
 431   return kTRUE;
 432 }
 433
 434 //________________________________________________________________________
 435 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
 436 {
 437   // Get array from event.
 438
 439   TClonesArray *arr = 0;
 440   TString sname(name);
 441   if (!sname.IsNull()) {
 442     arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
 443     if (!arr) {
 444       AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
 445       return 0;
 446     }
 447   } else {
 448     return 0;
 449   }
 450
 451   if (!clname)
 452     return arr;
 453
 454   TString objname(arr->GetClass()->GetName());
 455   TClass cls(objname);
 456   if (!cls.InheritsFrom(clname)) {
 457     AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
 458                     GetName(), cls.GetName(), name, clname));
 459     return 0;
 460   }
 461   return arr;
 462 }
 463
 464 //________________________________________________________________________
 465 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
 466 {
 467   // Retrieve objects from event.
 468
 469   fVertex[0] = 0;
 470   fVertex[1] = 0;
 471   fVertex[2] = 0;
 472   fNVertCont = 0;
 473
 474   const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
 475   if (vert) {
 476     vert->GetXYZ(fVertex);
 477     fNVertCont = vert->GetNContributors();
 478   }
 479
 480   fBeamType = GetBeamType();
 481
 482   if (fBeamType == kAA) {
 483     AliCentrality *aliCent = InputEvent()->GetCentrality();
 484     if (aliCent) {
 485       fCent = aliCent->GetCentralityPercentile("V0M");
 486       if      (fCent >=  0 && fCent <   10) fCentBin = 0;
 487       else if (fCent >= 10 && fCent <   30) fCentBin = 1;
 488       else if (fCent >= 30 && fCent <   50) fCentBin = 2;
 489       else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
 490       else {
 491         AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
 492         fCentBin = 3;
 493       }
 494     } else {
 495       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
 496       fCentBin = 3;
 497     }
 498     AliEventplane *aliEP = InputEvent()->GetEventplane();
 499     if (aliEP) {
 500       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
 501       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
 502       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
 503     } else {
 504       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
 505     }
 506   } else {
 507     fCent = 99;
 508     fCentBin = 0;
 509   }
 510
 511   return kTRUE;
 512 }