PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx

   1 // $Id$
   2
   3 #include "AliAnalysisTaskEMCALPi0PbPb.h"
   4 #include <TAxis.h>
   5 #include <TChain.h>
   6 #include <TClonesArray.h>
   7 #include <TFile.h>
   8 #include <TGeoGlobalMagField.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TNtuple.h>
  14 #include <TProfile.h>
  15 #include <TString.h>
  16 #include <TVector2.h>
  17 #include "AliAODEvent.h"
  18 #include "AliAODVertex.h"
  19 #include "AliAnalysisManager.h"
  20 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  21 #include "AliCDBManager.h"
  22 #include "AliCentrality.h"
  23 #include "AliEMCALGeoUtils.h"
  24 #include "AliEMCALGeometry.h"
  25 #include "AliEMCALRecoUtils.h"
  26 #include "AliESDEvent.h"
  27 #include "AliESDVertex.h"
  28 #include "AliESDtrack.h"
  29 #include "AliESDtrackCuts.h"
  30 #include "AliGeomManager.h"
  31 #include "AliInputEventHandler.h"
  32 #include "AliLog.h"
  33 #include "AliMagF.h"
  34 #include "AliTrackerBase.h"
  35
  36 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  40   : AliAnalysisTaskSE(name),
  41     fCentVar("V0M"),
  42     fCentFrom(0),
  43     fCentTo(100),
  44     fVtxZMin(-10),
  45     fVtxZMax(+10),
  46     fUseQualFlag(1),
  47     fClusName(),
  48     fDoNtuple(0),
  49     fDoAfterburner(0),
  50     fAsymMax(1),
  51     fNminCells(1),
  52     fMinE(0.100),
  53     fMinErat(0),
  54     fMinEcc(0),
  55     fGeoName("EMCAL_FIRSTYEARV1"),
  56     fMinNClustPerTrack(50),
  57     fMinPtPerTrack(1.0),
  58     fIsoDist(0.2),
  59     fTrClassNames(""),
  60     fTrCuts(0),
  61     fNEvs(0),
  62     fGeom(0),
  63     fReco(0),
  64     fOutput(0),
  65     fTrClassNamesArr(0),
  66     fEsdEv(0),
  67     fAodEv(0),
  68     fRecPoints(0),
  69     fEsdClusters(0),
  70     fEsdCells(0),
  71     fAodClusters(0),
  72     fAodCells(0),
  73     fPtRanges(0),
  74     fNtuple(0),
  75     fSelTracks(0),
  76     fHCuts(0x0),
  77     fHVertexZ(0x0),
  78     fHVertexZ2(0x0),
  79     fHCent(0x0),
  80     fHCentQual(0x0),
  81     fHTclsBeforeCuts(0x0),
  82     fHTclsAfterCuts(0x0),
  83     fHColuRow(0x0),
  84     fHColuRowE(0x0),
  85     fHCellMult(0x0),
  86     fHCellE(0x0),
  87     fHCellH(0x0),
  88     fHCellM(0x0),
  89     fHCellM2(0x0),
  90     fHCellFreqNoCut(0x0),
  91     fHCellFrequCut100M(0x0),
  92     fHCellFrequCut300M(0x0),
  93     fHCellCheckE(0x0),
  94     fHClustEccentricity(0),
  95     fHClustEtaPhi(0x0),
  96     fHClustEnergyPt(0x0),
  97     fHClustEnergySigma(0x0),
  98     fHClustSigmaSigma(0x0),
  99     fHClustNCellEnergyRatio(0x0),
 100     fHPionEtaPhi(0x0),
 101     fHPionMggPt(0x0),
 102     fHPionMggAsym(0x0),
 103     fHPionMggDgg(0x0)
 104 {
 105   // Constructor.
 106
 107   if (!name)
 108     return;
 109   SetName(name);
 110   DefineInput(0, TChain::Class());
 111   DefineOutput(1, TList::Class());
 112   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks "
 113                "AOD:header,vertices,emcalCells,tracks";
 114 }
 115
 116 //________________________________________________________________________
 117 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
 118 {
 119   // Destructor.
 120
 121   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 122     delete fOutput; fOutput = 0;
 123   }
 124   delete fPtRanges; fPtRanges = 0;
 125   delete fGeom; fGeom = 0;
 126   delete fReco; fReco = 0;
 127   delete fTrClassNamesArr;
 128   delete fSelTracks;
 129   delete [] fHColuRow;
 130   delete [] fHColuRowE;
 131   delete [] fHCellMult;
 132   delete [] fHCellFreqNoCut;
 133   delete [] fHCellFrequCut100M;
 134   delete [] fHCellFrequCut300M;
 135 }
 136
 137 //________________________________________________________________________
 138 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 139 {
 140   // Create user objects here.
 141
 142   fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
 143   fReco = new AliEMCALRecoUtils();
 144   fTrClassNamesArr = fTrClassNames.Tokenize(" ");
 145   fOutput = new TList();
 146   fOutput->SetOwner();
 147   fSelTracks = new TObjArray;
 148
 149   if (fDoNtuple) {
 150     TFile *f = OpenFile(1);
 151     if (f) {
 152       f->SetCompressionLevel(2);
 153       fNtuple = new TNtuple(Form("nt%.0fto%.0f",fCentFrom,fCentTo),"nt",
 154                             "run:evt:l0:tcls:cent:pt:eta:phi:e:emax:n:n1:idmax:nsm:db:disp:mn:ms:ecc:sig:tkdz:tkdr:tkep:tkiso:ceiso");
 155       fNtuple->SetDirectory(f);
 156       fNtuple->SetAutoFlush(-1024*1024*1024);
 157       fNtuple->SetAutoSave(-1024*1024*1024);
 158     }
 159   }
 160
 161   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 162   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
 163   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
 164
 165   // histograms
 166   TH1::SetDefaultSumw2(kTRUE);
 167   TH2::SetDefaultSumw2(kTRUE);
 168   fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
 169   fHCuts->GetXaxis()->SetBinLabel(1,"All");
 170   fHCuts->GetXaxis()->SetBinLabel(2,"PS");
 171   fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 172   fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
 173   fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 174   fOutput->Add(fHCuts);
 175   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 176   fHVertexZ->SetXTitle("z [cm]");
 177   fOutput->Add(fHVertexZ);
 178   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 179   fHVertexZ2->SetXTitle("z [cm]");
 180   fOutput->Add(fHVertexZ2);
 181   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 182   fHCent->SetXTitle(fCentVar.Data());
 183   fOutput->Add(fHCent);
 184   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 185   fHCentQual->SetXTitle(fCentVar.Data());
 186   fOutput->Add(fHCentQual);
 187   fHTclsBeforeCuts = new TH1F("fHTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 188   fHTclsAfterCuts = new TH1F("fHTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 189   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 190     const char *name = fTrClassNamesArr->At(i)->GetName();
 191     fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
 192     fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
 193   }
 194   fOutput->Add(fHTclsBeforeCuts);
 195   fOutput->Add(fHTclsAfterCuts);
 196
 197   // histograms for cells
 198   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 199   fHColuRow   = new TH2*[nsm];
 200   fHColuRowE  = new TH2*[nsm];
 201   fHCellMult  = new TH1*[nsm];
 202   for (Int_t i = 0; i<nsm; ++i) {
 203     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",49,0,49,25,0,25);
 204     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 205     fHColuRow[i]->SetXTitle("col (i#eta)");
 206     fHColuRow[i]->SetYTitle("row (i#phi)");
 207     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",49,0,49,25,0,25);
 208     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 209     fHColuRowE[i]->SetXTitle("col (i#eta)");
 210     fHColuRowE[i]->SetYTitle("row (i#phi)");
 211     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 212     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 213     fHCellMult[i]->SetXTitle("# of cells");
 214     fOutput->Add(fHColuRow[i]);
 215     fOutput->Add(fHColuRowE[i]);
 216     fOutput->Add(fHCellMult[i]);
 217   }
 218   fHCellE = new TH1F("hCellE","",150,0.,15.);
 219   fHCellE->SetXTitle("E_{cell} [GeV]");
 220   fOutput->Add(fHCellE);
 221   fHCellH = new TH1F ("fHCellHighestE","",150,0.,15.);
 222   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 223   fOutput->Add(fHCellH);
 224   fHCellM = new TH1F ("fHCellMeanEperHitCell","",250,0.,2.5);
 225   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 226   fOutput->Add(fHCellM);
 227   fHCellM2 = new TH1F ("fHCellMeanEperAllCells","",250,0.,1);
 228   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 229   fOutput->Add(fHCellM2);
 230
 231   fHCellFreqNoCut    = new TH1*[nsm];
 232   fHCellFrequCut100M = new TH1*[nsm];
 233   fHCellFrequCut300M = new TH1*[nsm];
 234   for (Int_t i = 0; i<nsm; ++i){
 235     Double_t lbin = i*24*48-0.5;
 236     Double_t hbin = lbin+24*48;
 237     fHCellFreqNoCut[i]    = new TH1F(Form("fHCellFreqNoCut_SM%d",i),
 238                                      Form("Frequency SM%d (no cut);id;#",i),  1152, lbin, hbin);
 239     fHCellFrequCut100M[i] = new TH1F(Form("fHCellFreqCut100M_SM%d",i),
 240                                      Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
 241     fHCellFrequCut300M[i] = new TH1F(Form("fHCellFreqCut300M_SM%d",i),
 242                                      Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
 243     fOutput->Add(fHCellFreqNoCut[i]);
 244     fOutput->Add(fHCellFrequCut100M[i]);
 245     fOutput->Add(fHCellFrequCut300M[i]);
 246   }
 247   if (1) {
 248     fHCellCheckE = new TH1*[24*48*nsm];
 249     memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
 250     Int_t tcs[1] = {4102};
 251     for (UInt_t i = 0; i<sizeof(tcs)/sizeof(Int_t); ++i){
 252       Int_t c=tcs[i];
 253       if (c<24*48*nsm) {
 254         fHCellCheckE[i] = new TH1F(Form("fHCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 500, 0, 8);
 255         fOutput->Add(fHCellCheckE[i]);
 256       }
 257     }
 258   }
 259
 260   // histograms for clusters
 261   fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 262   fHClustEccentricity->SetXTitle("#epsilon_{C}");
 263   fOutput->Add(fHClustEccentricity);
 264   fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
 265   fHClustEtaPhi->SetXTitle("#eta");
 266   fHClustEtaPhi->SetYTitle("#varphi");
 267   fOutput->Add(fHClustEtaPhi);
 268   fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 269   fHClustEnergyPt->SetXTitle("E [GeV]");
 270   fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 271   fOutput->Add(fHClustEnergyPt);
 272   fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 273   fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 274   fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 275   fOutput->Add(fHClustEnergySigma);
 276   fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 277   fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 278   fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 279   fOutput->Add(fHClustSigmaSigma);
 280   fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 281   fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 282   fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 283   fOutput->Add(fHClustNCellEnergyRatio);
 284
 285   // histograms for pion candidates
 286   fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
 287   fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 288   fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 289   fOutput->Add(fHPionEtaPhi);
 290   fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 291   fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 292   fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 293   fOutput->Add(fHPionMggPt);
 294   fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 295   fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 296   fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 297   fOutput->Add(fHPionMggAsym);
 298   fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 299   fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 300   fHPionMggDgg->SetYTitle("opening angle [grad]");
 301   fOutput->Add(fHPionMggDgg);
 302   const Int_t nbins = 20;
 303   Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
 304   fPtRanges = new TAxis(nbins-1,xbins);
 305   for (Int_t i = 0; i<=nbins; ++i) {
 306     fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 307     fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 308     if (i==0)
 309       fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 310     else if (i==nbins)
 311       fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 312     else
 313       fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 314     fOutput->Add(fHPionInvMasses[i]);
 315   }
 316
 317   PostData(1, fOutput);
 318 }
 319
 320 //________________________________________________________________________
 321 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 322 {
 323   // Called for each event.
 324
 325   if (!InputEvent())
 326     return;
 327
 328   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 329   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 330   if (fEsdEv) {
 331     am->LoadBranch("AliESDRun.");
 332     am->LoadBranch("AliESDHeader.");
 333   } else {
 334     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 335     if (!fAodEv) {
 336       AliFatal("Neither ESD nor AOD event found");
 337       return;
 338     }
 339     am->LoadBranch("header");
 340   }
 341
 342   Int_t cut = 1;
 343   fHCuts->Fill(cut++);
 344
 345   TString trgclasses;
 346   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 347   if (h) {
 348     trgclasses = fEsdEv->GetFiredTriggerClasses();
 349   } else {
 350     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 351     if (h2)
 352       trgclasses = h2->GetFiredTriggerClasses();
 353   }
 354   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 355     const char *name = fTrClassNamesArr->At(i)->GetName();
 356     if (trgclasses.Contains(name))
 357       fHTclsBeforeCuts->Fill(1+i);
 358   }
 359
 360   UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 361   if (res==0)
 362     return;
 363
 364   if (fHCuts->GetBinContent(2)==0) {
 365     if (!AliGeomManager::GetGeometry()) { // get geometry
 366       AliWarning("Accessing geometry from OCDB, this is not very efficient!");
 367       AliCDBManager *cdb = AliCDBManager::Instance();
 368       if (!cdb->IsDefaultStorageSet())
 369         cdb->SetDefaultStorage("raw://");
 370       Int_t runno = InputEvent()->GetRunNumber();
 371       if (runno != cdb->GetRun())
 372         cdb->SetRun(runno);
 373       AliGeomManager::LoadGeometry();
 374     }
 375
 376     if (!fGeom->GetMatrixForSuperModule(0)) { // set misalignment matrices (stored in first event)
 377       if (fEsdEv) {
 378         for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
 379           fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 380       } else {
 381         for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
 382           fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 383       }
 384     }
 385
 386     if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 387       if (fEsdEv) {
 388         const AliESDRun *erun = fEsdEv->GetESDRun();
 389         AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 390                                                  erun->GetCurrentDip(),
 391                                                  AliMagF::kConvLHC,
 392                                                  kFALSE,
 393                                                  erun->GetBeamEnergy(),
 394                                                  erun->GetBeamType());
 395         TGeoGlobalMagField::Instance()->SetField(field);
 396       } else {
 397         Double_t pol = -1; //polarity
 398         Double_t be = -1;  //beam energy
 399         AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 400         Int_t runno = fAodEv->GetRunNumber();
 401         if (runno>=136851 && runno<138275) {
 402           pol = -1;
 403           be = 2760;
 404           btype = AliMagF::kBeamTypeAA;
 405         } else if (runno>=138275 && runno<=139517) {
 406           pol = +1;
 407           be = 2760;
 408           btype = AliMagF::kBeamTypeAA;
 409         } else {
 410           AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 411         }
 412         TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 413       }
 414     }
 415   }
 416
 417   fHCuts->Fill(cut++);
 418
 419   const AliCentrality *centP = InputEvent()->GetCentrality();
 420   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 421   fHCent->Fill(cent);
 422   if (cent<fCentFrom||cent>fCentTo)
 423     return;
 424
 425   fHCuts->Fill(cut++);
 426
 427   if (fUseQualFlag) {
 428     if (centP->GetQuality()>0)
 429       return;
 430   }
 431
 432   fHCentQual->Fill(cent);
 433   fHCuts->Fill(cut++);
 434
 435   if (fEsdEv) {
 436     am->LoadBranch("PrimaryVertex.");
 437     am->LoadBranch("SPDVertex.");
 438     am->LoadBranch("TPCVertex.");
 439   } else {
 440     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 441     am->LoadBranch("vertices");
 442     if (!fAodEv) return;
 443   }
 444
 445   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 446   if (!vertex)
 447     return;
 448
 449   fHVertexZ->Fill(vertex->GetZ());
 450
 451   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 452     return;
 453
 454   fHCuts->Fill(cut++);
 455   fHVertexZ2->Fill(vertex->GetZ());
 456
 457   // count number of accepted events
 458   ++fNEvs;
 459
 460   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 461     const char *name = fTrClassNamesArr->At(i)->GetName();
 462     if (trgclasses.Contains(name))
 463       fHTclsAfterCuts->Fill(1+i);
 464   }
 465
 466   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 467   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
 468   fEsdCells    = 0; // will be set if ESD input used
 469   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
 470                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 471   fAodCells    = 0; // will be set if AOD input used
 472
 473   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 474   Bool_t clusattached = 0;
 475   Bool_t recalibrated = 0;
 476   if (1 && !fClusName.IsNull()) {
 477     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 478     TObjArray *ts = am->GetTasks();
 479     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 480     if (cltask && cltask->GetClusters()) {
 481       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 482       clusattached = cltask->GetAttachClusters();
 483       if (cltask->GetCalibData()!=0)
 484         recalibrated = kTRUE;
 485     }
 486   }
 487   if (1 && AODEvent() && !fClusName.IsNull()) {
 488     TList *l = AODEvent()->GetList();
 489     TClonesArray *clus = 0;
 490     if (l) {
 491       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 492       fAodClusters = clus;
 493     }
 494   }
 495
 496   if (fEsdEv) { // ESD input mode
 497     if (1 && (!fRecPoints||clusattached)) {
 498       if (!clusattached)
 499         am->LoadBranch("CaloClusters");
 500       TList *l = fEsdEv->GetList();
 501       TClonesArray *clus = 0;
 502       if (l) {
 503         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 504         fEsdClusters = clus;
 505       }
 506     }
 507     if (1) {
 508       if (!recalibrated)
 509         am->LoadBranch("EMCALCells.");
 510       fEsdCells = fEsdEv->GetEMCALCells();
 511     }
 512   } else if (fAodEv) { // AOD input mode
 513     if (1 && (!fAodClusters || clusattached)) {
 514       if (!clusattached)
 515         am->LoadBranch("caloClusters");
 516       TList *l = fAodEv->GetList();
 517       TClonesArray *clus = 0;
 518       if (l) {
 519         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 520         fAodClusters = clus;
 521       }
 522     }
 523     if (1) {
 524       if (!recalibrated)
 525         am->LoadBranch("emcalCells");
 526       fAodCells = fAodEv->GetEMCALCells();
 527     }
 528   } else {
 529     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 530   }
 531
 532   if (1) {
 533     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 534     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 535     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 536     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 537     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 538   }
 539
 540   if (fDoAfterburner)
 541     ClusterAfterburner();
 542
 543   CalcTracks();
 544   CalcClusterProps();
 545
 546   FillCellHists();
 547   FillClusHists();
 548   FillPionHists();
 549   FillOtherHists();
 550
 551   PostData(1, fOutput);
 552 }
 553
 554 //________________________________________________________________________
 555 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 556 {
 557   // Terminate called at the end of analysis.
 558
 559   if (fNtuple) {
 560     TFile *f = OpenFile(1);
 561     if (f)
 562       fNtuple->Write();
 563   }
 564
 565   AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
 566 }
 567
 568 //________________________________________________________________________
 569 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
 570 {
 571   // Calculate track properties.
 572
 573   fSelTracks->Clear();
 574
 575   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 576   TClonesArray *tracks = 0;
 577   if (fEsdEv) {
 578     am->LoadBranch("Tracks");
 579     TList *l = fEsdEv->GetList();
 580     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
 581   } else {
 582     am->LoadBranch("tracks");
 583     TList *l = fAodEv->GetList();
 584     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
 585   }
 586
 587   if (!tracks)
 588     return;
 589
 590   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 591   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
 592   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
 593
 594   if (fEsdEv) {
 595     fSelTracks->SetOwner(kTRUE);
 596     am->LoadBranch("PrimaryVertex.");
 597     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
 598     am->LoadBranch("SPDVertex.");
 599     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
 600     am->LoadBranch("Tracks");
 601     const Int_t Ntracks = tracks->GetEntries();
 602     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 603       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
 604       if (!track) {
 605         AliWarning(Form("Could not receive track %d\n", iTracks));
 606         continue;
 607       }
 608       if (fTrCuts && !fTrCuts->IsSelected(track))
 609         continue;
 610       Double_t eta = track->Eta();
 611       if (eta<-1||eta>1)
 612         continue;
 613       if (track->Phi()<phimin||track->Phi()>phimax)
 614         continue;
 615       if(track->GetTPCNcls()<fMinNClustPerTrack)
 616         continue;
 617       AliESDtrack copyt(*track);
 618       Double_t bfield[3];
 619       copyt.GetBxByBz(bfield);
 620       AliExternalTrackParam tParam;
 621       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
 622       if (!relate)
 623         continue;
 624       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
 625       Double_t p[3]      = { 0. };
 626       copyt.GetPxPyPz(p);
 627       Double_t pos[3]    = { 0. };
 628       copyt.GetXYZ(pos);
 629       Double_t covTr[21] = { 0. };
 630       copyt.GetCovarianceXYZPxPyPz(covTr);
 631       Double_t pid[10]   = { 0. };
 632       copyt.GetESDpid(pid);
 633       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
 634                                             copyt.GetLabel(),
 635                                             p,
 636                                             kTRUE,
 637                                             pos,
 638                                             kFALSE,
 639                                             covTr,
 640                                             (Short_t)copyt.GetSign(),
 641                                             copyt.GetITSClusterMap(),
 642                                             pid,
 643                                             0,/*fPrimaryVertex,*/
 644                                             kTRUE, // check if this is right
 645                                             vtx->UsesTrack(copyt.GetID()));
 646       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
 647       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
 648       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
 649       if(ndf>0)
 650         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
 651       else
 652         aTrack->SetChi2perNDF(-1);
 653       aTrack->SetFlags(copyt.GetStatus());
 654       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
 655       fSelTracks->Add(aTrack);
 656     }
 657   } else {
 658     Int_t ntracks = tracks->GetEntries();
 659     for (Int_t i=0; i<ntracks; ++i) {
 660       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
 661       if (!track)
 662         continue;
 663       Double_t eta = track->Eta();
 664       if (eta<-1||eta>1)
 665         continue;
 666       if (track->Phi()<phimin||track->Phi()>phimax)
 667         continue;
 668       if(track->GetTPCNcls()<fMinNClustPerTrack)
 669         continue;
 670
 671       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
 672         AliExternalTrackParam tParam(track);
 673         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
 674           Double_t trkPos[3];
 675           tParam.GetXYZ(trkPos);
 676           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
 677         }
 678       }
 679       fSelTracks->Add(track);
 680     }
 681   }
 682 }
 683
 684 //________________________________________________________________________
 685 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 686 {
 687   // Calculate cluster properties
 688
 689   TObjArray *clusters = fEsdClusters;
 690   if (!clusters)
 691     clusters = fAodClusters;
 692   if (!clusters)
 693     return;
 694
 695   Int_t nclus = clusters->GetEntries();
 696   Int_t ntrks = fSelTracks->GetEntries();
 697
 698   for(Int_t i = 0; i<nclus; ++i) {
 699     fClusProps[i].Reset();
 700
 701     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 702     if (!clus)
 703       continue;
 704     if (!clus->IsEMCAL())
 705       continue;
 706     if (clus->E()<fMinE)
 707       continue;
 708
 709     Float_t  clsPos[3] = {0};
 710     clus->GetPosition(clsPos);
 711     TVector3 clsVec(clsPos);
 712     Double_t vertex[3] = {0};
 713     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 714     TLorentzVector clusterVec;
 715     clus->GetMomentum(clusterVec,vertex);
 716     Double_t clsEta = clusterVec.Eta();
 717
 718     Double_t mind2 = 1e10;
 719     for(Int_t j = 0; j<ntrks; ++j) {
 720       AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
 721       if (!track)
 722         continue;
 723       Double_t pt  = track->Pt();
 724       if (pt<fMinPtPerTrack)
 725         continue;
 726       if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
 727         continue;
 728
 729       AliExternalTrackParam tParam(track);
 730       Float_t tmpR=-1, tmpZ=-1;
 731 #if 1
 732       if (1) {
 733         TVector3 vec(clsPos);
 734         Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
 735         vec.RotateZ(-alpha);  //Rotate the cluster to the local extrapolation coordinate system
 736         tParam.Rotate(alpha); //Rotate the track to the same local extrapolation system
 737         if (!AliTrackerBase::PropagateTrackToBxByBz(&tParam, vec.X(), 0.139, 1, kFALSE))
 738           continue;
 739         Double_t trkPos[3];
 740         tParam.GetXYZ(trkPos); //Get the extrapolated global position
 741         tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
 742         tmpZ = clsPos[2]-trkPos[2];
 743       }
 744 #else
 745       if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
 746         continue;
 747 #endif
 748       Double_t d2 = tmpR;
 749       if (mind2>d2) {
 750         mind2=d2;
 751         fClusProps[i].fTrIndex = j;
 752         fClusProps[i].fTrDz    = tmpZ;
 753         fClusProps[i].fTrDr    = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
 754         fClusProps[i].fTrDist  = d2;
 755         fClusProps[i].fTrEp    = clus->E()/track->P();
 756       }
 757     }
 758
 759     if (0 && (fClusProps[i].fTrIndex>=0)) {
 760       cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
 761     }
 762     fClusProps[i].fTrIso      = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
 763     fClusProps[i].fTrLowPtIso = 0;
 764     fClusProps[i].fCellIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
 765   }
 766 }
 767
 768 //________________________________________________________________________
 769 void  AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
 770 {
 771   // Run custer reconstruction afterburner.
 772
 773   AliVCaloCells *cells = fEsdCells;
 774   if (!cells)
 775     cells = fAodCells;
 776
 777   if (!cells)
 778     return;
 779
 780   Int_t ncells = cells->GetNumberOfCells();
 781   if (ncells<=0)
 782     return;
 783
 784   Double_t cellMeanE = 0, cellSigE = 0;
 785   for (Int_t i = 0; i<ncells; ++i) {
 786     Double_t cellE = cells->GetAmplitude(i);
 787     cellMeanE += cellE;
 788     cellSigE += cellE*cellE;
 789   }
 790   cellMeanE /= ncells;
 791   cellSigE /= ncells;
 792   cellSigE -= cellMeanE*cellMeanE;
 793   if (cellSigE<0)
 794     cellSigE = 0;
 795   cellSigE = TMath::Sqrt(cellSigE / ncells);
 796
 797   Double_t subE = cellMeanE - 7*cellSigE;
 798   if (subE<0)
 799     return;
 800
 801   for (Int_t i = 0; i<ncells; ++i) {
 802     Short_t id=-1;
 803     Double_t amp=0,time=0;
 804     if (!cells->GetCell(i, id, amp, time))
 805       continue;
 806     amp -= cellMeanE;
 807     if (amp<0.001)
 808       amp = 0;
 809     cells->SetCell(i, id, amp, time);
 810   }
 811
 812   TObjArray *clusters = fEsdClusters;
 813   if (!clusters)
 814     clusters = fAodClusters;
 815   if (!clusters)
 816     return;
 817
 818   Int_t nclus = clusters->GetEntries();
 819   for (Int_t i = 0; i<nclus; ++i) {
 820     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 821     if (!clus->IsEMCAL())
 822       continue;
 823     Int_t nc = clus->GetNCells();
 824     Double_t clusE = 0;
 825     UShort_t ids[100] = {0};
 826     Double_t fra[100] = {0};
 827     for (Int_t j = 0; j<nc; ++j) {
 828       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
 829       Double_t cen = cells->GetCellAmplitude(id);
 830       clusE += cen;
 831       if (cen>0) {
 832         ids[nc] = id;
 833         ++nc;
 834       }
 835     }
 836     if (clusE<=0) {
 837       clusters->RemoveAt(i);
 838       continue;
 839     }
 840
 841     for (Int_t j = 0; j<nc; ++j) {
 842       Short_t id = ids[j];
 843       Double_t cen = cells->GetCellAmplitude(id);
 844       fra[j] = cen/clusE;
 845     }
 846     clus->SetE(clusE);
 847     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
 848     if (aodclus) {
 849       aodclus->Clear("");
 850       aodclus->SetNCells(nc);
 851       aodclus->SetCellsAmplitudeFraction(fra);
 852       aodclus->SetCellsAbsId(ids);
 853     }
 854   }
 855   clusters->Compress();
 856 }
 857
 858 //________________________________________________________________________
 859 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
 860 {
 861   // Fill histograms related to cell properties.
 862
 863   AliVCaloCells *cells = fEsdCells;
 864   if (!cells)
 865     cells = fAodCells;
 866
 867   if (!cells)
 868     return;
 869
 870   Int_t cellModCount[12] = {0};
 871   Double_t cellMaxE = 0;
 872   Double_t cellMeanE = 0;
 873   Int_t ncells = cells->GetNumberOfCells();
 874   if (ncells==0)
 875     return;
 876
 877   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 878
 879   for (Int_t i = 0; i<ncells; ++i) {
 880     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
 881     Double_t cellE = cells->GetAmplitude(i);
 882     fHCellE->Fill(cellE);
 883     if (cellE>cellMaxE)
 884       cellMaxE = cellE;
 885     cellMeanE += cellE;
 886
 887     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
 888     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
 889     if (!ret) {
 890       AliError(Form("Could not get cell index for %d", absID));
 891       continue;
 892     }
 893     ++cellModCount[iSM];
 894     Int_t iPhi=-1, iEta=-1;
 895     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
 896     fHColuRow[iSM]->Fill(iEta,iPhi,1);
 897     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
 898     fHCellFreqNoCut[iSM]->Fill(absID);
 899     if (cellE > 0.1) fHCellFrequCut100M[iSM]->Fill(absID);
 900     if (cellE > 0.3) fHCellFrequCut300M[iSM]->Fill(absID);
 901     if (fHCellCheckE && fHCellCheckE[absID])
 902       fHCellCheckE[absID]->Fill(cellE);
 903   }
 904   fHCellH->Fill(cellMaxE);
 905   cellMeanE /= ncells;
 906   fHCellM->Fill(cellMeanE);
 907   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
 908   for (Int_t i=0; i<nsm; ++i)
 909     fHCellMult[i]->Fill(cellModCount[i]);
 910 }
 911
 912 //________________________________________________________________________
 913 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
 914 {
 915   // Fill histograms related to cluster properties.
 916
 917   TObjArray *clusters = fEsdClusters;
 918   if (!clusters)
 919     clusters = fAodClusters;
 920   if (!clusters)
 921     return;
 922
 923   Double_t vertex[3] = {0};
 924   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 925
 926   Int_t nclus = clusters->GetEntries();
 927   for(Int_t i = 0; i<nclus; ++i) {
 928     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 929     if (!clus)
 930       continue;
 931     if (!clus->IsEMCAL())
 932       continue;
 933     TLorentzVector clusterVec;
 934     clus->GetMomentum(clusterVec,vertex);
 935     Double_t maxAxis = 0, minAxis = 0;
 936     GetSigma(clus,maxAxis,minAxis);
 937     Double_t clusterEcc = 0;
 938     if (maxAxis > 0)
 939       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 940     clus->SetChi2(clusterEcc); // store ecc in chi2
 941     fHClustEccentricity->Fill(clusterEcc);
 942     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
 943     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
 944     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
 945     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
 946     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
 947     if (fNtuple) {
 948       if (clus->E()<fMinE)
 949         continue;
 950       Float_t vals[25];
 951       TString trgclasses;
 952       vals[0]  = InputEvent()->GetRunNumber();
 953       vals[1]  = (((UInt_t)InputEvent()->GetOrbitNumber()  << 12) | (UInt_t)InputEvent()->GetBunchCrossNumber());
 954       if (vals[1]<=0)
 955         vals[1] = fNEvs;
 956       AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 957       if (h) {
 958         vals[2]  = h->GetL0TriggerInputs();
 959         trgclasses = fEsdEv->GetFiredTriggerClasses();
 960       }
 961       else {
 962         AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 963         if (h2) {
 964           vals[2]  = h2->GetL0TriggerInputs();
 965           trgclasses = h2->GetFiredTriggerClasses();
 966         } else
 967           vals[2] = 0;
 968       }
 969       vals[3]  = 0;
 970
 971       for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
 972         const char *name = fTrClassNamesArr->At(j)->GetName();
 973         if (trgclasses.Contains(name))
 974           vals[3] += TMath::Power(2,j);
 975       }
 976       vals[4]  = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar);
 977       vals[5]  = clusterVec.Pt();
 978       vals[6]  = clusterVec.Eta();
 979       vals[7]  = clusterVec.Phi();
 980       vals[8]  = clusterVec.E();
 981       vals[9]  = GetMaxCellEnergy(clus);
 982       vals[10] = clus->GetNCells();
 983       vals[11] = GetNCells(clus,0.100);
 984       vals[12] = clus->GetCellAbsId(0);
 985       vals[13] = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0));
 986       vals[14] = clus->GetDistanceToBadChannel();
 987       vals[15] = clus->GetDispersion();
 988       vals[16] = clus->GetM20();
 989       vals[17] = clus->GetM02();
 990       vals[18] = clusterEcc;
 991       vals[19] = maxAxis;
 992       vals[20] = fClusProps[i].fTrDz;
 993       vals[21] = fClusProps[i].fTrDr;
 994       vals[22] = fClusProps[i].fTrEp;
 995       vals[23] = fClusProps[i].fTrIso;
 996       vals[24] = fClusProps[i].fCellIso;
 997       fNtuple->Fill(vals);
 998     }
 999   }
1000 }
1001
1002 //________________________________________________________________________
1003 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1004 {
1005   // Fill histograms related to pions.
1006
1007   Double_t vertex[3] = {0};
1008   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1009
1010   TObjArray *clusters = fEsdClusters;
1011   if (!clusters)
1012     clusters = fAodClusters;
1013
1014   if (clusters) {
1015     TLorentzVector clusterVec1;
1016     TLorentzVector clusterVec2;
1017     TLorentzVector pionVec;
1018
1019     Int_t nclus = clusters->GetEntries();
1020     for (Int_t i = 0; i<nclus; ++i) {
1021       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1022       if (!clus1)
1023         continue;
1024       if (!clus1->IsEMCAL())
1025         continue;
1026       if (clus1->E()<fMinE)
1027         continue;
1028       if (clus1->GetNCells()<fNminCells)
1029         continue;
1030       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1031         continue;
1032       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1033         continue;
1034       clus1->GetMomentum(clusterVec1,vertex);
1035       for (Int_t j = i+1; j<nclus; ++j) {
1036         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1037         if (!clus2)
1038           continue;
1039         if (!clus2->IsEMCAL())
1040           continue;
1041         if (clus2->E()<fMinE)
1042           continue;
1043         if (clus2->GetNCells()<fNminCells)
1044           continue;
1045         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1046           continue;
1047         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1048           continue;
1049         clus2->GetMomentum(clusterVec2,vertex);
1050         pionVec = clusterVec1 + clusterVec2;
1051         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1052         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1053         if (pionZgg < fAsymMax) {
1054           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1055           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1056           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1057           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1058           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1059           fHPionInvMasses[bin]->Fill(pionVec.M());
1060         }
1061       }
1062     }
1063   }
1064 }
1065
1066 //________________________________________________________________________
1067 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1068 {
1069   // Fill histograms related to cell properties.
1070 }
1071
1072 //________________________________________________________________________
1073 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1074 {
1075   // Compute isolation based on cell content.
1076
1077   AliVCaloCells *cells = fEsdCells;
1078   if (!cells)
1079     cells = fAodCells;
1080   if (!cells)
1081     return 0;
1082
1083   Double_t cellIsolation = 0;
1084   Double_t rad2 = radius*radius;
1085   Int_t ncells = cells->GetNumberOfCells();
1086   for (Int_t i = 0; i<ncells; ++i) {
1087     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1088     Double_t cellE = cells->GetAmplitude(i);
1089     Float_t eta, phi;
1090     fGeom->EtaPhiFromIndex(absID,eta,phi);
1091     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1092     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1093     if(dist>rad2)
1094       continue;
1095     cellIsolation += cellE;
1096   }
1097   return cellIsolation;
1098 }
1099
1100 //________________________________________________________________________
1101 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster) const
1102 {
1103   // Get maximum energy of attached cell.
1104
1105   Double_t maxe = 0;
1106   Int_t ncells = cluster->GetNCells();
1107   if (fEsdCells) {
1108     for (Int_t i=0; i<ncells; i++) {
1109       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1110       if (e>maxe) {
1111         maxe = e;
1112       }
1113     }
1114   } else {
1115     for (Int_t i=0; i<ncells; i++) {
1116       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1117       if (e>maxe)
1118         maxe = e;
1119     }
1120   }
1121   return maxe;
1122 }
1123
1124 //________________________________________________________________________
1125 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
1126 {
1127   // Calculate the (E) weighted variance along the longer (eigen) axis.
1128
1129   sigmaMax = 0;          // cluster variance along its longer axis
1130   sigmaMin = 0;          // cluster variance along its shorter axis
1131   Double_t Ec  = c->E(); // cluster energy
1132   if(Ec<=0)
1133     return;
1134   Double_t Xc  = 0 ;     // cluster first moment along X
1135   Double_t Yc  = 0 ;     // cluster first moment along Y
1136   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
1137   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
1138   Double_t Syy = 0 ;     // cluster covariance^2
1139
1140   AliVCaloCells *cells = fEsdCells;
1141   if (!cells)
1142     cells = fAodCells;
1143
1144   if (!cells)
1145     return;
1146
1147   Int_t ncells = c->GetNCells();
1148   if (ncells==1)
1149     return;
1150
1151   TVector3 pos;
1152   for(Int_t j=0; j<ncells; ++j) {
1153     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1154     fGeom->GetGlobal(id,pos);
1155     Double_t cellen = cells->GetCellAmplitude(id);
1156     Xc  += cellen*pos.X();
1157     Yc  += cellen*pos.Y();
1158     Sxx += cellen*pos.X()*pos.X();
1159     Syy += cellen*pos.Y()*pos.Y();
1160     Sxy += cellen*pos.X()*pos.Y();
1161   }
1162   Xc  /= Ec;
1163   Yc  /= Ec;
1164   Sxx /= Ec;
1165   Syy /= Ec;
1166   Sxy /= Ec;
1167   Sxx -= Xc*Xc;
1168   Syy -= Yc*Yc;
1169   Sxy -= Xc*Yc;
1170   Sxx = TMath::Abs(Sxx);
1171   Syy = TMath::Abs(Syy);
1172   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1173   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1174   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1175   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
1176 }
1177
1178 //________________________________________________________________________
1179 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
1180 {
1181   // Calculate number of attached cells above emin.
1182
1183   AliVCaloCells *cells = fEsdCells;
1184   if (!cells)
1185     cells = fAodCells;
1186   if (!cells)
1187     return 0;
1188
1189   Int_t n = 0;
1190   Int_t ncells = c->GetNCells();
1191   for(Int_t j=0; j<ncells; ++j) {
1192     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1193     Double_t cellen = cells->GetCellAmplitude(id);
1194     if (cellen>=emin)
1195       ++n;
1196   }
1197   return n;
1198 }
1199
1200 //________________________________________________________________________
1201 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1202 {
1203   // Compute isolation based on tracks.
1204
1205   Double_t trkIsolation = 0;
1206   Double_t rad2 = radius*radius;
1207   Int_t ntrks = fSelTracks->GetEntries();
1208   for(Int_t j = 0; j<ntrks; ++j) {
1209     AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1210     if (!track)
1211       continue;
1212     Float_t eta = track->Eta();
1213     Float_t phi = track->Phi();
1214     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1215     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1216     if(dist>rad2)
1217       continue;
1218     trkIsolation += track->Pt();
1219   }
1220   return trkIsolation;
1221 }