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 "AliESDCaloTrigger.h"
  27 #include "AliESDEvent.h"
  28 #include "AliESDVertex.h"
  29 #include "AliESDtrack.h"
  30 #include "AliESDtrackCuts.h"
  31 #include "AliEventplane.h"
  32 #include "AliGeomManager.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliLog.h"
  35 #include "AliMagF.h"
  36 #include "AliTrackerBase.h"
  37
  38 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  42   : AliAnalysisTaskSE(name),
  43     fCentVar("V0M"),
  44     fCentFrom(0),
  45     fCentTo(100),
  46     fVtxZMin(-10),
  47     fVtxZMax(+10),
  48     fUseQualFlag(1),
  49     fClusName(),
  50     fDoNtuple(0),
  51     fDoAfterburner(0),
  52     fAsymMax(1),
  53     fNminCells(1),
  54     fMinE(0.100),
  55     fMinErat(0),
  56     fMinEcc(0),
  57     fGeoName("EMCAL_FIRSTYEARV1"),
  58     fMinNClusPerTr(50),
  59     fIsoDist(0.2),
  60     fTrClassNames(""),
  61     fTrCuts(0),
  62     fPrimTrCuts(0),
  63     fDoTrMatGeom(0),
  64     fTrainMode(0),
  65     fMarkCells(),
  66     fMinL0Time(-1),
  67     fMaxL0Time(1024),
  68     fIsGeoMatsSet(0),
  69     fNEvs(0),
  70     fGeom(0),
  71     fReco(0),
  72     fOutput(0),
  73     fTrClassNamesArr(0),
  74     fEsdEv(0),
  75     fAodEv(0),
  76     fRecPoints(0),
  77     fEsdClusters(0),
  78     fEsdCells(0),
  79     fAodClusters(0),
  80     fAodCells(0),
  81     fPtRanges(0),
  82     fSelTracks(0),
  83     fSelPrimTracks(0),
  84     fNAmpInTrigger(0),
  85     fAmpInTrigger(0),
  86     fNtuple(0),
  87     fHeader(0),
  88     fPrimVert(0),
  89     fSpdVert(0),
  90     fTpcVert(0),
  91     fClusters(0),
  92     fTriggers(0),
  93     fHCuts(0x0),
  94     fHVertexZ(0x0),
  95     fHVertexZ2(0x0),
  96     fHCent(0x0),
  97     fHCentQual(0x0),
  98     fHTclsBeforeCuts(0x0),
  99     fHTclsAfterCuts(0x0),
 100     fHColuRow(0x0),
 101     fHColuRowE(0x0),
 102     fHCellMult(0x0),
 103     fHCellE(0x0),
 104     fHCellH(0x0),
 105     fHCellM(0x0),
 106     fHCellM2(0x0),
 107     fHCellFreqNoCut(0x0),
 108     fHCellFreqCut100M(0x0),
 109     fHCellFreqCut300M(0x0),
 110     fHCellFreqE(0x0),
 111     fHCellCheckE(0x0),
 112     fHClustEccentricity(0),
 113     fHClustEtaPhi(0x0),
 114     fHClustEnergyPt(0x0),
 115     fHClustEnergySigma(0x0),
 116     fHClustSigmaSigma(0x0),
 117     fHClustNCellEnergyRatio(0x0),
 118     fHMatchDr(0x0),
 119     fHMatchDz(0x0),
 120     fHMatchEp(0x0),
 121     fHPionEtaPhi(0x0),
 122     fHPionMggPt(0x0),
 123     fHPionMggAsym(0x0),
 124     fHPionMggDgg(0x0)
 125 {
 126   // Constructor.
 127
 128   if (!name)
 129     return;
 130   SetName(name);
 131   DefineInput(0, TChain::Class());
 132   DefineOutput(1, TList::Class());
 133   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks  EMCALTrigger."
 134                "AOD:header,vertices,emcalCells,tracks";
 135 }
 136
 137 //________________________________________________________________________
 138 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
 139 {
 140   // Destructor.
 141
 142   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 143     delete fOutput; fOutput = 0;
 144   }
 145   delete fPtRanges; fPtRanges = 0;
 146   delete fGeom; fGeom = 0;
 147   delete fReco; fReco = 0;
 148   delete fTrClassNamesArr;
 149   delete fSelTracks;
 150   delete fSelPrimTracks;
 151   delete [] fAmpInTrigger;
 152   delete [] fHColuRow;
 153   delete [] fHColuRowE;
 154   delete [] fHCellMult;
 155   delete [] fHCellFreqNoCut;
 156   delete [] fHCellFreqCut100M;
 157   delete [] fHCellFreqCut300M;
 158 }
 159
 160 //________________________________________________________________________
 161 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 162 {
 163   // Create user objects here.
 164
 165   fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
 166   fReco = new AliEMCALRecoUtils();
 167   fTrClassNamesArr = fTrClassNames.Tokenize(" ");
 168   fOutput = new TList();
 169   fOutput->SetOwner();
 170   fSelTracks = new TObjArray;
 171   fSelPrimTracks = new TObjArray;
 172
 173   if (fDoNtuple) {
 174     TFile *f = OpenFile(1);
 175     if (f) {
 176       f->SetCompressionLevel(2);
 177       fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
 178       fNtuple->SetDirectory(f);
 179       if (fTrainMode) {
 180         fNtuple->SetAutoFlush(-2*1024*1024);
 181         fNtuple->SetAutoSave(0);
 182       } else {
 183         fNtuple->SetAutoFlush(-32*1024*1024);
 184         fNtuple->SetAutoSave(0);
 185       }
 186
 187       fHeader = new AliStaHeader;
 188       fNtuple->Branch("header", &fHeader, 16*1024, 99);
 189       fPrimVert = new AliStaVertex;
 190       fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
 191       fSpdVert = new AliStaVertex;
 192       fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
 193       fTpcVert = new AliStaVertex;
 194       fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
 195       if (TClass::GetClass("AliStaCluster"))
 196         TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
 197       fClusters = new TClonesArray("AliStaCluster");
 198       fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
 199       if (TClass::GetClass("AliStaTrigger"))
 200         TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
 201       fTriggers = new TClonesArray("AliStaTrigger");
 202       fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
 203     }
 204   }
 205
 206   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 207   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
 208   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
 209
 210   // histograms
 211   TH1::SetDefaultSumw2(kTRUE);
 212   TH2::SetDefaultSumw2(kTRUE);
 213   fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
 214   fHCuts->GetXaxis()->SetBinLabel(1,"All");
 215   fHCuts->GetXaxis()->SetBinLabel(2,"PS");
 216   fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 217   fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
 218   fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 219   fOutput->Add(fHCuts);
 220   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 221   fHVertexZ->SetXTitle("z [cm]");
 222   fOutput->Add(fHVertexZ);
 223   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 224   fHVertexZ2->SetXTitle("z [cm]");
 225   fOutput->Add(fHVertexZ2);
 226   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 227   fHCent->SetXTitle(fCentVar.Data());
 228   fOutput->Add(fHCent);
 229   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 230   fHCentQual->SetXTitle(fCentVar.Data());
 231   fOutput->Add(fHCentQual);
 232   fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 233   fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 234   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 235     const char *name = fTrClassNamesArr->At(i)->GetName();
 236     fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
 237     fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
 238   }
 239   fOutput->Add(fHTclsBeforeCuts);
 240   fOutput->Add(fHTclsAfterCuts);
 241
 242   // histograms for cells
 243   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 244   fHColuRow   = new TH2*[nsm];
 245   fHColuRowE  = new TH2*[nsm];
 246   fHCellMult  = new TH1*[nsm];
 247   for (Int_t i = 0; i<nsm; ++i) {
 248     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
 249     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 250     fHColuRow[i]->SetXTitle("col (i#eta)");
 251     fHColuRow[i]->SetYTitle("row (i#phi)");
 252     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
 253     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 254     fHColuRowE[i]->SetXTitle("col (i#eta)");
 255     fHColuRowE[i]->SetYTitle("row (i#phi)");
 256     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 257     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 258     fHCellMult[i]->SetXTitle("# of cells");
 259     fOutput->Add(fHColuRow[i]);
 260     fOutput->Add(fHColuRowE[i]);
 261     fOutput->Add(fHCellMult[i]);
 262   }
 263   fHCellE = new TH1F("hCellE","",250,0.,25.);
 264   fHCellE->SetXTitle("E_{cell} [GeV]");
 265   fOutput->Add(fHCellE);
 266   fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
 267   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 268   fOutput->Add(fHCellH);
 269   fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
 270   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 271   fOutput->Add(fHCellM);
 272   fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
 273   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 274   fOutput->Add(fHCellM2);
 275
 276   fHCellFreqNoCut   = new TH1*[nsm];
 277   fHCellFreqCut100M = new TH1*[nsm];
 278   fHCellFreqCut300M = new TH1*[nsm];
 279   fHCellFreqE       = new TH1*[nsm];
 280   for (Int_t i = 0; i<nsm; ++i){
 281     Double_t lbin = i*24*48-0.5;
 282     Double_t hbin = lbin+24*48;
 283     fHCellFreqNoCut[i]   = new TH1F(Form("hCellFreqNoCut_SM%d",i),
 284                                     Form("Frequency SM%d (no cut);id;#",i),  1152, lbin, hbin);
 285     fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
 286                                     Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
 287     fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
 288                                     Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
 289     fHCellFreqE[i]       = new TH1F(Form("hCellFreqE_SM%d",i),
 290                                     Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
 291     fOutput->Add(fHCellFreqNoCut[i]);
 292     fOutput->Add(fHCellFreqCut100M[i]);
 293     fOutput->Add(fHCellFreqCut300M[i]);
 294     fOutput->Add(fHCellFreqE[i]);
 295   }
 296   if (!fMarkCells.IsNull()) {
 297     fHCellCheckE = new TH1*[24*48*nsm];
 298     memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
 299     TObjArray *cells = fMarkCells.Tokenize(" ");
 300     Int_t n = cells->GetEntries();
 301     Int_t *tcs = new Int_t[n];
 302     for (Int_t i=0;i<n;++i) {
 303       TString name(cells->At(i)->GetName());
 304       tcs[i]=name.Atoi();
 305     }
 306     for (Int_t i = 0; i<n; ++i) {
 307       Int_t c=tcs[i];
 308       if (c<24*48*nsm) {
 309         fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
 310         fOutput->Add(fHCellCheckE[i]);
 311       }
 312     }
 313     delete cells;
 314     delete [] tcs;
 315   }
 316
 317   // histograms for clusters
 318   if (!fTrainMode) {
 319     fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 320     fHClustEccentricity->SetXTitle("#epsilon_{C}");
 321     fOutput->Add(fHClustEccentricity);
 322     fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
 323     fHClustEtaPhi->SetXTitle("#eta");
 324     fHClustEtaPhi->SetYTitle("#varphi");
 325     fOutput->Add(fHClustEtaPhi);
 326     fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 327     fHClustEnergyPt->SetXTitle("E [GeV]");
 328     fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 329     fOutput->Add(fHClustEnergyPt);
 330     fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 331     fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 332     fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 333     fOutput->Add(fHClustEnergySigma);
 334     fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 335     fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 336     fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 337     fOutput->Add(fHClustSigmaSigma);
 338     fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 339     fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 340     fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 341     fOutput->Add(fHClustNCellEnergyRatio);
 342   }
 343
 344   // histograms for track matching
 345   fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
 346   fOutput->Add(fHMatchDr);
 347   fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
 348   fOutput->Add(fHMatchDz);
 349   fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
 350   fOutput->Add(fHMatchEp);
 351
 352   // histograms for pion candidates
 353   if (!fTrainMode) {
 354     fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
 355     fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 356     fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 357     fOutput->Add(fHPionEtaPhi);
 358     fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 359     fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 360     fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 361     fOutput->Add(fHPionMggPt);
 362     fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 363     fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 364     fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 365     fOutput->Add(fHPionMggAsym);
 366     fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 367     fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 368     fHPionMggDgg->SetYTitle("opening angle [grad]");
 369     fOutput->Add(fHPionMggDgg);
 370     const Int_t nbins = 20;
 371     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};
 372     fPtRanges = new TAxis(nbins-1,xbins);
 373     for (Int_t i = 0; i<=nbins; ++i) {
 374       fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 375       fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 376       if (i==0)
 377         fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 378       else if (i==nbins)
 379         fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 380       else
 381         fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 382       fOutput->Add(fHPionInvMasses[i]);
 383     }
 384   }
 385
 386   PostData(1, fOutput);
 387 }
 388
 389 //________________________________________________________________________
 390 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 391 {
 392   // Called for each event.
 393
 394   if (!InputEvent())
 395     return;
 396
 397   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 398   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 399   if (fEsdEv) {
 400     am->LoadBranch("AliESDRun.");
 401     am->LoadBranch("AliESDHeader.");
 402   } else {
 403     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 404     if (!fAodEv) {
 405       AliFatal("Neither ESD nor AOD event found");
 406       return;
 407     }
 408     am->LoadBranch("header");
 409   }
 410
 411   if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
 412     AliWarning("Accessing geometry from OCDB, this is not very efficient!");
 413     AliCDBManager *cdb = AliCDBManager::Instance();
 414     if (!cdb->IsDefaultStorageSet())
 415       cdb->SetDefaultStorage("raw://");
 416     Int_t runno = InputEvent()->GetRunNumber();
 417     if (runno != cdb->GetRun())
 418       cdb->SetRun(runno);
 419     AliGeomManager::LoadGeometry();
 420   }
 421
 422   if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
 423     Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 424     if (fEsdEv) {
 425       for (Int_t i=0; i<nsm; ++i)
 426         fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 427     } else {
 428       for (Int_t i=0; i<nsm; ++i)
 429         fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 430     }
 431     fIsGeoMatsSet = kTRUE;
 432   }
 433
 434   if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 435     if (fEsdEv) {
 436       const AliESDRun *erun = fEsdEv->GetESDRun();
 437       AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 438                                                erun->GetCurrentDip(),
 439                                                AliMagF::kConvLHC,
 440                                                kFALSE,
 441                                                erun->GetBeamEnergy(),
 442                                                erun->GetBeamType());
 443       TGeoGlobalMagField::Instance()->SetField(field);
 444     } else {
 445       Double_t pol = -1; //polarity
 446       Double_t be = -1;  //beam energy
 447       AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 448       Int_t runno = fAodEv->GetRunNumber();
 449       if (runno>=136851 && runno<138275) {
 450         pol = -1;
 451         be = 2760;
 452         btype = AliMagF::kBeamTypeAA;
 453       } else if (runno>=138275 && runno<=139517) {
 454         pol = +1;
 455         be = 2760;
 456         btype = AliMagF::kBeamTypeAA;
 457       } else {
 458         AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 459       }
 460       TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 461     }
 462   }
 463
 464   Int_t cut = 1;
 465   fHCuts->Fill(cut++);
 466
 467   TString trgclasses;
 468   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 469   if (h) {
 470     trgclasses = fEsdEv->GetFiredTriggerClasses();
 471   } else {
 472     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 473     if (h2)
 474       trgclasses = h2->GetFiredTriggerClasses();
 475   }
 476   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 477     const char *name = fTrClassNamesArr->At(i)->GetName();
 478     if (trgclasses.Contains(name))
 479       fHTclsBeforeCuts->Fill(1+i);
 480   }
 481
 482   UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 483   if (res==0)
 484     return;
 485
 486   fHCuts->Fill(cut++);
 487
 488   const AliCentrality *centP = InputEvent()->GetCentrality();
 489   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 490   fHCent->Fill(cent);
 491   if (cent<fCentFrom||cent>fCentTo)
 492     return;
 493
 494   fHCuts->Fill(cut++);
 495
 496   if (fUseQualFlag) {
 497     if (centP->GetQuality()>0)
 498       return;
 499   }
 500
 501   fHCentQual->Fill(cent);
 502   fHCuts->Fill(cut++);
 503
 504   if (fEsdEv) {
 505     am->LoadBranch("PrimaryVertex.");
 506     am->LoadBranch("SPDVertex.");
 507     am->LoadBranch("TPCVertex.");
 508   } else {
 509     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 510     am->LoadBranch("vertices");
 511     if (!fAodEv) return;
 512   }
 513
 514   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 515   if (!vertex)
 516     return;
 517
 518   fHVertexZ->Fill(vertex->GetZ());
 519
 520   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 521     return;
 522
 523   fHCuts->Fill(cut++);
 524   fHVertexZ2->Fill(vertex->GetZ());
 525
 526   // count number of accepted events
 527   ++fNEvs;
 528
 529   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 530     const char *name = fTrClassNamesArr->At(i)->GetName();
 531     if (trgclasses.Contains(name))
 532       fHTclsAfterCuts->Fill(1+i);
 533   }
 534
 535   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 536   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
 537   fEsdCells    = 0; // will be set if ESD input used
 538   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
 539                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 540   fAodCells    = 0; // will be set if AOD input used
 541
 542   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 543   Bool_t clusattached = 0;
 544   Bool_t recalibrated = 0;
 545   if (1 && !fClusName.IsNull()) {
 546     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 547     TObjArray *ts = am->GetTasks();
 548     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 549     if (cltask && cltask->GetClusters()) {
 550       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 551       clusattached = cltask->GetAttachClusters();
 552       if (cltask->GetCalibData()!=0)
 553         recalibrated = kTRUE;
 554     }
 555   }
 556   if (1 && AODEvent() && !fClusName.IsNull()) {
 557     TList *l = AODEvent()->GetList();
 558     TClonesArray *clus = 0;
 559     if (l) {
 560       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 561       fAodClusters = clus;
 562     }
 563   }
 564
 565   if (fEsdEv) { // ESD input mode
 566     if (1 && (!fRecPoints||clusattached)) {
 567       if (!clusattached)
 568         am->LoadBranch("CaloClusters");
 569       TList *l = fEsdEv->GetList();
 570       TClonesArray *clus = 0;
 571       if (l) {
 572         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 573         fEsdClusters = clus;
 574       }
 575     }
 576     if (1) {
 577       if (!recalibrated)
 578         am->LoadBranch("EMCALCells.");
 579       fEsdCells = fEsdEv->GetEMCALCells();
 580     }
 581   } else if (fAodEv) { // AOD input mode
 582     if (1 && (!fAodClusters || clusattached)) {
 583       if (!clusattached)
 584         am->LoadBranch("caloClusters");
 585       TList *l = fAodEv->GetList();
 586       TClonesArray *clus = 0;
 587       if (l) {
 588         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 589         fAodClusters = clus;
 590       }
 591     }
 592     if (1) {
 593       if (!recalibrated)
 594         am->LoadBranch("emcalCells");
 595       fAodCells = fAodEv->GetEMCALCells();
 596     }
 597   } else {
 598     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 599   }
 600
 601   if (1) {
 602     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 603     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 604     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 605     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 606     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 607   }
 608
 609   if (fDoAfterburner)
 610     ClusterAfterburner();
 611
 612   CalcCaloTriggers();
 613   CalcPrimTracks();
 614   CalcTracks();
 615   CalcClusterProps();
 616
 617   FillCellHists();
 618   if (!fTrainMode) {
 619     FillClusHists();
 620     FillPionHists();
 621     FillOtherHists();
 622   }
 623   FillNtuple();
 624
 625   if (fTrainMode) {
 626     fSelTracks->Clear();
 627     fSelPrimTracks->Clear();
 628     fTriggers->Clear();
 629   }
 630
 631   PostData(1, fOutput);
 632 }
 633
 634 //________________________________________________________________________
 635 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 636 {
 637   // Terminate called at the end of analysis.
 638
 639   if (fNtuple) {
 640     TFile *f = OpenFile(1);
 641     if (f)
 642       fNtuple->Write();
 643   }
 644
 645   AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
 646 }
 647
 648 //________________________________________________________________________
 649 void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
 650 {
 651   // Calculate triggers
 652
 653   fTriggers->Clear();
 654
 655   AliVCaloCells *cells = fEsdCells;
 656   if (!cells)
 657     cells = fAodCells;
 658   if (!cells)
 659     return;
 660
 661   Int_t ncells = cells->GetNumberOfCells();
 662   if (ncells<=0)
 663     return;
 664
 665   if (ncells>fNAmpInTrigger) {
 666     delete [] fAmpInTrigger;
 667     fAmpInTrigger = new Float_t[ncells];
 668     fNAmpInTrigger = ncells;
 669   }
 670   for (Int_t i=0;i<ncells;++i)
 671     fAmpInTrigger[i] = 0;
 672
 673   std::map<Short_t,Short_t> map;
 674   for (Short_t pos=0;pos<ncells;++pos) {
 675     Short_t id = cells->GetCellNumber(pos);
 676     map[id]=pos;
 677   }
 678
 679   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 680   am->LoadBranch("EMCALTrigger.");
 681
 682   AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
 683   if (!triggers)
 684     return;
 685   if (triggers->GetEntries()<=0)
 686     return;
 687
 688   triggers->Reset();
 689   Int_t ntrigs=0;
 690   while (triggers->Next()) {
 691     Int_t gCol=0, gRow=0, ntimes=0;
 692     triggers->GetPosition(gCol,gRow);
 693     triggers->GetNL0Times(ntimes);
 694     if (ntimes<1)
 695       continue;
 696     Float_t amp=0;
 697     triggers->GetAmplitude(amp);
 698     Int_t find = -1;
 699     fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
 700     if (find<0)
 701       continue;
 702     Int_t cidx[4] = {-1};
 703     Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
 704     if (!ret)
 705       continue;
 706     Int_t trgtimes[25];
 707     triggers->GetL0Times(trgtimes);
 708     Int_t mintime = trgtimes[0];
 709     Int_t maxtime = trgtimes[0];
 710     Bool_t trigInTimeWindow = 0;
 711     for (Int_t i=0;i<ntimes;++i) {
 712       if (trgtimes[i]<mintime)
 713         mintime = trgtimes[i];
 714       if (maxtime<trgtimes[i])
 715         maxtime = trgtimes[i];
 716       if ((fMinL0Time<trgtimes[i]) && (fMaxL0Time>trgtimes[i]))
 717         trigInTimeWindow = 1;
 718     }
 719
 720     Double_t tenergy = 0;
 721     Double_t tphi=0;
 722     Double_t teta=0;
 723     for (Int_t i=0;i<3;++i) {
 724       Short_t pos = -1;
 725       std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
 726       if (it!=map.end())
 727         pos = it->second;
 728       if (pos<0)
 729         continue;
 730       if (trigInTimeWindow)
 731         fAmpInTrigger[pos] = amp; // save for usage in CalcClusProperties
 732       Float_t eta=-1, phi=-1;
 733       fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
 734       Double_t en= cells->GetAmplitude(pos);
 735       tenergy+=en;
 736       teta+=eta*en;
 737       tphi+=phi*en;
 738     }
 739     teta/=tenergy;
 740     tphi/=tenergy;
 741     AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
 742     trignew->fE       = tenergy;
 743     trignew->fEta     = teta;
 744     trignew->fPhi     = tphi;
 745     trignew->fAmp     = amp;
 746     trignew->fMinTime = mintime;
 747     trignew->fMaxTime = maxtime;
 748   }
 749 }
 750
 751 //________________________________________________________________________
 752 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 753 {
 754   // Calculate cluster properties
 755
 756   fClusters->Clear();
 757
 758   AliVCaloCells *cells = fEsdCells;
 759   if (!cells)
 760     cells = fAodCells;
 761   if (!cells)
 762     return;
 763
 764   TObjArray *clusters = fEsdClusters;
 765   if (!clusters)
 766     clusters = fAodClusters;
 767   if (!clusters)
 768     return;
 769
 770   Int_t ncells = cells->GetNumberOfCells();
 771   Int_t nclus  = clusters->GetEntries();
 772   Int_t ntrks  = fSelTracks->GetEntries();
 773   Bool_t btracks[6][ntrks];
 774   memset(btracks,0,sizeof(btracks));
 775
 776   std::map<Short_t,Short_t> map;
 777   for (Short_t pos=0;pos<ncells;++pos) {
 778     Short_t id = cells->GetCellNumber(pos);
 779     map[id]=pos;
 780   }
 781
 782   for(Int_t i=0, ncl=0; i<nclus; ++i) {
 783     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 784     if (!clus)
 785       continue;
 786     if (!clus->IsEMCAL())
 787       continue;
 788     if (clus->E()<fMinE)
 789       continue;
 790
 791     Float_t clsPos[3] = {0};
 792     clus->GetPosition(clsPos);
 793     TVector3 clsVec(clsPos);
 794     Double_t vertex[3] = {0};
 795     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 796     TLorentzVector clusterVec;
 797     clus->GetMomentum(clusterVec,vertex);
 798     Double_t clsEta = clusterVec.Eta();
 799
 800     AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
 801     cl->fE        = clus->E();
 802     cl->fR        = clsVec.Perp();
 803     cl->fEta      = clsVec.Eta();
 804     cl->fPhi      = clsVec.Phi();
 805     cl->fN        = clus->GetNCells();
 806     cl->fN1       = GetNCells(clus,0.100);
 807     cl->fN3       = GetNCells(clus,0.300);
 808     Short_t id    = -1;
 809     Double_t emax = GetMaxCellEnergy(clus, id);
 810     cl->fIdMax    = id;
 811     cl->fEmax     = emax;
 812     if (clus->GetDistanceToBadChannel()<10000)
 813       cl->fDbc    = clus->GetDistanceToBadChannel();
 814     if (!TMath::IsNaN(clus->GetDispersion()))
 815       cl->fDisp   = clus->GetDispersion();
 816     if (!TMath::IsNaN(clus->GetM20()))
 817       cl->fM20      = clus->GetM20();
 818     if (!TMath::IsNaN(clus->GetM02()))
 819       cl->fM02      = clus->GetM02();
 820     Double_t maxAxis = 0, minAxis = 0;
 821     GetSigma(clus,maxAxis,minAxis);
 822     clus->SetTOF(maxAxis);     // store sigma in TOF
 823     cl->fSig      = maxAxis;
 824     Double_t clusterEcc = 0;
 825     if (maxAxis > 0)
 826       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 827     clus->SetChi2(clusterEcc); // store ecc in chi2
 828     cl->fEcc      = clusterEcc;
 829     cl->fTrIso    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
 830     cl->fTrIso1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
 831     cl->fTrIso2   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
 832     cl->fCeCore   = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
 833     cl->fCeIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
 834     Double_t trigpen = 0;
 835     Double_t trignen = 0;
 836     for(Int_t j=0; j<cl->fN; ++j) {
 837       Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
 838       Short_t pos = -1;
 839       std::map<Short_t,Short_t>::iterator it = map.find(cid);
 840       if (it!=map.end())
 841         pos = it->second;
 842       if (pos<0)
 843         continue;
 844       if (fAmpInTrigger[pos]>0)
 845         trigpen += cells->GetAmplitude(pos);
 846       else if (fAmpInTrigger[pos]<0)
 847         trignen += cells->GetAmplitude(pos);
 848     }
 849     if (trigpen>0) {
 850       cl->fTrigM = 1;
 851       cl->fTrigE = trigpen;
 852     }
 853     if (trignen>0) {
 854       cl->fTrigM     = 1;
 855       cl->fTrigMaskE = trignen;
 856     }
 857
 858     // track matching
 859     Double_t mind2 = 1e10;
 860     for(Int_t j = 0; j<ntrks; ++j) {
 861       AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
 862       if (!track)
 863         continue;
 864
 865       if (TMath::Abs(clsEta-track->Eta())>0.5)
 866         continue;
 867
 868       TVector3 vec(clsPos);
 869       Int_t index =  (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
 870       if (btracks[index-4][j]) {
 871         continue;
 872       }
 873
 874       Float_t tmpR=-1, tmpZ=-1;
 875       if (!fDoTrMatGeom) {
 876         AliExternalTrackParam *tParam = 0;
 877         if (fEsdEv) {
 878           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
 879           tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
 880         } else
 881           tParam = new AliExternalTrackParam(track);
 882
 883         Double_t bfield[3] = {0};
 884         track->GetBxByBz(bfield);
 885         Double_t alpha = (index+0.5)*20*TMath::DegToRad();
 886         vec.RotateZ(-alpha);   //Rotate the cluster to the local extrapolation coordinate system
 887         tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
 888         Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
 889         if (!ret) {
 890           btracks[index-4][j]=1;
 891           delete tParam;
 892           continue;
 893         }
 894         Double_t trkPos[3] = {0};
 895         tParam->GetXYZ(trkPos); //Get the extrapolated global position
 896         tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
 897         tmpZ = clsPos[2]-trkPos[2];
 898         delete tParam;
 899       } else {
 900         if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
 901           continue;
 902         AliExternalTrackParam tParam(track);
 903         if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
 904           continue;
 905       }
 906
 907       Double_t d2 = tmpR;
 908       if (mind2>d2) {
 909         mind2=d2;
 910         cl->fTrDz   = tmpZ;
 911         cl->fTrDr   = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
 912         cl->fTrEp   = clus->E()/track->P();
 913         cl->fTrackM = 1;
 914       }
 915     }
 916
 917     if (cl->fTrackM) {
 918       fHMatchDr->Fill(cl->fTrDr);
 919       fHMatchDz->Fill(cl->fTrDz);
 920       fHMatchEp->Fill(cl->fTrEp);
 921     }
 922   }
 923 }
 924
 925 //________________________________________________________________________
 926 void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
 927 {
 928   // Calculate track properties.
 929
 930   fSelPrimTracks->Clear();
 931
 932   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 933   TClonesArray *tracks = 0;
 934   if (fEsdEv) {
 935     am->LoadBranch("Tracks");
 936     TList *l = fEsdEv->GetList();
 937     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
 938   } else {
 939     am->LoadBranch("tracks");
 940     TList *l = fAodEv->GetList();
 941     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
 942   }
 943
 944   if (!tracks)
 945     return;
 946
 947   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 948   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
 949   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
 950
 951   if (fEsdEv) {
 952     fSelPrimTracks->SetOwner(kTRUE);
 953     am->LoadBranch("PrimaryVertex.");
 954     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
 955     am->LoadBranch("SPDVertex.");
 956     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
 957     am->LoadBranch("Tracks");
 958     const Int_t Ntracks = tracks->GetEntries();
 959     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 960       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
 961       if (!track) {
 962         AliWarning(Form("Could not receive track %d\n", iTracks));
 963         continue;
 964       }
 965       if (fTrCuts && !fTrCuts->IsSelected(track))
 966         continue;
 967       Double_t eta = track->Eta();
 968       if (eta<-1||eta>1)
 969         continue;
 970       if (track->Phi()<phimin||track->Phi()>phimax)
 971         continue;
 972
 973       AliESDtrack copyt(*track);
 974       Double_t bfield[3];
 975       copyt.GetBxByBz(bfield);
 976       AliExternalTrackParam tParam;
 977       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
 978       if (!relate)
 979         continue;
 980       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
 981
 982       Double_t p[3]      = { 0. };
 983       copyt.GetPxPyPz(p);
 984       Double_t pos[3]    = { 0. };
 985       copyt.GetXYZ(pos);
 986       Double_t covTr[21] = { 0. };
 987       copyt.GetCovarianceXYZPxPyPz(covTr);
 988       Double_t pid[10]   = { 0. };
 989       copyt.GetESDpid(pid);
 990       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
 991                                             copyt.GetLabel(),
 992                                             p,
 993                                             kTRUE,
 994                                             pos,
 995                                             kFALSE,
 996                                             covTr,
 997                                             (Short_t)copyt.GetSign(),
 998                                             copyt.GetITSClusterMap(),
 999                                             pid,
1000                                             0,/*fPrimaryVertex,*/
1001                                             kTRUE, // check if this is right
1002                                             vtx->UsesTrack(copyt.GetID()));
1003       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1004       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1005       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1006       if(ndf>0)
1007         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1008       else
1009         aTrack->SetChi2perNDF(-1);
1010       aTrack->SetFlags(copyt.GetStatus());
1011       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
1012       fSelPrimTracks->Add(aTrack);
1013     }
1014   } else {
1015     Int_t ntracks = tracks->GetEntries();
1016     for (Int_t i=0; i<ntracks; ++i) {
1017       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1018       if (!track)
1019         continue;
1020       Double_t eta = track->Eta();
1021       if (eta<-1||eta>1)
1022         continue;
1023       if (track->Phi()<phimin||track->Phi()>phimax)
1024         continue;
1025       if(track->GetTPCNcls()<fMinNClusPerTr)
1026         continue;
1027       //todo: Learn how to set/filter AODs for prim/sec tracks
1028       fSelPrimTracks->Add(track);
1029     }
1030   }
1031 }
1032
1033 //________________________________________________________________________
1034 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
1035 {
1036   // Calculate track properties (including secondaries).
1037
1038   fSelTracks->Clear();
1039
1040   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1041   TClonesArray *tracks = 0;
1042   if (fEsdEv) {
1043     am->LoadBranch("Tracks");
1044     TList *l = fEsdEv->GetList();
1045     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1046   } else {
1047     am->LoadBranch("tracks");
1048     TList *l = fAodEv->GetList();
1049     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1050   }
1051
1052   if (!tracks)
1053     return;
1054
1055   if (fEsdEv) {
1056     const Int_t Ntracks = tracks->GetEntries();
1057     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1058       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1059       if (!track) {
1060         AliWarning(Form("Could not receive track %d\n", iTracks));
1061         continue;
1062       }
1063       if (fTrCuts && !fTrCuts->IsSelected(track))
1064         continue;
1065       Double_t eta = track->Eta();
1066       if (eta<-1||eta>1)
1067         continue;
1068       fSelTracks->Add(track);
1069     }
1070   } else {
1071     Int_t ntracks = tracks->GetEntries();
1072     for (Int_t i=0; i<ntracks; ++i) {
1073       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1074       if (!track)
1075         continue;
1076       Double_t eta = track->Eta();
1077       if (eta<-1||eta>1)
1078         continue;
1079       if(track->GetTPCNcls()<fMinNClusPerTr)
1080         continue;
1081
1082       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
1083         AliExternalTrackParam tParam(track);
1084         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1085           Double_t trkPos[3];
1086           tParam.GetXYZ(trkPos);
1087           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1088         }
1089       }
1090       fSelTracks->Add(track);
1091     }
1092   }
1093 }
1094
1095 //________________________________________________________________________
1096 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
1097 {
1098   // Run custer reconstruction afterburner.
1099
1100   AliVCaloCells *cells = fEsdCells;
1101   if (!cells)
1102     cells = fAodCells;
1103
1104   if (!cells)
1105     return;
1106
1107   Int_t ncells = cells->GetNumberOfCells();
1108   if (ncells<=0)
1109     return;
1110
1111   Double_t cellMeanE = 0, cellSigE = 0;
1112   for (Int_t i = 0; i<ncells; ++i) {
1113     Double_t cellE = cells->GetAmplitude(i);
1114     cellMeanE += cellE;
1115     cellSigE += cellE*cellE;
1116   }
1117   cellMeanE /= ncells;
1118   cellSigE /= ncells;
1119   cellSigE -= cellMeanE*cellMeanE;
1120   if (cellSigE<0)
1121     cellSigE = 0;
1122   cellSigE = TMath::Sqrt(cellSigE / ncells);
1123
1124   Double_t subE = cellMeanE - 7*cellSigE;
1125   if (subE<0)
1126     return;
1127
1128   for (Int_t i = 0; i<ncells; ++i) {
1129     Short_t id=-1;
1130     Double_t amp=0,time=0;
1131     if (!cells->GetCell(i, id, amp, time))
1132       continue;
1133     amp -= cellMeanE;
1134     if (amp<0.001)
1135       amp = 0;
1136     cells->SetCell(i, id, amp, time);
1137   }
1138
1139   TObjArray *clusters = fEsdClusters;
1140   if (!clusters)
1141     clusters = fAodClusters;
1142   if (!clusters)
1143     return;
1144
1145   Int_t nclus = clusters->GetEntries();
1146   for (Int_t i = 0; i<nclus; ++i) {
1147     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1148     if (!clus->IsEMCAL())
1149       continue;
1150     Int_t nc = clus->GetNCells();
1151     Double_t clusE = 0;
1152     UShort_t ids[100] = {0};
1153     Double_t fra[100] = {0};
1154     for (Int_t j = 0; j<nc; ++j) {
1155       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1156       Double_t cen = cells->GetCellAmplitude(id);
1157       clusE += cen;
1158       if (cen>0) {
1159         ids[nc] = id;
1160         ++nc;
1161       }
1162     }
1163     if (clusE<=0) {
1164       clusters->RemoveAt(i);
1165       continue;
1166     }
1167
1168     for (Int_t j = 0; j<nc; ++j) {
1169       Short_t id = ids[j];
1170       Double_t cen = cells->GetCellAmplitude(id);
1171       fra[j] = cen/clusE;
1172     }
1173     clus->SetE(clusE);
1174     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1175     if (aodclus) {
1176       aodclus->Clear("");
1177       aodclus->SetNCells(nc);
1178       aodclus->SetCellsAmplitudeFraction(fra);
1179       aodclus->SetCellsAbsId(ids);
1180     }
1181   }
1182   clusters->Compress();
1183 }
1184
1185 //________________________________________________________________________
1186 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1187 {
1188   // Fill histograms related to cell properties.
1189
1190   AliVCaloCells *cells = fEsdCells;
1191   if (!cells)
1192     cells = fAodCells;
1193
1194   if (!cells)
1195     return;
1196
1197   Int_t cellModCount[12] = {0};
1198   Double_t cellMaxE = 0;
1199   Double_t cellMeanE = 0;
1200   Int_t ncells = cells->GetNumberOfCells();
1201   if (ncells==0)
1202     return;
1203
1204   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1205
1206   for (Int_t i = 0; i<ncells; ++i) {
1207     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1208     Double_t cellE = cells->GetAmplitude(i);
1209     fHCellE->Fill(cellE);
1210     if (cellE>cellMaxE)
1211       cellMaxE = cellE;
1212     cellMeanE += cellE;
1213
1214     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1215     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1216     if (!ret) {
1217       AliError(Form("Could not get cell index for %d", absID));
1218       continue;
1219     }
1220     ++cellModCount[iSM];
1221     Int_t iPhi=-1, iEta=-1;
1222     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1223     fHColuRow[iSM]->Fill(iEta,iPhi,1);
1224     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1225     fHCellFreqNoCut[iSM]->Fill(absID);
1226     if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1227     if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
1228     if (fHCellCheckE && fHCellCheckE[absID])
1229       fHCellCheckE[absID]->Fill(cellE);
1230     fHCellFreqE[iSM]->Fill(absID, cellE);
1231     if (1) {
1232       //todo
1233       TVector3 pos;
1234       fGeom->GetGlobal(absID,pos);
1235       Int_t id2 = fGeom-> GetAbsCellIdFromCellIndexes(iSM, iPhi+1,iEta+1);
1236       TVector3 pos2;
1237       fGeom->GetGlobal(id2,pos2);
1238       cout << "delta phi " << pos.DeltaPhi(pos2) << " " << TMath::Abs(pos2.Eta()-pos.Eta()) << endl;
1239     }
1240   }
1241   fHCellH->Fill(cellMaxE);
1242   cellMeanE /= ncells;
1243   fHCellM->Fill(cellMeanE);
1244   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1245   for (Int_t i=0; i<nsm; ++i)
1246     fHCellMult[i]->Fill(cellModCount[i]);
1247 }
1248
1249 //________________________________________________________________________
1250 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1251 {
1252   // Fill histograms related to cluster properties.
1253
1254   TObjArray *clusters = fEsdClusters;
1255   if (!clusters)
1256     clusters = fAodClusters;
1257   if (!clusters)
1258     return;
1259
1260   Double_t vertex[3] = {0};
1261   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1262
1263   Int_t nclus = clusters->GetEntries();
1264   for(Int_t i = 0; i<nclus; ++i) {
1265     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1266     if (!clus)
1267       continue;
1268     if (!clus->IsEMCAL())
1269       continue;
1270     TLorentzVector clusterVec;
1271     clus->GetMomentum(clusterVec,vertex);
1272     Double_t maxAxis    = clus->GetTOF(); //sigma
1273     Double_t clusterEcc = clus->Chi2();   //eccentricity
1274     fHClustEccentricity->Fill(clusterEcc);
1275     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1276     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1277     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1278     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1279     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1280   }
1281 }
1282
1283 //________________________________________________________________________
1284 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1285 {
1286   // Fill ntuple.
1287
1288   if (!fNtuple)
1289     return;
1290
1291   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1292   if (fAodEv) {
1293     fHeader->fRun            = fAodEv->GetRunNumber();
1294     fHeader->fOrbit          = fAodEv->GetHeader()->GetOrbitNumber();
1295     fHeader->fPeriod         = fAodEv->GetHeader()->GetPeriodNumber();
1296     fHeader->fBx             = fAodEv->GetHeader()->GetBunchCrossNumber();
1297     fHeader->fL0             = fAodEv->GetHeader()->GetL0TriggerInputs();
1298     fHeader->fL1             = fAodEv->GetHeader()->GetL1TriggerInputs();
1299     fHeader->fL2             = fAodEv->GetHeader()->GetL2TriggerInputs();
1300     fHeader->fTrClassMask    = fAodEv->GetHeader()->GetTriggerMask();
1301     fHeader->fTrCluster      = fAodEv->GetHeader()->GetTriggerCluster();
1302     fHeader->fOffTriggers    = fAodEv->GetHeader()->GetOfflineTrigger();
1303     fHeader->fFiredTriggers  = fAodEv->GetHeader()->GetFiredTriggerClasses();
1304   } else {
1305     fHeader->fRun            = fEsdEv->GetRunNumber();
1306     fHeader->fOrbit          = fEsdEv->GetHeader()->GetOrbitNumber();
1307     fHeader->fPeriod         = fEsdEv->GetESDRun()->GetPeriodNumber();
1308     fHeader->fBx             = fEsdEv->GetHeader()->GetBunchCrossNumber();
1309     fHeader->fL0             = fEsdEv->GetHeader()->GetL0TriggerInputs();
1310     fHeader->fL1             = fEsdEv->GetHeader()->GetL1TriggerInputs();
1311     fHeader->fL2             = fEsdEv->GetHeader()->GetL2TriggerInputs();
1312     fHeader->fTrClassMask    = fEsdEv->GetHeader()->GetTriggerMask();
1313     fHeader->fTrCluster      = fEsdEv->GetHeader()->GetTriggerCluster();
1314     fHeader->fOffTriggers    = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1315     fHeader->fFiredTriggers  = fEsdEv->GetFiredTriggerClasses();
1316   }
1317   AliCentrality *cent = InputEvent()->GetCentrality();
1318   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
1319   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
1320   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
1321   fHeader->fCqual     = cent->GetQuality();
1322
1323   AliEventplane *ep = InputEvent()->GetEventplane();
1324   if (ep) {
1325     if (ep->GetQVector())
1326       fHeader->fPsi     = ep->GetQVector()->Phi()/2. ;
1327     fHeader->fPsi = -1;
1328     if (ep->GetQsub1()&&ep->GetQsub2())
1329       fHeader->fPsiRes  = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1330     else
1331       fHeader->fPsiRes = 0;
1332   }
1333
1334   Double_t val = 0;
1335   TString trgclasses(fHeader->fFiredTriggers);
1336   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1337     const char *name = fTrClassNamesArr->At(j)->GetName();
1338     if (trgclasses.Contains(name))
1339       val += TMath::Power(2,j);
1340   }
1341   fHeader->fTcls = (UInt_t)val;
1342
1343   if (fAodEv) {
1344     am->LoadBranch("vertices");
1345     AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1346     FillVertex(fPrimVert, pv);
1347     AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1348     FillVertex(fSpdVert, sv);
1349   } else {
1350     am->LoadBranch("PrimaryVertex.");
1351     const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1352     FillVertex(fPrimVert, pv);
1353     am->LoadBranch("SPDVertex.");
1354     const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1355     FillVertex(fSpdVert, sv);
1356     am->LoadBranch("TPCVertex.");
1357     const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1358     FillVertex(fTpcVert, tv);
1359   }
1360
1361   fNtuple->Fill();
1362 }
1363
1364 //________________________________________________________________________
1365 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1366 {
1367   // Fill histograms related to pions.
1368
1369   Double_t vertex[3] = {0};
1370   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1371
1372   TObjArray *clusters = fEsdClusters;
1373   if (!clusters)
1374     clusters = fAodClusters;
1375
1376   if (clusters) {
1377     TLorentzVector clusterVec1;
1378     TLorentzVector clusterVec2;
1379     TLorentzVector pionVec;
1380
1381     Int_t nclus = clusters->GetEntries();
1382     for (Int_t i = 0; i<nclus; ++i) {
1383       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1384       if (!clus1)
1385         continue;
1386       if (!clus1->IsEMCAL())
1387         continue;
1388       if (clus1->E()<fMinE)
1389         continue;
1390       if (clus1->GetNCells()<fNminCells)
1391         continue;
1392       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1393         continue;
1394       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1395         continue;
1396       clus1->GetMomentum(clusterVec1,vertex);
1397       for (Int_t j = i+1; j<nclus; ++j) {
1398         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1399         if (!clus2)
1400           continue;
1401         if (!clus2->IsEMCAL())
1402           continue;
1403         if (clus2->E()<fMinE)
1404           continue;
1405         if (clus2->GetNCells()<fNminCells)
1406           continue;
1407         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1408           continue;
1409         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1410           continue;
1411         clus2->GetMomentum(clusterVec2,vertex);
1412         pionVec = clusterVec1 + clusterVec2;
1413         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1414         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1415         if (pionZgg < fAsymMax) {
1416           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1417           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1418           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1419           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1420           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1421           fHPionInvMasses[bin]->Fill(pionVec.M());
1422         }
1423       }
1424     }
1425   }
1426 }
1427
1428 //________________________________________________________________________
1429 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1430 {
1431   // Fill other histograms.
1432 }
1433
1434 //__________________________________________________________________________________________________
1435 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1436 {
1437   // Fill vertex from ESD vertex info.
1438
1439   v->fVx   = esdv->GetX();
1440   v->fVy   = esdv->GetY();
1441   v->fVz   = esdv->GetZ();
1442   v->fVc   = esdv->GetNContributors();
1443   v->fDisp = esdv->GetDispersion();
1444   v->fZres = esdv->GetZRes();
1445   v->fChi2 = esdv->GetChi2();
1446   v->fSt   = esdv->GetStatus();
1447   v->fIs3D = esdv->IsFromVertexer3D();
1448   v->fIsZ  = esdv->IsFromVertexerZ();
1449 }
1450
1451 //__________________________________________________________________________________________________
1452 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1453 {
1454   // Fill vertex from AOD vertex info.
1455
1456   v->fVx   = aodv->GetX();
1457   v->fVy   = aodv->GetY();
1458   v->fVz   = aodv->GetZ();
1459   v->fVc   = aodv->GetNContributors();
1460   v->fChi2 = aodv->GetChi2();
1461 }
1462
1463 //________________________________________________________________________
1464 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1465 {
1466   // Compute isolation based on cell content.
1467
1468   AliVCaloCells *cells = fEsdCells;
1469   if (!cells)
1470     cells = fAodCells;
1471   if (!cells)
1472     return 0;
1473
1474   Double_t cellIsolation = 0;
1475   Double_t rad2 = radius*radius;
1476   Int_t ncells = cells->GetNumberOfCells();
1477   for (Int_t i = 0; i<ncells; ++i) {
1478     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1479     Double_t cellE = cells->GetAmplitude(i);
1480     Float_t eta, phi;
1481     fGeom->EtaPhiFromIndex(absID,eta,phi);
1482     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1483     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1484     if(dist>rad2)
1485       continue;
1486     cellIsolation += cellE;
1487   }
1488   return cellIsolation;
1489 }
1490
1491 //________________________________________________________________________
1492 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster, Short_t &id) const
1493 {
1494   // Get maximum energy of attached cell.
1495
1496   id = -1;
1497   Double_t maxe = 0;
1498   Int_t ncells = cluster->GetNCells();
1499   if (fEsdCells) {
1500     for (Int_t i=0; i<ncells; i++) {
1501       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1502       if (e>maxe) {
1503         maxe = e;
1504         id   = cluster->GetCellAbsId(i);
1505       }
1506     }
1507   } else {
1508     for (Int_t i=0; i<ncells; i++) {
1509       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1510       if (e>maxe)
1511         maxe = e;
1512         id   = cluster->GetCellAbsId(i);
1513     }
1514   }
1515   return maxe;
1516 }
1517
1518 //________________________________________________________________________
1519 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
1520 {
1521   // Calculate the (E) weighted variance along the longer (eigen) axis.
1522
1523   sigmaMax = 0;          // cluster variance along its longer axis
1524   sigmaMin = 0;          // cluster variance along its shorter axis
1525   Double_t Ec  = c->E(); // cluster energy
1526   if(Ec<=0)
1527     return;
1528   Double_t Xc  = 0 ;     // cluster first moment along X
1529   Double_t Yc  = 0 ;     // cluster first moment along Y
1530   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
1531   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
1532   Double_t Syy = 0 ;     // cluster covariance^2
1533
1534   AliVCaloCells *cells = fEsdCells;
1535   if (!cells)
1536     cells = fAodCells;
1537
1538   if (!cells)
1539     return;
1540
1541   Int_t ncells = c->GetNCells();
1542   if (ncells==1)
1543     return;
1544
1545   for(Int_t j=0; j<ncells; ++j) {
1546     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1547     Double_t cellen = cells->GetCellAmplitude(id);
1548     TVector3 pos;
1549     fGeom->GetGlobal(id,pos);
1550     Xc  += cellen*pos.X();
1551     Yc  += cellen*pos.Y();
1552     Sxx += cellen*pos.X()*pos.X();
1553     Syy += cellen*pos.Y()*pos.Y();
1554     Sxy += cellen*pos.X()*pos.Y();
1555   }
1556   Xc  /= Ec;
1557   Yc  /= Ec;
1558   Sxx /= Ec;
1559   Syy /= Ec;
1560   Sxy /= Ec;
1561   Sxx -= Xc*Xc;
1562   Syy -= Yc*Yc;
1563   Sxy -= Xc*Yc;
1564   Sxx = TMath::Abs(Sxx);
1565   Syy = TMath::Abs(Syy);
1566   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1567   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1568   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1569   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
1570 }
1571
1572 //________________________________________________________________________
1573 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
1574 {
1575   // Calculate number of attached cells above emin.
1576
1577   AliVCaloCells *cells = fEsdCells;
1578   if (!cells)
1579     cells = fAodCells;
1580   if (!cells)
1581     return 0;
1582
1583   Int_t n = 0;
1584   Int_t ncells = c->GetNCells();
1585   for(Int_t j=0; j<ncells; ++j) {
1586     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1587     Double_t cellen = cells->GetCellAmplitude(id);
1588     if (cellen>=emin)
1589       ++n;
1590   }
1591   return n;
1592 }
1593
1594 //________________________________________________________________________
1595 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1596 {
1597   // Compute isolation based on tracks.
1598
1599   Double_t trkIsolation = 0;
1600   Double_t rad2 = radius*radius;
1601   Int_t ntrks = fSelPrimTracks->GetEntries();
1602   for(Int_t j = 0; j<ntrks; ++j) {
1603     AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1604     if (!track)
1605       continue;
1606     if (track->Pt()<pt)
1607       continue;
1608     Float_t eta = track->Eta();
1609     Float_t phi = track->Phi();
1610     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1611     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1612     if(dist>rad2)
1613       continue;
1614     trkIsolation += track->Pt();
1615   }
1616   return trkIsolation;
1617 }
1618
1619