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