PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0PbPb.cxx

   1 // $Id$
   2 //
   3 // Analysis task for neutral pions (into two gammas).
   4 //
   5 // Author: C.Loizides, E.Braidot
   6
   7 #include "AliAnalysisTaskEMCALPi0PbPb.h"
   8 #include <TAxis.h>
   9 #include <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TFile.h>
  12 #include <TGeoGlobalMagField.h>
  13 #include <TH1F.h>
  14 #include <TH2F.h>
  15 #include <TList.h>
  16 #include <TLorentzVector.h>
  17 #include <TNtuple.h>
  18 #include <TProfile.h>
  19 #include <TRegexp.h>
  20 #include <TString.h>
  21 #include <TVector2.h>
  22 #include "AliAODEvent.h"
  23 #include "AliAODMCParticle.h"
  24 #include "AliAODVertex.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  27 #include "AliCDBManager.h"
  28 #include "AliCentrality.h"
  29 #include "AliEMCALDigit.h"
  30 #include "AliEMCALGeometry.h"
  31 #include "AliEMCALRecPoint.h"
  32 #include "AliEMCALRecoUtils.h"
  33 #include "AliESDCaloTrigger.h"
  34 #include "AliESDEvent.h"
  35 #include "AliESDUtils.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDtrackCuts.h"
  39 #include "AliEventplane.h"
  40 #include "AliGeomManager.h"
  41 #include "AliInputEventHandler.h"
  42 #include "AliLog.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCParticle.h"
  45 #include "AliMagF.h"
  46 #include "AliMultiplicity.h"
  47 #include "AliStack.h"
  48 #include "AliTrackerBase.h"
  49 #include "AliTriggerAnalysis.h"
  50
  51 using std::cout;
  52 using std::endl;
  53 using std::max;
  54
  55 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  56
  57 //________________________________________________________________________
  58 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
  59   : AliAnalysisTaskSE(),
  60     fCentVar("V0M"),
  61     fCentFrom(0),
  62     fCentTo(100),
  63     fVtxZMin(-10),
  64     fVtxZMax(+10),
  65     fUseQualFlag(1),
  66     fClusName(),
  67     fDoNtuple(0),
  68     fDoAfterburner(0),
  69     fAsymMax(1),
  70     fNminCells(1),
  71     fMinE(0.100),
  72     fMinErat(0),
  73     fMinEcc(0),
  74     fGeoName("EMCAL_FIRSTYEARV1"),
  75     fMinNClusPerTr(50),
  76     fIsoDist(0.2),
  77     fTrClassNames(""),
  78     fTrCuts(0),
  79     fPrimTrCuts(0),
  80     fPrimTracksName(""),
  81     fDoTrMatGeom(0),
  82     fTrainMode(0),
  83     fMarkCells(),
  84     fMinL0Time(-1),
  85     fMaxL0Time(1024),
  86     fMcMode(0),
  87     fEmbedMode(0),
  88     fGeom(0),
  89     fReco(0),
  90     fTrigName(),
  91     fDoPSel(kTRUE),
  92     fIsGeoMatsSet(0),
  93     fNEvs(0),
  94     fOutput(0),
  95     fTrClassNamesArr(0),
  96     fEsdEv(0),
  97     fAodEv(0),
  98     fRecPoints(0),
  99     fDigits(0),
 100     fEsdClusters(0),
 101     fEsdCells(0),
 102     fAodClusters(0),
 103     fAodCells(0),
 104     fPtRanges(0),
 105     fSelTracks(0),
 106     fSelPrimTracks(0),
 107     fNtuple(0),
 108     fHeader(0),
 109     fPrimVert(0),
 110     fSpdVert(0),
 111     fTpcVert(0),
 112     fClusters(0),
 113     fTriggers(0),
 114     fMcParts(0),
 115     fHCuts(0x0),
 116     fHVertexZ(0x0),
 117     fHVertexZ2(0x0),
 118     fHCent(0x0),
 119     fHCentQual(0x0),
 120     fHTclsBeforeCuts(0x0),
 121     fHTclsAfterCuts(0x0),
 122     fHColuRow(0x0),
 123     fHColuRowE(0x0),
 124     fHCellMult(0x0),
 125     fHCellE(0x0),
 126     fHCellH(0x0),
 127     fHCellM(0x0),
 128     fHCellM2(0x0),
 129     fHCellFreqNoCut(0x0),
 130     fHCellFreqCut100M(0x0),
 131     fHCellFreqCut300M(0x0),
 132     fHCellFreqE(0x0),
 133     fHCellCheckE(0x0),
 134     fHClustEccentricity(0),
 135     fHClustEtaPhi(0x0),
 136     fHClustEnergyPt(0x0),
 137     fHClustEnergySigma(0x0),
 138     fHClustSigmaSigma(0x0),
 139     fHClustNCellEnergyRatio(0x0),
 140     fHClustEnergyNCell(0x0),
 141     fHPrimTrackPt(0x0),
 142     fHPrimTrackEta(0x0),
 143     fHPrimTrackPhi(0x0),
 144     fHMatchDr(0x0),
 145     fHMatchDz(0x0),
 146     fHMatchEp(0x0),
 147     fHPionEtaPhi(0x0),
 148     fHPionMggPt(0x0),
 149     fHPionMggAsym(0x0),
 150     fHPionMggDgg(0x0),
 151     fHPionInvMasses()
 152 {
 153   // Constructor.
 154 }
 155
 156 //________________________________________________________________________
 157 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
 158   : AliAnalysisTaskSE(name),
 159     fCentVar("V0M"),
 160     fCentFrom(0),
 161     fCentTo(100),
 162     fVtxZMin(-10),
 163     fVtxZMax(+10),
 164     fUseQualFlag(1),
 165     fClusName(),
 166     fDoNtuple(0),
 167     fDoAfterburner(0),
 168     fAsymMax(1),
 169     fNminCells(1),
 170     fMinE(0.100),
 171     fMinErat(0),
 172     fMinEcc(0),
 173     fGeoName("EMCAL_FIRSTYEARV1"),
 174     fMinNClusPerTr(50),
 175     fIsoDist(0.2),
 176     fTrClassNames(""),
 177     fTrCuts(0),
 178     fPrimTrCuts(0),
 179     fPrimTracksName(""),
 180     fDoTrMatGeom(0),
 181     fTrainMode(0),
 182     fMarkCells(),
 183     fMinL0Time(-1),
 184     fMaxL0Time(1024),
 185     fMcMode(0),
 186     fEmbedMode(0),
 187     fGeom(0),
 188     fReco(0),
 189     fTrigName(),
 190     fDoPSel(kTRUE),
 191     fIsGeoMatsSet(0),
 192     fNEvs(0),
 193     fOutput(0),
 194     fTrClassNamesArr(0),
 195     fEsdEv(0),
 196     fAodEv(0),
 197     fRecPoints(0),
 198     fDigits(0),
 199     fEsdClusters(0),
 200     fEsdCells(0),
 201     fAodClusters(0),
 202     fAodCells(0),
 203     fPtRanges(0),
 204     fSelTracks(0),
 205     fSelPrimTracks(0),
 206     fNtuple(0),
 207     fHeader(0),
 208     fPrimVert(0),
 209     fSpdVert(0),
 210     fTpcVert(0),
 211     fClusters(0),
 212     fTriggers(0),
 213     fMcParts(0),
 214     fHCuts(0x0),
 215     fHVertexZ(0x0),
 216     fHVertexZ2(0x0),
 217     fHCent(0x0),
 218     fHCentQual(0x0),
 219     fHTclsBeforeCuts(0x0),
 220     fHTclsAfterCuts(0x0),
 221     fHColuRow(0x0),
 222     fHColuRowE(0x0),
 223     fHCellMult(0x0),
 224     fHCellE(0x0),
 225     fHCellH(0x0),
 226     fHCellM(0x0),
 227     fHCellM2(0x0),
 228     fHCellFreqNoCut(0x0),
 229     fHCellFreqCut100M(0x0),
 230     fHCellFreqCut300M(0x0),
 231     fHCellFreqE(0x0),
 232     fHCellCheckE(0x0),
 233     fHClustEccentricity(0),
 234     fHClustEtaPhi(0x0),
 235     fHClustEnergyPt(0x0),
 236     fHClustEnergySigma(0x0),
 237     fHClustSigmaSigma(0x0),
 238     fHClustNCellEnergyRatio(0x0),
 239     fHClustEnergyNCell(0x0),
 240     fHPrimTrackPt(0x0),
 241     fHPrimTrackEta(0x0),
 242     fHPrimTrackPhi(0x0),
 243     fHMatchDr(0x0),
 244     fHMatchDz(0x0),
 245     fHMatchEp(0x0),
 246     fHPionEtaPhi(0x0),
 247     fHPionMggPt(0x0),
 248     fHPionMggAsym(0x0),
 249     fHPionMggDgg(0x0),
 250     fHPionInvMasses()
 251 {
 252   // Constructor.
 253
 254   DefineOutput(1, TList::Class());
 255   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks,EMCALTrigger.,SPDPileupVertices,TrkPileupVertices "
 256                "AOD:header,vertices,emcalCells,tracks";
 257 }
 258
 259 //________________________________________________________________________
 260 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
 261 {
 262   // Destructor.
 263
 264   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 265     delete fOutput; fOutput = 0;
 266   }
 267   delete fPtRanges; fPtRanges = 0;
 268   fGeom = 0; // do not delete geometry when using instance
 269   delete fReco; fReco = 0;
 270   delete fTrClassNamesArr;
 271   delete fSelTracks;
 272   delete fSelPrimTracks;
 273   delete [] fHColuRow;
 274   delete [] fHColuRowE;
 275   delete [] fHCellMult;
 276   delete [] fHCellFreqNoCut;
 277   delete [] fHCellFreqCut100M;
 278   delete [] fHCellFreqCut300M;
 279 }
 280
 281 //________________________________________________________________________
 282 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 283 {
 284   // Create user objects here.
 285
 286   cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
 287   cout << " fCentVar:       " << fCentVar << endl;
 288   cout << " fCentFrom:      " << fCentFrom << endl;
 289   cout << " fCentTo:        " << fCentTo << endl;
 290   cout << " fVtxZMin:       " << fVtxZMin << endl;
 291   cout << " fVtxZMax:       " << fVtxZMax << endl;
 292   cout << " fUseQualFlag:   " << fUseQualFlag << endl;
 293   cout << " fClusName:      \"" << fClusName << "\"" << endl;
 294   cout << " fDoNtuple:      " << fDoNtuple << endl;
 295   cout << " fDoAfterburner: " << fDoAfterburner << endl;
 296   cout << " fAsymMax:       " << fAsymMax << endl;
 297   cout << " fNminCells:     " << fNminCells << endl;
 298   cout << " fMinE:          " << fMinE << endl;
 299   cout << " fMinErat:       " << fMinErat << endl;
 300   cout << " fMinEcc:        " << fMinEcc << endl;
 301   cout << " fGeoName:       \"" << fGeoName << "\"" << endl;
 302   cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
 303   cout << " fIsoDist:       " << fIsoDist << endl;
 304   cout << " fTrClassNames:  \"" << fTrClassNames << "\"" << endl;
 305   cout << " fTrCuts:        " << fTrCuts << endl;
 306   cout << " fPrimTrCuts:    " << fPrimTrCuts << endl;
 307   cout << " fDoTrMatGeom:   " << fDoTrMatGeom << endl;
 308   cout << " fTrainMode:     " << fTrainMode << endl;
 309   cout << " fMarkCells:     " << fMarkCells << endl;
 310   cout << " fMinL0Time:     " << fMinL0Time << endl;
 311   cout << " fMaxL0Time:     " << fMaxL0Time << endl;
 312   cout << " fMcMode:        " << fMcMode << endl;
 313   cout << " fEmbedMode:     " << fEmbedMode << endl;
 314   cout << " fGeom:          " << fGeom << endl;
 315   cout << " fReco:          " << fReco << endl;
 316   cout << " fTrigName:      " << fTrigName << endl;
 317   cout << " fDoPSel:        " << fDoPSel << endl;
 318
 319   if (!fGeom)
 320     fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 321   else {
 322     if (fGeom->GetMatrixForSuperModule(0))
 323       fIsGeoMatsSet = kTRUE;
 324   }
 325   if (!fReco)
 326     fReco = new AliEMCALRecoUtils();
 327   fTrClassNamesArr = fTrClassNames.Tokenize(" ");
 328   fOutput = new TList();
 329   fOutput->SetOwner(1);
 330   fSelTracks = new TObjArray;
 331   fSelPrimTracks = new TObjArray;
 332
 333   if (fDoNtuple) {
 334     TFile *f = OpenFile(1);
 335     TDirectory::TContext context(f);
 336     if (f) {
 337       f->SetCompressionLevel(2);
 338       fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
 339       fNtuple->SetDirectory(f);
 340       if (fTrainMode) {
 341         fNtuple->SetAutoFlush(-2*1024*1024);
 342         fNtuple->SetAutoSave(0);
 343       } else {
 344         fNtuple->SetAutoFlush(-32*1024*1024);
 345         fNtuple->SetAutoSave(0);
 346       }
 347
 348       fHeader = new AliStaHeader;
 349       fNtuple->Branch("header", &fHeader, 16*1024, 99);
 350       fPrimVert = new AliStaVertex;
 351       fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
 352       fSpdVert = new AliStaVertex;
 353       fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
 354       fTpcVert = new AliStaVertex;
 355       fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
 356       if (TClass::GetClass("AliStaCluster"))
 357         TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
 358       fClusters = new TClonesArray("AliStaCluster");
 359       fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
 360       if (TClass::GetClass("AliStaTrigger"))
 361         TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
 362       fTriggers = new TClonesArray("AliStaTrigger");
 363       fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
 364       if (fMcMode||fEmbedMode) {
 365         if (TClass::GetClass("AliStaPart"))
 366           TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
 367         fMcParts = new TClonesArray("AliStaPart");
 368         fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
 369       }
 370     }
 371   }
 372
 373   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 374   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
 375   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
 376
 377   // histograms
 378   Bool_t th1 =   TH1::GetDefaultSumw2();
 379   TH1::SetDefaultSumw2(kTRUE);
 380   Bool_t th2 =   TH2::GetDefaultSumw2();
 381   TH2::SetDefaultSumw2(kTRUE);
 382   fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
 383   fHCuts->GetXaxis()->SetBinLabel(1,"All");
 384   fHCuts->GetXaxis()->SetBinLabel(2,"PS");
 385   fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 386   fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
 387   fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 388   fOutput->Add(fHCuts);
 389   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 390   fHVertexZ->SetXTitle("z [cm]");
 391   fOutput->Add(fHVertexZ);
 392   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 393   fHVertexZ2->SetXTitle("z [cm]");
 394   fOutput->Add(fHVertexZ2);
 395   fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
 396   fHCent->SetXTitle(fCentVar.Data());
 397   fOutput->Add(fHCent);
 398   fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
 399   fHCentQual->SetXTitle(fCentVar.Data());
 400   fOutput->Add(fHCentQual);
 401   fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 402   fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 403   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 404     const char *name = fTrClassNamesArr->At(i)->GetName();
 405     fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
 406     fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
 407   }
 408   fOutput->Add(fHTclsBeforeCuts);
 409   fOutput->Add(fHTclsAfterCuts);
 410
 411   // histograms for cells
 412   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 413   fHColuRow   = new TH2*[nsm];
 414   fHColuRowE  = new TH2*[nsm];
 415   fHCellMult  = new TH1*[nsm];
 416   for (Int_t i = 0; i<nsm; ++i) {
 417     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
 418     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 419     fHColuRow[i]->SetXTitle("col (i#eta)");
 420     fHColuRow[i]->SetYTitle("row (i#phi)");
 421     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
 422     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 423     fHColuRowE[i]->SetXTitle("col (i#eta)");
 424     fHColuRowE[i]->SetYTitle("row (i#phi)");
 425     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 426     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 427     fHCellMult[i]->SetXTitle("# of cells");
 428     fOutput->Add(fHColuRow[i]);
 429     fOutput->Add(fHColuRowE[i]);
 430     fOutput->Add(fHCellMult[i]);
 431   }
 432   fHCellE = new TH1F("hCellE","",250,0.,25.);
 433   fHCellE->SetXTitle("E_{cell} [GeV]");
 434   fOutput->Add(fHCellE);
 435   fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
 436   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 437   fOutput->Add(fHCellH);
 438   fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
 439   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 440   fOutput->Add(fHCellM);
 441   fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
 442   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 443   fOutput->Add(fHCellM2);
 444
 445   fHCellFreqNoCut   = new TH1*[nsm];
 446   fHCellFreqCut100M = new TH1*[nsm];
 447   fHCellFreqCut300M = new TH1*[nsm];
 448   fHCellFreqE       = new TH1*[nsm];
 449   for (Int_t i = 0; i<nsm; ++i){
 450     Double_t lbin = i*24*48-0.5;
 451     Double_t hbin = lbin+24*48;
 452     fHCellFreqNoCut[i]   = new TH1F(Form("hCellFreqNoCut_SM%d",i),
 453                                     Form("Frequency SM%d (no cut);id;#",i),  1152, lbin, hbin);
 454     fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
 455                                     Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
 456     fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
 457                                     Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
 458     fHCellFreqE[i]       = new TH1F(Form("hCellFreqE_SM%d",i),
 459                                     Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
 460     fOutput->Add(fHCellFreqNoCut[i]);
 461     fOutput->Add(fHCellFreqCut100M[i]);
 462     fOutput->Add(fHCellFreqCut300M[i]);
 463     fOutput->Add(fHCellFreqE[i]);
 464   }
 465   if (!fMarkCells.IsNull()) {
 466     fHCellCheckE = new TH1*[24*48*nsm];
 467     memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
 468     TObjArray *cells = fMarkCells.Tokenize(" ");
 469     Int_t n = cells->GetEntries();
 470     Int_t *tcs = new Int_t[n];
 471     for (Int_t i=0;i<n;++i) {
 472       TString name(cells->At(i)->GetName());
 473       tcs[i]=name.Atoi();
 474     }
 475     for (Int_t i = 0; i<n; ++i) {
 476       Int_t c=tcs[i];
 477       if (c<24*48*nsm) {
 478         fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
 479         fOutput->Add(fHCellCheckE[i]);
 480       }
 481     }
 482     delete cells;
 483     delete [] tcs;
 484   }
 485
 486   // histograms for clusters
 487   if (!fTrainMode) {
 488     fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 489     fHClustEccentricity->SetXTitle("#epsilon_{C}");
 490     fOutput->Add(fHClustEccentricity);
 491     fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
 492     fHClustEtaPhi->SetXTitle("#eta");
 493     fHClustEtaPhi->SetYTitle("#varphi");
 494     fOutput->Add(fHClustEtaPhi);
 495     fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 496     fHClustEnergyPt->SetXTitle("E [GeV]");
 497     fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 498     fOutput->Add(fHClustEnergyPt);
 499     fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 500     fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 501     fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 502     fOutput->Add(fHClustEnergySigma);
 503     fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 504     fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 505     fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 506     fOutput->Add(fHClustSigmaSigma);
 507     fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 508     fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 509     fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 510     fOutput->Add(fHClustNCellEnergyRatio);
 511     fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
 512     fHClustEnergyNCell->SetXTitle("E_{clus}");
 513     fHClustEnergyNCell->SetYTitle("N_{cells}");
 514     fOutput->Add(fHClustEnergyNCell);
 515   }
 516
 517   // histograms for primary tracks
 518   fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
 519   fOutput->Add(fHPrimTrackPt);
 520   fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
 521   fOutput->Add(fHPrimTrackEta);
 522   fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
 523   fOutput->Add(fHPrimTrackPhi);
 524
 525   // histograms for track matching
 526   if (fDoTrMatGeom) {
 527     fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
 528     fOutput->Add(fHMatchDr);
 529     fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
 530     fOutput->Add(fHMatchDz);
 531     fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
 532     fOutput->Add(fHMatchEp);
 533   }
 534
 535   // histograms for pion candidates
 536   if (!fTrainMode) {
 537     fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
 538     fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 539     fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 540     fOutput->Add(fHPionEtaPhi);
 541     fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 542     fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 543     fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 544     fOutput->Add(fHPionMggPt);
 545     fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 546     fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 547     fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 548     fOutput->Add(fHPionMggAsym);
 549     fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 550     fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 551     fHPionMggDgg->SetYTitle("opening angle [grad]");
 552     fOutput->Add(fHPionMggDgg);
 553     const Int_t nbins = 20;
 554     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};
 555     fPtRanges = new TAxis(nbins-1,xbins);
 556     for (Int_t i = 0; i<=nbins; ++i) {
 557       fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 558       fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 559       if (i==0)
 560         fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 561       else if (i==nbins)
 562         fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 563       else
 564         fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 565       fOutput->Add(fHPionInvMasses[i]);
 566     }
 567   }
 568
 569   TH1::SetDefaultSumw2(th1);
 570   TH2::SetDefaultSumw2(th2);
 571   PostData(1, fOutput);
 572 }
 573
 574 //________________________________________________________________________
 575 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 576 {
 577   // Called for each event.
 578
 579   if (!InputEvent())
 580     return;
 581
 582   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 583   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 584   UInt_t offtrigger = 0;
 585   if (fEsdEv) {
 586     am->LoadBranch("AliESDRun.");
 587     am->LoadBranch("AliESDHeader.");
 588     UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
 589     UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
 590     Bool_t desc1 = (mask1 >> 18) & 0x1;
 591     Bool_t desc2 = (mask2 >> 18) & 0x1;
 592     if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(18)=="EMCAL"
 593       AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
 594                     mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
 595                     mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
 596       return;
 597     }
 598     offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
 599   } else {
 600     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 601     if (!fAodEv) {
 602       AliFatal("Neither ESD nor AOD event found");
 603       return;
 604     }
 605     am->LoadBranch("header");
 606     offtrigger =  ((AliVAODHeader*)fAodEv->GetHeader())->GetOfflineTrigger();
 607   }
 608   if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
 609     AliWarning(Form("EMCAL not in fast only partition"));
 610     return;
 611   }
 612   if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
 613     AliWarning("Accessing geometry from OCDB, this is not very efficient!");
 614     AliCDBManager *cdb = AliCDBManager::Instance();
 615     if (!cdb->IsDefaultStorageSet())
 616       cdb->SetDefaultStorage("raw://");
 617     Int_t runno = InputEvent()->GetRunNumber();
 618     if (runno != cdb->GetRun())
 619       cdb->SetRun(runno);
 620     AliGeomManager::LoadGeometry();
 621   }
 622
 623   if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
 624     Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 625     for (Int_t i=0; i<nsm; ++i) {
 626       const TGeoHMatrix *geom = 0;
 627       if (fEsdEv)
 628         geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
 629       else {
 630         AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(fAodEv->GetHeader());
 631         if(!aodheader) AliFatal("Not a standard AOD");
 632         geom = aodheader->GetEMCALMatrix(i);
 633       }
 634       if (!geom)
 635         continue;
 636       geom->Print();
 637       fGeom->SetMisalMatrix(geom,i);
 638     }
 639     fIsGeoMatsSet = kTRUE;
 640   }
 641
 642   if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 643     if (fEsdEv) {
 644       const AliESDRun *erun = fEsdEv->GetESDRun();
 645       AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 646                                                erun->GetCurrentDip(),
 647                                                AliMagF::kConvLHC,
 648                                                kFALSE,
 649                                                erun->GetBeamEnergy(),
 650                                                erun->GetBeamType());
 651       TGeoGlobalMagField::Instance()->SetField(field);
 652     } else {
 653       Double_t pol = -1; //polarity
 654       Double_t be = -1;  //beam energy
 655       AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 656       Int_t runno = fAodEv->GetRunNumber();
 657       if (runno>=136851 && runno<138275) {
 658         pol = -1;
 659         be = 2760;
 660         btype = AliMagF::kBeamTypeAA;
 661       } else if (runno>=138275 && runno<=139517) {
 662         pol = +1;
 663         be = 2760;
 664         btype = AliMagF::kBeamTypeAA;
 665       } else {
 666         AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 667       }
 668       TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 669     }
 670   }
 671
 672   Int_t cut = 1;
 673   fHCuts->Fill(cut++);
 674
 675   TString trgclasses;
 676   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 677   if (h) {
 678     trgclasses = fEsdEv->GetFiredTriggerClasses();
 679   } else {
 680     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 681     if (h2)
 682       trgclasses = h2->GetFiredTriggerClasses();
 683   }
 684   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 685     const char *name = fTrClassNamesArr->At(i)->GetName();
 686     TRegexp regexp(name);
 687     if (trgclasses.Contains(regexp))
 688       fHTclsBeforeCuts->Fill(1+i);
 689   }
 690
 691   if (fDoPSel && offtrigger==0)
 692     return;
 693
 694   fHCuts->Fill(cut++);
 695
 696   const AliCentrality *centP = InputEvent()->GetCentrality();
 697   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 698   fHCent->Fill(cent);
 699   if (cent<fCentFrom||cent>fCentTo)
 700     return;
 701
 702   fHCuts->Fill(cut++);
 703
 704   if (fUseQualFlag) {
 705     if (centP->GetQuality()>0)
 706       return;
 707   }
 708
 709   fHCentQual->Fill(cent);
 710   fHCuts->Fill(cut++);
 711
 712   if (fEsdEv) {
 713     am->LoadBranch("PrimaryVertex.");
 714     am->LoadBranch("SPDVertex.");
 715     am->LoadBranch("TPCVertex.");
 716   } else {
 717     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 718     am->LoadBranch("vertices");
 719     if (!fAodEv) return;
 720   }
 721
 722   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 723   if (!vertex)
 724     return;
 725
 726   fHVertexZ->Fill(vertex->GetZ());
 727
 728   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 729     return;
 730
 731   fHCuts->Fill(cut++);
 732   fHVertexZ2->Fill(vertex->GetZ());
 733
 734   // count number of accepted events
 735   ++fNEvs;
 736
 737   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 738     const char *name = fTrClassNamesArr->At(i)->GetName();
 739     TRegexp regexp(name);
 740     if (trgclasses.Contains(regexp))
 741       fHTclsAfterCuts->Fill(1+i);
 742   }
 743
 744   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 745   fDigits      = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 746   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
 747   fEsdCells    = 0; // will be set if ESD input used
 748   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
 749                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 750   fAodCells    = 0; // will be set if AOD input used
 751
 752   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 753   Bool_t overwrite    = 0;
 754   Bool_t clusattached = 0;
 755   Bool_t recalibrated = 0;
 756   if (1 && !fClusName.IsNull()) {
 757     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 758     TObjArray *ts = am->GetTasks();
 759     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 760     if (cltask && cltask->GetClusters()) {
 761       fRecPoints   = cltask->GetClusters();
 762       fDigits      = cltask->GetDigits();
 763       clusattached = cltask->GetAttachClusters();
 764       overwrite    = cltask->GetOverwrite();
 765       if (cltask->GetCalibData()!=0)
 766         recalibrated = kTRUE;
 767     }
 768   }
 769   if (1 && !fClusName.IsNull()) {
 770     TList *l = 0;
 771     if (AODEvent())
 772       l = AODEvent()->GetList();
 773     else if (fAodEv)
 774       l = fAodEv->GetList();
 775     if (l) {
 776       fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 777     }
 778   }
 779
 780   if (fEsdEv) { // ESD input mode
 781     if (1 && (!fRecPoints||clusattached)) {
 782       if (!clusattached && !overwrite)
 783         am->LoadBranch("CaloClusters");
 784       TList *l = fEsdEv->GetList();
 785       if (clusattached) {
 786         fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 787       } else {
 788         fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 789       }
 790     }
 791     if (1) {
 792       if (!recalibrated)
 793         am->LoadBranch("EMCALCells.");
 794       fEsdCells = fEsdEv->GetEMCALCells();
 795     }
 796   } else if (fAodEv) { // AOD input mode
 797     if (1 && (!fAodClusters || clusattached)) {
 798       if (!clusattached)
 799         am->LoadBranch("caloClusters");
 800       TList *l = fAodEv->GetList();
 801       if (l) {
 802         fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 803       }
 804     }
 805     if (1) {
 806       if (!recalibrated)
 807         am->LoadBranch("emcalCells");
 808       fAodCells = fAodEv->GetEMCALCells();
 809     }
 810   } else {
 811     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 812   }
 813
 814   if (1) {
 815     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 816     AliDebug(2,Form("fDigits      set: %p", fDigits));
 817     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 818     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 819     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 820     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 821   }
 822
 823   if (fDoAfterburner)
 824     ClusterAfterburner();
 825
 826   CalcMcInfo();
 827   CalcCaloTriggers();
 828   CalcPrimTracks();
 829   CalcTracks();
 830   CalcClusterProps();
 831
 832   FillCellHists();
 833   if (!fTrainMode) {
 834     FillClusHists();
 835     FillPionHists();
 836     FillOtherHists();
 837   }
 838   FillMcHists();
 839   FillNtuple();
 840
 841   if (fTrainMode) {
 842     fSelTracks->Clear();
 843     fSelPrimTracks->Clear();
 844     if (fMcParts)
 845       fMcParts->Clear();
 846     if (fTriggers)
 847       fTriggers->Clear();
 848     if (fClusters)
 849       fClusters->Clear();
 850   }
 851
 852   PostData(1, fOutput);
 853 }
 854
 855 //________________________________________________________________________
 856 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 857 {
 858   // Terminate called at the end of analysis.
 859
 860   if (fNtuple) {
 861     TFile *f = OpenFile(1);
 862     TDirectory::TContext context(f);
 863     if (f)
 864       fNtuple->Write();
 865   }
 866
 867   AliInfo(Form("%s: Accepted %lld events          ", GetName(), fNEvs));
 868 }
 869
 870 //________________________________________________________________________
 871 void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
 872 {
 873   // Calculate triggers
 874
 875   if (fAodEv)
 876     return; // information not available in AOD
 877
 878   if (!fTriggers)
 879     return;
 880
 881   fTriggers->Clear();
 882
 883   if (fTrigName.Length()<=0)
 884     return;
 885
 886   TClonesArray *arr = dynamic_cast<TClonesArray*>(fEsdEv->FindListObject(fTrigName));
 887   if (!arr) {
 888     AliError(Form("Could not get array with name %s", fTrigName.Data()));
 889     return;
 890   }
 891
 892   Int_t nNumberOfCaloClusters = arr->GetEntries();
 893   for(Int_t j = 0, ntrigs = 0; j < nNumberOfCaloClusters; ++j) {
 894     AliVCluster *cl = dynamic_cast<AliVCluster*>(arr->At(j));
 895     if (!cl)
 896       continue;
 897     if (!cl->IsEMCAL())
 898       continue;
 899     if (cl->E()<1)
 900       continue;
 901     AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
 902     Float_t pos[3] = {0,0,0};
 903     cl->GetPosition(pos);
 904     TVector3 vpos(pos);
 905     trignew->fE       = cl->E();
 906     trignew->fEta     = vpos.Eta();
 907     trignew->fPhi     = vpos.Phi();
 908     Short_t  id    = -1;
 909     GetMaxCellEnergy(cl, id);
 910     trignew->fIdMax    = id;
 911   }
 912 }
 913
 914 //________________________________________________________________________
 915 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 916 {
 917   // Calculate cluster properties
 918
 919   if (!fClusters)
 920     return;
 921
 922   fClusters->Clear();
 923
 924   AliVCaloCells *cells = fEsdCells;
 925   if (!cells)
 926     cells = fAodCells;
 927   if (!cells)
 928     return;
 929
 930   TObjArray *clusters = fEsdClusters;
 931   if (!clusters)
 932     clusters = fAodClusters;
 933   if (!clusters)
 934     return;
 935
 936   const Int_t ncells = cells->GetNumberOfCells();
 937   const Int_t nclus  = clusters->GetEntries();
 938   const Int_t ntrks  = fSelTracks->GetEntries();
 939
 940   std::map<Short_t,Short_t> map;
 941   for (Short_t pos=0; pos<ncells; ++pos) {
 942     Short_t id = cells->GetCellNumber(pos);
 943     if (id>24*48*fGeom->GetNumberOfSuperModules()) {
 944       AliError(Form("Id of cell found to be too large: %d", id));
 945       continue;
 946     }
 947     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
 948     if (digit->GetId() != id) {
 949       AliError(Form("Ids should be equal: %d %d", id, digit->GetId()));
 950       return;
 951     }
 952     map[id]    = pos;
 953   }
 954
 955   TObjArray filtMcParts;
 956   if (fMcParts) {
 957     Int_t nmc = fMcParts->GetEntries();
 958     for (Int_t i=0; i<nmc; ++i) {
 959       AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
 960       if (pa->OnEmcal())
 961         filtMcParts.Add(pa);
 962     }
 963   }
 964
 965   Double_t vertex[3] = {0,0,0};
 966   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 967
 968   for(Int_t i=0, ncl=0; i<nclus; ++i) {
 969     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 970
 971     if (!clus)
 972       continue;
 973     if (!clus->IsEMCAL())
 974       continue;
 975     if (clus->E()<fMinE)
 976       continue;
 977
 978     Float_t clsPos[3] = {0,0,0};
 979     clus->GetPosition(clsPos);
 980     TVector3 clsVec(clsPos);
 981     TLorentzVector clusterVec;
 982     clus->GetMomentum(clusterVec,vertex);
 983     Double_t clsEta = clusterVec.Eta();
 984
 985     AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
 986     cl->fE        = clus->E();
 987     cl->fR        = clsVec.Perp();
 988     cl->fEta      = clsVec.Eta();
 989     cl->fPhi      = clsVec.Phi();
 990     cl->fN        = clus->GetNCells();
 991     cl->fN1       = GetNCells(clus,0.100);
 992     cl->fN3       = GetNCells(clus,0.300);
 993     Short_t id    = -1;
 994     Double_t emax = GetMaxCellEnergy(clus, id);
 995     cl->fIdMax    = id;
 996     cl->fSM       = fGeom->GetSuperModuleNumber(id);
 997     cl->fEmax     = emax;
 998     Short_t id2   = -1;
 999     cl->fE2max    = GetSecondMaxCellEnergy(clus,id2);
1000     cl->fTmax     = cells->GetCellTime(id);
1001     if (clus->GetDistanceToBadChannel()<10000)
1002       cl->fDbc    = clus->GetDistanceToBadChannel();
1003     if (!TMath::IsNaN(clus->GetDispersion()))
1004       cl->fDisp   = clus->GetDispersion();
1005     if (!TMath::IsNaN(clus->GetM20()))
1006       cl->fM20    = clus->GetM20();
1007     if (!TMath::IsNaN(clus->GetM02()))
1008       cl->fM02    = clus->GetM02();
1009     Double_t maxAxis = -1, minAxis = -1;
1010     GetSigma(clus,maxAxis,minAxis);
1011     clus->SetTOF(maxAxis);     // store sigma in TOF for later plotting
1012     cl->fSig      = maxAxis;
1013     Double_t sEtaEta = -1;
1014     Double_t sPhiPhi = -1;
1015     GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
1016     cl->fSigEtaEta = sEtaEta;
1017     cl->fSigPhiPhi = sPhiPhi;
1018     Double_t clusterEcc = -1;
1019     if (maxAxis > 0)
1020       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
1021     clus->SetChi2(clusterEcc); // store ecc in chi2 for later plotting
1022     cl->fEcc        = clusterEcc;
1023     cl->fTrIso      = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1024     cl->fTrIso1     = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1025     cl->fTrIso2     = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
1026     cl->fTrIsoD1    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
1027     cl->fTrIso1D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
1028     cl->fTrIso2D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
1029     cl->fTrIsoD3    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
1030     cl->fTrIso1D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
1031     cl->fTrIso2D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
1032     cl->fTrIsoD4    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2);
1033     cl->fTrIso1D4   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 1);
1034     cl->fTrIso2D4   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 2);
1035     cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
1036     cl->fCeCore     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1037     cl->fCeIso      = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1038     cl->fCeIso1     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
1039     cl->fCeIso3     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
1040     cl->fCeIso4     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.40);
1041     cl->fCeIso3x3   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 3);
1042     cl->fCeIso4x4   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
1043     cl->fCeIso5x5   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
1044     cl->fCeIso3x22  = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
1045     cl->fIsShared   = IsShared(clus);
1046     cl->fTrigId     = -1;
1047     cl->fTrigE      = 0;
1048
1049     if (fTriggers) {
1050       Int_t ntrig = fTriggers->GetEntries();
1051       for (Int_t j = 0; j<ntrig; ++j) {
1052         AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1053         if (!sta)
1054           continue;
1055         Short_t idmax = sta->fIdMax;
1056         Bool_t inc = IsIdPartOfCluster(clus, idmax);
1057         if (inc) {
1058           cl->fTrigId     = j;
1059           cl->fTrigE      = sta->fE;
1060           break;
1061         }
1062       }
1063     }
1064
1065     // track matching
1066     if (fDoTrMatGeom) {
1067       Double_t mind2 = 1e10;
1068       for(Int_t j = 0; j<ntrks; ++j) {
1069         AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1070         if (!track)
1071           continue;
1072         if (TMath::Abs(clsEta-track->Eta())>0.5)
1073           continue;
1074         TVector3 vec(clsPos);
1075         Float_t tmpEta=-1, tmpPhi=-1, tmpR2=1e10;
1076         Double_t dedx = -1;
1077         if (fEsdEv) {
1078           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1079           dedx = esdTrack->GetTPCsignal();
1080         }
1081         if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1082           continue;
1083         AliExternalTrackParam tParam;
1084         tParam.CopyFromVTrack(track);
1085         if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpEta, tmpPhi))
1086           continue;
1087         tmpR2 = tmpEta*tmpEta + tmpPhi*tmpPhi;
1088         Double_t d2 = tmpR2;
1089         if (mind2>d2) {
1090           mind2=d2;
1091           cl->fTrDz   = tmpEta;
1092           cl->fTrDr   = tmpPhi;
1093           cl->fTrEp   = clus->E()/track->P();
1094           cl->fTrDedx = dedx;
1095           cl->fIsTrackM = 1;
1096         }
1097       }
1098       if (cl->fIsTrackM) {
1099         fHMatchDr->Fill(cl->fTrDr);
1100         fHMatchDz->Fill(cl->fTrDz);
1101         fHMatchEp->Fill(cl->fTrEp);
1102       }
1103     }
1104
1105     //mc matching
1106     if (fMcParts) {
1107       Int_t nmc = filtMcParts.GetEntries();
1108       Double_t diffR2 = 1e9;
1109       AliStaPart *msta = 0;
1110       for (Int_t j=0; j<nmc; ++j) {
1111         AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1112         Double_t t1=clsVec.Eta()-pa->fVEta;
1113         Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1114         Double_t tmp = t1*t1+t2*t2;
1115         if (tmp<diffR2) {
1116           diffR2 = tmp;
1117           msta   = pa;
1118         }
1119       }
1120       if (diffR2<10 && msta!=0) {
1121         cl->fMcLabel = msta->fLab;
1122       }
1123     }
1124
1125     cl->fEmbE = 0;
1126     if (fDigits && fEmbedMode) {
1127       for(Int_t j=0; j<cl->fN; ++j) {
1128         Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1129         Short_t pos = -1;
1130         std::map<Short_t,Short_t>::iterator it = map.find(cid);
1131         if (it!=map.end())
1132           pos = it->second;
1133         if (pos<0)
1134           continue;
1135         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1136         if (!digit)
1137           continue;
1138         if (digit->GetType()<-1) {
1139           cl->fEmbE += digit->GetChi2();
1140         }
1141       }
1142     }
1143   }
1144 }
1145
1146 //________________________________________________________________________
1147 void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1148 {
1149   // Calculate track properties for primary tracks.
1150
1151   fSelPrimTracks->Clear();
1152
1153   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1154   TClonesArray *tracks = 0;
1155   Bool_t doConstrain = kFALSE;
1156   TString trkname(fPrimTracksName);
1157   if (fEsdEv) {
1158     if (trkname.IsNull()) {
1159       trkname = "Tracks";
1160       am->LoadBranch("Tracks");
1161       fSelPrimTracks->SetOwner(kTRUE);
1162       doConstrain = kTRUE;
1163     }
1164     TList *l = fEsdEv->GetList();
1165     tracks = dynamic_cast<TClonesArray*>(l->FindObject(trkname));
1166   } else {
1167     trkname = "tracks";
1168     am->LoadBranch("tracks");
1169     TList *l = fAodEv->GetList();
1170     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1171   }
1172
1173   if (!tracks) {
1174     AliError(Form("Pointer to tracks %s == 0", trkname.Data() ));
1175     return;
1176   }
1177
1178   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1179   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1180   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
1181
1182   if (fEsdEv) {
1183     am->LoadBranch("PrimaryVertex.");
1184     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1185     am->LoadBranch("SPDVertex.");
1186     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1187     am->LoadBranch("Tracks");
1188     const Int_t Ntracks = tracks->GetEntries();
1189     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1190       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1191       if (!track) {
1192         AliWarning(Form("Could not receive track %d\n", iTracks));
1193         continue;
1194       }
1195       if (fPrimTrCuts && !fPrimTrCuts->IsSelected(track))
1196         continue;
1197       Double_t eta = track->Eta();
1198       if (eta<-1||eta>1)
1199         continue;
1200       if (track->Phi()<phimin||track->Phi()>phimax)
1201         continue;
1202       if (!doConstrain) {
1203         fSelPrimTracks->Add(track);
1204         continue;
1205       }
1206       AliESDtrack copyt(*track);
1207       Double_t bfield[3];
1208       copyt.GetBxByBz(bfield);
1209       AliExternalTrackParam tParam;
1210       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1211       if (!relate)
1212         continue;
1213       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1214       Double_t p[3]      = { 0. };
1215       copyt.GetPxPyPz(p);
1216       Double_t pos[3]    = { 0. };
1217       copyt.GetXYZ(pos);
1218       Double_t covTr[21] = { 0. };
1219       copyt.GetCovarianceXYZPxPyPz(covTr);
1220       //      Double_t pid[10]   = { 0. };
1221       //      copyt.GetESDpid(pid);
1222       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1223                                             copyt.GetLabel(),
1224                                             p,
1225                                             kTRUE,
1226                                             pos,
1227                                             kFALSE,
1228                                             covTr,
1229                                             (Short_t)copyt.GetSign(),
1230                                             copyt.GetITSClusterMap(),
1231                                             //pid,
1232                                             0,/*fPrimaryVertex,*/
1233                                             kTRUE, // check if this is right
1234                                             vtx->UsesTrack(copyt.GetID()));
1235       aTrack->SetPIDForTracking(copyt.GetPIDForTracking());
1236       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1237       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1238       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1239       if(ndf>0)
1240         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1241       else
1242         aTrack->SetChi2perNDF(-1);
1243       aTrack->SetFlags(copyt.GetStatus());
1244       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
1245       fSelPrimTracks->Add(aTrack);
1246     }
1247   } else {
1248     Int_t ntracks = tracks->GetEntries();
1249     for (Int_t i=0; i<ntracks; ++i) {
1250       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1251       if (!track)
1252         continue;
1253       Double_t eta = track->Eta();
1254       if (eta<-1||eta>1)
1255         continue;
1256       if (track->Phi()<phimin||track->Phi()>phimax)
1257         continue;
1258       if(track->GetTPCNcls()<fMinNClusPerTr)
1259         continue;
1260       //todo: Learn how to set/filter AODs for prim/sec tracks
1261       fSelPrimTracks->Add(track);
1262     }
1263   }
1264 }
1265
1266 //________________________________________________________________________
1267 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
1268 {
1269   // Calculate track properties (including secondaries).
1270
1271   fSelTracks->Clear();
1272
1273   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1274   TClonesArray *tracks = 0;
1275   if (fEsdEv) {
1276     am->LoadBranch("Tracks");
1277     TList *l = fEsdEv->GetList();
1278     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1279   } else {
1280     am->LoadBranch("tracks");
1281     TList *l = fAodEv->GetList();
1282     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1283   }
1284
1285   if (!tracks)
1286     return;
1287
1288   if (fEsdEv) {
1289     const Int_t Ntracks = tracks->GetEntries();
1290     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1291       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1292       if (!track) {
1293         AliWarning(Form("Could not receive track %d\n", iTracks));
1294         continue;
1295       }
1296       if (fTrCuts && !fTrCuts->IsSelected(track))
1297         continue;
1298       Double_t eta = track->Eta();
1299       if (eta<-1||eta>1)
1300         continue;
1301       fSelTracks->Add(track);
1302     }
1303   } else {
1304     Int_t ntracks = tracks->GetEntries();
1305     for (Int_t i=0; i<ntracks; ++i) {
1306       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1307       if (!track)
1308         continue;
1309       Double_t eta = track->Eta();
1310       if (eta<-1||eta>1)
1311         continue;
1312       if(track->GetTPCNcls()<fMinNClusPerTr)
1313         continue;
1314
1315       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
1316         AliExternalTrackParam tParam;
1317         tParam.CopyFromVTrack(track);
1318         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 430, 0.139, 1, kTRUE)) {
1319           Double_t trkPos[3];
1320           tParam.GetXYZ(trkPos);
1321           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1322         }
1323       }
1324       fSelTracks->Add(track);
1325     }
1326   }
1327 }
1328
1329 //________________________________________________________________________
1330 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
1331 {
1332   // Run custer reconstruction afterburner.
1333
1334   AliVCaloCells *cells = fEsdCells;
1335   if (!cells)
1336     cells = fAodCells;
1337
1338   if (!cells)
1339     return;
1340
1341   Int_t ncells = cells->GetNumberOfCells();
1342   if (ncells<=0)
1343     return;
1344
1345   Double_t cellMeanE = 0, cellSigE = 0;
1346   for (Int_t i = 0; i<ncells; ++i) {
1347     Double_t cellE = cells->GetAmplitude(i);
1348     cellMeanE += cellE;
1349     cellSigE += cellE*cellE;
1350   }
1351   cellMeanE /= ncells;
1352   cellSigE /= ncells;
1353   cellSigE -= cellMeanE*cellMeanE;
1354   if (cellSigE<0)
1355     cellSigE = 0;
1356   cellSigE = TMath::Sqrt(cellSigE / ncells);
1357
1358   Double_t subE = cellMeanE - 7*cellSigE;
1359   if (subE<0)
1360     return;
1361
1362   for (Int_t i = 0; i<ncells; ++i) {
1363     Short_t id=-1;
1364     Int_t mclabel = -1;
1365     Double_t amp=0,time=0, efrac = 0;
1366     if (!cells->GetCell(i, id, amp, time,mclabel,efrac))
1367       continue;
1368     amp -= cellMeanE;
1369     if (amp<0.001)
1370       amp = 0;
1371     cells->SetCell(i, id, amp, time);
1372   }
1373
1374   TObjArray *clusters = fEsdClusters;
1375   if (!clusters)
1376     clusters = fAodClusters;
1377   if (!clusters)
1378     return;
1379
1380   Int_t nclus = clusters->GetEntries();
1381   for (Int_t i = 0; i<nclus; ++i) {
1382     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1383     if (!clus->IsEMCAL())
1384       continue;
1385     Int_t nc = clus->GetNCells();
1386     Double_t clusE = 0;
1387     UShort_t ids[100] = {0};
1388     Double_t fra[100] = {0};
1389     for (Int_t j = 0; j<nc; ++j) {
1390       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1391       Double_t cen = cells->GetCellAmplitude(id);
1392       clusE += cen;
1393       if (cen>0) {
1394         ids[nc] = id;
1395         ++nc;
1396       }
1397     }
1398     if (clusE<=0) {
1399       clusters->RemoveAt(i);
1400       continue;
1401     }
1402
1403     for (Int_t j = 0; j<nc; ++j) {
1404       Short_t id = ids[j];
1405       Double_t cen = cells->GetCellAmplitude(id);
1406       fra[j] = cen/clusE;
1407     }
1408     clus->SetE(clusE);
1409     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1410     if (aodclus) {
1411       aodclus->Clear("");
1412       aodclus->SetNCells(nc);
1413       aodclus->SetCellsAmplitudeFraction(fra);
1414       aodclus->SetCellsAbsId(ids);
1415     }
1416   }
1417   clusters->Compress();
1418 }
1419
1420 //________________________________________________________________________
1421 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1422 {
1423   // Fill histograms related to cell properties.
1424
1425   AliVCaloCells *cells = fEsdCells;
1426   if (!cells)
1427     cells = fAodCells;
1428
1429   if (!cells)
1430     return;
1431
1432   Int_t cellModCount[12] = {0};
1433   Double_t cellMaxE = 0;
1434   Double_t cellMeanE = 0;
1435   Int_t ncells = cells->GetNumberOfCells();
1436   if (ncells==0)
1437     return;
1438
1439   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1440
1441   for (Int_t i = 0; i<ncells; ++i) {
1442     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1443     Double_t cellE = cells->GetAmplitude(i);
1444     fHCellE->Fill(cellE);
1445     if (cellE>cellMaxE)
1446       cellMaxE = cellE;
1447     cellMeanE += cellE;
1448
1449     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1450     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1451     if (!ret) {
1452       AliError(Form("Could not get cell index for %d", absID));
1453       continue;
1454     }
1455     ++cellModCount[iSM];
1456     Int_t iPhi=-1, iEta=-1;
1457     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1458     fHColuRow[iSM]->Fill(iEta,iPhi,1);
1459     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1460     fHCellFreqNoCut[iSM]->Fill(absID);
1461     if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1462     if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
1463     if (fHCellCheckE && fHCellCheckE[absID])
1464       fHCellCheckE[absID]->Fill(cellE);
1465     fHCellFreqE[iSM]->Fill(absID, cellE);
1466   }
1467   fHCellH->Fill(cellMaxE);
1468   cellMeanE /= ncells;
1469   fHCellM->Fill(cellMeanE);
1470   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1471   for (Int_t i=0; i<nsm; ++i)
1472     fHCellMult[i]->Fill(cellModCount[i]);
1473 }
1474
1475 //________________________________________________________________________
1476 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1477 {
1478   // Fill histograms related to cluster properties.
1479
1480   TObjArray *clusters = fEsdClusters;
1481   if (!clusters)
1482     clusters = fAodClusters;
1483   if (!clusters)
1484     return;
1485
1486   Double_t vertex[3] = {0};
1487   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1488
1489   Int_t nclus = clusters->GetEntries();
1490   for(Int_t i = 0; i<nclus; ++i) {
1491     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1492     if (!clus)
1493       continue;
1494     if (!clus->IsEMCAL())
1495       continue;
1496     TLorentzVector clusterVec;
1497     clus->GetMomentum(clusterVec,vertex);
1498     Double_t maxAxis    = clus->GetTOF(); //sigma
1499     Double_t clusterEcc = clus->Chi2();   //eccentricity
1500     fHClustEccentricity->Fill(clusterEcc);
1501     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1502     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1503     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1504     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1505     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1506     fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
1507   }
1508 }
1509
1510 //________________________________________________________________________
1511 void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1512 {
1513   // Get Mc truth particle information.
1514
1515   if (!fMcParts)
1516     return;
1517
1518   fMcParts->Clear();
1519
1520   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1521   Double_t etamin = -0.7;
1522   Double_t etamax = +0.7;
1523   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1524   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1525
1526   if (fAodEv) {
1527     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1528     am->LoadBranch(AliAODMCParticle::StdBranchName());
1529     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1530     if (!tca)
1531       return;
1532
1533     Int_t nents = tca->GetEntries();
1534     for(int it=0; it<nents; ++it) {
1535       AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
1536       part->Print();
1537
1538       // pion or eta meson or direct photon
1539       if(part->GetPdgCode() == 111) {
1540       } else if(part->GetPdgCode() == 221) {
1541       } else if(part->GetPdgCode() == 22 ) {
1542       } else
1543         continue;
1544
1545       // primary particle
1546       Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1547       if(dR > 1.0)
1548         continue;
1549
1550       // kinematic cuts
1551       Double_t pt = part->Pt() ;
1552       if (pt<0.5)
1553         continue;
1554       Double_t eta = part->Eta();
1555       if (eta<etamin||eta>etamax)
1556         continue;
1557       Double_t phi  = part->Phi();
1558       if (phi<phimin||phi>phimax)
1559         continue;
1560
1561       ProcessDaughters(part, it, tca);
1562     }
1563     return;
1564   }
1565
1566   AliMCEvent *mcEvent = MCEvent();
1567   if (!mcEvent)
1568     return;
1569
1570   const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1571   if (!evtVtx)
1572     return;
1573
1574   mcEvent->PreReadAll();
1575
1576   Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1577   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1578     AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1579     if (!mcP)
1580       continue;
1581
1582     // pion or eta meson or direct photon
1583     if(mcP->PdgCode() == 111) {
1584     } else if(mcP->PdgCode() == 221) {
1585     } else if(mcP->PdgCode() == 22 ) {
1586     } else
1587       continue;
1588
1589     // primary particle
1590     Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1591                               (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1592     if(dR > 1.0)
1593       continue;
1594
1595     // kinematic cuts
1596     Double_t pt = mcP->Pt() ;
1597     if (pt<0.5)
1598       continue;
1599     Double_t eta = mcP->Eta();
1600     if (eta<etamin||eta>etamax)
1601       continue;
1602     Double_t phi  = mcP->Phi();
1603     if (phi<phimin||phi>phimax)
1604       continue;
1605
1606     ProcessDaughters(mcP, iTrack, mcEvent);
1607   }
1608 }
1609
1610 //________________________________________________________________________
1611 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1612 {
1613   // Fill ntuple.
1614
1615   if (!fNtuple)
1616     return;
1617
1618   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1619   if (fAodEv) {
1620     AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(fAodEv->GetHeader());
1621     if(!aodheader) AliFatal("Not a standard AOD");
1622
1623     fHeader->fRun            = fAodEv->GetRunNumber();
1624     fHeader->fOrbit          = aodheader->GetOrbitNumber();
1625     fHeader->fPeriod         = aodheader->GetPeriodNumber();
1626     fHeader->fBx             = aodheader->GetBunchCrossNumber();
1627     fHeader->fL0             = aodheader->GetL0TriggerInputs();
1628     fHeader->fL1             = aodheader->GetL1TriggerInputs();
1629     fHeader->fL2             = aodheader->GetL2TriggerInputs();
1630     fHeader->fTrClassMask    = aodheader->GetTriggerMask();
1631     fHeader->fTrCluster      = aodheader->GetTriggerCluster();
1632     fHeader->fOffTriggers    = aodheader->GetOfflineTrigger();
1633     fHeader->fFiredTriggers  = aodheader->GetFiredTriggerClasses();
1634   } else {
1635     fHeader->fRun            = fEsdEv->GetRunNumber();
1636     fHeader->fOrbit          = fEsdEv->GetHeader()->GetOrbitNumber();
1637     fHeader->fPeriod         = fEsdEv->GetESDRun()->GetPeriodNumber();
1638     fHeader->fBx             = fEsdEv->GetHeader()->GetBunchCrossNumber();
1639     fHeader->fL0             = fEsdEv->GetHeader()->GetL0TriggerInputs();
1640     fHeader->fL1             = fEsdEv->GetHeader()->GetL1TriggerInputs();
1641     fHeader->fL2             = fEsdEv->GetHeader()->GetL2TriggerInputs();
1642     fHeader->fTrClassMask    = fEsdEv->GetHeader()->GetTriggerMask();
1643     fHeader->fTrCluster      = fEsdEv->GetHeader()->GetTriggerCluster();
1644     fHeader->fOffTriggers    = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1645     fHeader->fFiredTriggers  = fEsdEv->GetFiredTriggerClasses();
1646     Float_t v0CorrR = 0;
1647     fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1648     const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1649     if (mult)
1650       fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1651     fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
1652     AliTriggerAnalysis trAn; /// Trigger Analysis
1653     Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1654     Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
1655     fHeader->fV0And        = v0A && v0B;
1656     fHeader->fIsHT         = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1657     am->LoadBranch("SPDPileupVertices");
1658     am->LoadBranch("TrkPileupVertices");
1659     fHeader->fIsPileup     = fEsdEv->IsPileupFromSPD(3,0.8);
1660     fHeader->fIsPileup2    = fEsdEv->IsPileupFromSPD(3,0.4);
1661     fHeader->fIsPileup4    = fEsdEv->IsPileupFromSPD(3,0.2);
1662     fHeader->fIsPileup8    = fEsdEv->IsPileupFromSPD(3,0.1);
1663     fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1664     fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
1665   }
1666
1667   AliCentrality *cent = InputEvent()->GetCentrality();
1668   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
1669   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
1670   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
1671   fHeader->fCqual     = cent->GetQuality();
1672
1673   AliEventplane *ep = InputEvent()->GetEventplane();
1674   if (ep) {
1675     if (ep->GetQVector())
1676       fHeader->fPsi     = ep->GetQVector()->Phi()/2. ;
1677     else
1678       fHeader->fPsi = -1;
1679     if (ep->GetQsub1()&&ep->GetQsub2())
1680       fHeader->fPsiRes  = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1681     else
1682       fHeader->fPsiRes = 0;
1683   }
1684
1685   Double_t val = 0;
1686   TString trgclasses(fHeader->fFiredTriggers);
1687   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1688     const char *name = fTrClassNamesArr->At(j)->GetName();
1689     TRegexp regexp(name);
1690     if (trgclasses.Contains(regexp))
1691       val += TMath::Power(2,j);
1692   }
1693   fHeader->fTcls = (UInt_t)val;
1694
1695   fHeader->fNSelTr     = fSelTracks->GetEntries();
1696   fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1697   fHeader->fNSelPrimTr1   = 0;
1698   fHeader->fNSelPrimTr2   = 0;
1699   for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1700     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1701     if(track->Pt()>1)
1702       ++fHeader->fNSelPrimTr1;
1703     if(track->Pt()>2)
1704       ++fHeader->fNSelPrimTr2;
1705   }
1706
1707   fHeader->fNCells   = 0;
1708   fHeader->fNCells0  = 0;
1709   fHeader->fNCells01 = 0;
1710   fHeader->fNCells03 = 0;
1711   fHeader->fNCells1  = 0;
1712   fHeader->fNCells2  = 0;
1713   fHeader->fNCells5  = 0;
1714   fHeader->fMaxCellE = 0;
1715
1716   AliVCaloCells *cells = fEsdCells;
1717   if (!cells)
1718     cells = fAodCells;
1719
1720   if (cells) {
1721     Int_t ncells = cells->GetNumberOfCells();
1722     for(Int_t j=0; j<ncells; ++j) {
1723       Double_t cellen = cells->GetAmplitude(j);
1724       if (cellen>0.045)
1725         ++fHeader->fNCells0;
1726       if (cellen>0.1)
1727         ++fHeader->fNCells01;
1728       if (cellen>0.3)
1729         ++fHeader->fNCells03;
1730       if (cellen>1)
1731         ++fHeader->fNCells1;
1732       if (cellen>2)
1733         ++fHeader->fNCells2;
1734       if (cellen>5)
1735         ++fHeader->fNCells5;
1736       if (cellen>fHeader->fMaxCellE)
1737         fHeader->fMaxCellE = cellen;
1738     }
1739     fHeader->fNCells = ncells;
1740   }
1741
1742   fHeader->fNClus      = 0;
1743   fHeader->fNClus1     = 0;
1744   fHeader->fNClus2     = 0;
1745   fHeader->fNClus5     = 0;
1746   fHeader->fMaxClusE   = 0;
1747
1748   TObjArray *clusters = fEsdClusters;
1749   if (!clusters)
1750     clusters = fAodClusters;
1751
1752   if (clusters) {
1753     Int_t nclus = clusters->GetEntries();
1754     for(Int_t j=0; j<nclus; ++j) {
1755       AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1756       if (!clus->IsEMCAL())
1757         continue;
1758       Double_t clusen = clus->E();
1759       if (clusen>1)
1760         ++fHeader->fNClus1;
1761       if (clusen>2)
1762         ++fHeader->fNClus2;
1763       if (clusen>5)
1764         ++fHeader->fNClus5;
1765       if (clusen>fHeader->fMaxClusE)
1766         fHeader->fMaxClusE = clusen;
1767     }
1768     fHeader->fNClus = nclus;
1769   }
1770
1771   fHeader->fMaxTrE     = 0;
1772   if (fTriggers) {
1773     Int_t ntrig = fTriggers->GetEntries();
1774     for (Int_t j = 0; j<ntrig; ++j) {
1775       AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1776       if (!sta)
1777         continue;
1778       if (sta->fE>fHeader->fMaxTrE)
1779         fHeader->fMaxTrE = sta->fE;
1780     }
1781   }
1782
1783   // count cells above 100 MeV on super modules
1784   fHeader->fNcSM0 = GetNCells(0, 0.100);
1785   fHeader->fNcSM1 = GetNCells(1, 0.100);
1786   fHeader->fNcSM2 = GetNCells(2, 0.100);
1787   fHeader->fNcSM3 = GetNCells(3, 0.100);
1788   fHeader->fNcSM4 = GetNCells(4, 0.100);
1789   fHeader->fNcSM5 = GetNCells(5, 0.100);
1790   fHeader->fNcSM6 = GetNCells(6, 0.100);
1791   fHeader->fNcSM7 = GetNCells(7, 0.100);
1792   fHeader->fNcSM8 = GetNCells(8, 0.100);
1793   fHeader->fNcSM9 = GetNCells(9, 0.100);
1794
1795   if (fAodEv) {
1796     am->LoadBranch("vertices");
1797     AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1798     FillVertex(fPrimVert, pv);
1799     AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1800     FillVertex(fSpdVert, sv);
1801   } else {
1802     am->LoadBranch("PrimaryVertex.");
1803     const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1804     FillVertex(fPrimVert, pv);
1805     am->LoadBranch("SPDVertex.");
1806     const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1807     FillVertex(fSpdVert, sv);
1808     am->LoadBranch("TPCVertex.");
1809     const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1810     FillVertex(fTpcVert, tv);
1811   }
1812
1813   fNtuple->Fill();
1814 }
1815
1816 //________________________________________________________________________
1817 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1818 {
1819   // Fill histograms related to pions.
1820
1821   Double_t vertex[3] = {0};
1822   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1823
1824   TObjArray *clusters = fEsdClusters;
1825   if (!clusters)
1826     clusters = fAodClusters;
1827
1828   if (clusters) {
1829     TLorentzVector clusterVec1;
1830     TLorentzVector clusterVec2;
1831     TLorentzVector pionVec;
1832
1833     Int_t nclus = clusters->GetEntries();
1834     for (Int_t i = 0; i<nclus; ++i) {
1835       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1836       if (!clus1)
1837         continue;
1838       if (!clus1->IsEMCAL())
1839         continue;
1840       if (clus1->E()<fMinE)
1841         continue;
1842       if (clus1->GetNCells()<fNminCells)
1843         continue;
1844       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1845         continue;
1846       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1847         continue;
1848       clus1->GetMomentum(clusterVec1,vertex);
1849       for (Int_t j = i+1; j<nclus; ++j) {
1850         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1851         if (!clus2)
1852           continue;
1853         if (!clus2->IsEMCAL())
1854           continue;
1855         if (clus2->E()<fMinE)
1856           continue;
1857         if (clus2->GetNCells()<fNminCells)
1858           continue;
1859         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1860           continue;
1861         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1862           continue;
1863         clus2->GetMomentum(clusterVec2,vertex);
1864         pionVec = clusterVec1 + clusterVec2;
1865         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1866         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1867         if (pionZgg < fAsymMax) {
1868           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1869           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1870           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1871           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1872           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1873           fHPionInvMasses[bin]->Fill(pionVec.M());
1874         }
1875       }
1876     }
1877   }
1878 }
1879
1880 //________________________________________________________________________
1881 void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1882 {
1883   // Fill additional MC information histograms.
1884
1885   if (!fMcParts)
1886     return;
1887
1888   // check if aod or esd mc mode and the fill histos
1889 }
1890
1891 //________________________________________________________________________
1892 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1893 {
1894   // Fill other histograms.
1895
1896   for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1897     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1898     if(!track)
1899       continue;
1900     fHPrimTrackPt->Fill(track->Pt());
1901     fHPrimTrackEta->Fill(track->Eta());
1902     fHPrimTrackPhi->Fill(track->Phi());
1903   }
1904 }
1905
1906 //________________________________________________________________________
1907 void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1908 {
1909   // Fill track histograms.
1910
1911   if (fSelPrimTracks) {
1912     for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1913       AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1914       if(!track)
1915         continue;
1916       fHPrimTrackPt->Fill(track->Pt());
1917       fHPrimTrackEta->Fill(track->Eta());
1918       fHPrimTrackPhi->Fill(track->Phi());
1919     }
1920   }
1921 }
1922
1923 //__________________________________________________________________________________________________
1924 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1925 {
1926   // Fill vertex from ESD vertex info.
1927
1928   v->fVx   = esdv->GetX();
1929   v->fVy   = esdv->GetY();
1930   v->fVz   = esdv->GetZ();
1931   v->fVc   = esdv->GetNContributors();
1932   v->fDisp = esdv->GetDispersion();
1933   v->fZres = esdv->GetZRes();
1934   v->fChi2 = esdv->GetChi2();
1935   v->fSt   = esdv->GetStatus();
1936   v->fIs3D = esdv->IsFromVertexer3D();
1937   v->fIsZ  = esdv->IsFromVertexerZ();
1938 }
1939
1940 //__________________________________________________________________________________________________
1941 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1942 {
1943   // Fill vertex from AOD vertex info.
1944
1945   v->fVx   = aodv->GetX();
1946   v->fVy   = aodv->GetY();
1947   v->fVz   = aodv->GetZ();
1948   v->fVc   = aodv->GetNContributors();
1949   v->fChi2 = aodv->GetChi2();
1950 }
1951
1952 //________________________________________________________________________
1953 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1954 {
1955   // Compute isolation based on cell content.
1956
1957   AliVCaloCells *cells = fEsdCells;
1958   if (!cells)
1959     cells = fAodCells;
1960   if (!cells)
1961     return 0;
1962
1963   Double_t cellIsolation = 0;
1964   Double_t rad2 = radius*radius;
1965   Int_t ncells = cells->GetNumberOfCells();
1966   for (Int_t i = 0; i<ncells; ++i) {
1967     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1968     Float_t eta=-1, phi=-1;
1969     fGeom->EtaPhiFromIndex(absID,eta,phi);
1970     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1971     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1972     if(dist>rad2)
1973       continue;
1974     Double_t cellE = cells->GetAmplitude(i);
1975     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1976     Double_t cellEt = cellE*sin(theta);
1977     cellIsolation += cellEt;
1978   }
1979   return cellIsolation;
1980 }
1981
1982 //________________________________________________________________________
1983 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1984 {
1985   // Compute isolation based on cell content, in a NxM rectangle.
1986
1987   AliVCaloCells *cells = fEsdCells;
1988   if (!cells)
1989     cells = fAodCells;
1990   if (!cells)
1991     return 0;
1992
1993   Double_t cellIsolation = 0;
1994   Int_t ncells = cells->GetNumberOfCells();
1995   for (Int_t i = 0; i<ncells; ++i) {
1996     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1997     Float_t eta=-1, phi=-1;
1998     fGeom->EtaPhiFromIndex(absID,eta,phi);
1999     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2000     Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2001     if(TMath::Abs(etadiff)/0.014>N)
2002       continue;
2003     if(TMath::Abs(phidiff)/0.014>M)
2004       continue;
2005     Double_t cellE = cells->GetAmplitude(i);
2006     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
2007     Double_t cellEt = cellE*sin(theta);
2008     cellIsolation += cellEt;
2009   }
2010   return cellIsolation;
2011 }
2012
2013 //________________________________________________________________________
2014 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
2015 {
2016   // Get maximum energy of attached cell.
2017
2018   Double_t ret = 0;
2019   Int_t ncells = cluster->GetNCells();
2020   if (fEsdCells) {
2021     for (Int_t i=0; i<ncells; i++) {
2022       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2023       ret += e;
2024       }
2025   } else {
2026     for (Int_t i=0; i<ncells; i++) {
2027       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2028       ret += e;
2029     }
2030   }
2031   return ret;
2032 }
2033
2034 //________________________________________________________________________
2035 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
2036 {
2037   // Get maximum energy of attached cell.
2038
2039   id = -1;
2040   Double_t maxe = 0;
2041   Int_t ncells = cluster->GetNCells();
2042   if (fEsdCells) {
2043     for (Int_t i=0; i<ncells; i++) {
2044       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2045       if (e>maxe) {
2046         maxe = e;
2047         id   = cluster->GetCellAbsId(i);
2048       }
2049     }
2050   } else {
2051     for (Int_t i=0; i<ncells; i++) {
2052       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2053       if (e>maxe)
2054         maxe = e;
2055         id   = cluster->GetCellAbsId(i);
2056     }
2057   }
2058   return maxe;
2059 }
2060
2061 //________________________________________________________________________
2062 Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
2063 {
2064   // Get second maximum cell.
2065
2066   AliVCaloCells *cells = fEsdCells;
2067   if (!cells)
2068     cells = fAodCells;
2069   if (!cells)
2070     return -1;
2071
2072   Double_t secondEmax=0, firstEmax=0;
2073   Double_t cellen;
2074   for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2075     Int_t absId = clus->GetCellAbsId(iCell);
2076     cellen = cells->GetCellAmplitude(absId);
2077     if(cellen > firstEmax)
2078       firstEmax = cellen;
2079   }
2080   for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2081     Int_t absId = clus->GetCellAbsId(iCell);
2082     cellen = cells->GetCellAmplitude(absId);
2083     if(cellen < firstEmax && cellen > secondEmax) {
2084       secondEmax = cellen;
2085       id = absId;
2086     }
2087   }
2088   return secondEmax;
2089 }
2090
2091 //________________________________________________________________________
2092 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
2093 {
2094   // Calculate the (E) weighted variance along the longer (eigen) axis.
2095
2096   sigmaMax = 0;          // cluster variance along its longer axis
2097   sigmaMin = 0;          // cluster variance along its shorter axis
2098   Double_t Ec  = c->E(); // cluster energy
2099   if(Ec<=0)
2100     return;
2101   Double_t Xc  = 0 ;     // cluster first moment along X
2102   Double_t Yc  = 0 ;     // cluster first moment along Y
2103   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
2104   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
2105   Double_t Syy = 0 ;     // cluster covariance^2
2106
2107   AliVCaloCells *cells = fEsdCells;
2108   if (!cells)
2109     cells = fAodCells;
2110
2111   if (!cells)
2112     return;
2113
2114   Int_t ncells = c->GetNCells();
2115   if (ncells==1)
2116     return;
2117
2118   for(Int_t j=0; j<ncells; ++j) {
2119     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2120     Double_t cellen = cells->GetCellAmplitude(id);
2121     TVector3 pos;
2122     fGeom->GetGlobal(id,pos);
2123     Xc  += cellen*pos.X();
2124     Yc  += cellen*pos.Y();
2125     Sxx += cellen*pos.X()*pos.X();
2126     Syy += cellen*pos.Y()*pos.Y();
2127     Sxy += cellen*pos.X()*pos.Y();
2128   }
2129   Xc  /= Ec;
2130   Yc  /= Ec;
2131   Sxx /= Ec;
2132   Syy /= Ec;
2133   Sxy /= Ec;
2134   Sxx -= Xc*Xc;
2135   Syy -= Yc*Yc;
2136   Sxy -= Xc*Yc;
2137   Sxx = TMath::Abs(Sxx);
2138   Syy = TMath::Abs(Syy);
2139   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2140   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2141   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2142   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
2143 }
2144
2145 //________________________________________________________________________
2146 void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2147 {
2148   // Calculate the (E) weighted variance along the pseudorapidity.
2149
2150   sEtaEta = 0;
2151   sPhiPhi = 0;
2152
2153   Double_t Ec  = c->E(); // cluster energy
2154   if(Ec<=0)
2155     return;
2156
2157   const Int_t ncells = c->GetNCells();
2158
2159   Double_t EtaC    = 0;  // cluster first moment along eta
2160   Double_t PhiC    = 0;  // cluster first moment along phi
2161   Double_t Setaeta = 0;  // cluster second central moment along eta
2162   Double_t Sphiphi = 0;  // cluster second central moment along phi
2163   Double_t w[ncells];    // weight max(0,4.5*log(E_i/Ec))
2164   Double_t sumw = 0;
2165   Int_t id[ncells];
2166
2167   AliVCaloCells *cells = fEsdCells;
2168   if (!cells)
2169     cells = fAodCells;
2170
2171   if (!cells)
2172     return;
2173
2174   if (ncells==1)
2175     return;
2176
2177   for(Int_t j=0; j<ncells; ++j) {
2178     id[j] = TMath::Abs(c->GetCellAbsId(j));
2179     Double_t cellen = cells->GetCellAmplitude(id[j]);
2180     w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2181     TVector3 pos;
2182     fGeom->GetGlobal(id[j],pos);
2183     EtaC += w[j]*pos.Eta();
2184     PhiC += w[j]*pos.Phi();
2185     sumw += w[j];
2186   }
2187   EtaC /= sumw;
2188   PhiC /= sumw;
2189
2190   for(Int_t j=0; j<ncells; ++j) {
2191     TVector3 pos;
2192     fGeom->GetGlobal(id[j],pos);
2193     Setaeta =  w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2194     Sphiphi =  w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2195   }
2196   Setaeta /= sumw;
2197   sEtaEta = TMath::Sqrt(Setaeta);
2198   Sphiphi /= sumw;
2199   sPhiPhi = TMath::Sqrt(Sphiphi);
2200 }
2201
2202 //________________________________________________________________________
2203 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
2204 {
2205   // Calculate number of attached cells above emin.
2206
2207   AliVCaloCells *cells = fEsdCells;
2208   if (!cells)
2209     cells = fAodCells;
2210   if (!cells)
2211     return 0;
2212
2213   Int_t n = 0;
2214   Int_t ncells = c->GetNCells();
2215   for(Int_t j=0; j<ncells; ++j) {
2216     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2217     Double_t cellen = cells->GetCellAmplitude(id);
2218     if (cellen>=emin)
2219       ++n;
2220   }
2221   return n;
2222 }
2223
2224 //________________________________________________________________________
2225 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2226 {
2227   // Calculate number of cells per SM above emin.
2228
2229   AliVCaloCells *cells = fEsdCells;
2230   if (!cells)
2231     cells = fAodCells;
2232   if (!cells)
2233     return 0;
2234
2235   Int_t n = 0;
2236   Int_t ncells = cells->GetNumberOfCells();
2237   for(Int_t j=0; j<ncells; ++j) {
2238     Int_t id = TMath::Abs(cells->GetCellNumber(j));
2239     Double_t cellen = cells->GetCellAmplitude(id);
2240     if (cellen<emin)
2241       continue;
2242     Int_t fsm = fGeom->GetSuperModuleNumber(id);
2243     if (fsm != sm)
2244       continue;
2245     ++n;
2246   }
2247   return n;
2248 }
2249
2250 //________________________________________________________________________
2251 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
2252 {
2253   // Compute isolation based on tracks.
2254
2255   Double_t trkIsolation = 0;
2256   Double_t rad2 = radius*radius;
2257   Int_t ntrks = fSelPrimTracks->GetEntries();
2258   for(Int_t j = 0; j<ntrks; ++j) {
2259     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
2260     if (!track)
2261       continue;
2262     if (track->Pt()<pt)
2263       continue;
2264     Float_t eta = track->Eta();
2265     Float_t phi = track->Phi();
2266     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2267     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2268     if(dist>rad2)
2269       continue;
2270     trkIsolation += track->Pt();
2271   }
2272   return trkIsolation;
2273 }
2274
2275 //________________________________________________________________________
2276 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2277 {
2278   // Compute isolation based on tracks.
2279
2280   Double_t trkIsolation = 0;
2281   Int_t ntrks = fSelPrimTracks->GetEntries();
2282   for(Int_t j = 0; j<ntrks; ++j) {
2283     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
2284     if (!track)
2285       continue;
2286     if (track->Pt()<pt)
2287       continue;
2288     Float_t eta = track->Eta();
2289     Float_t phi = track->Phi();
2290     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2291     Double_t etadiff = (eta-cEta);
2292     if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2293       continue;
2294     trkIsolation += track->Pt();
2295   }
2296   return trkIsolation;
2297 }
2298
2299 //________________________________________________________________________
2300 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2301 {
2302   // Returns if cluster shared across super modules.
2303
2304   if (!c)
2305     return 0;
2306
2307   Int_t n = -1;
2308   Int_t ncells = c->GetNCells();
2309   for(Int_t j=0; j<ncells; ++j) {
2310     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2311     Int_t got = id / (24*48);
2312     if (n==-1) {
2313       n = got;
2314       continue;
2315     }
2316     if (got!=n)
2317       return 1;
2318   }
2319   return 0;
2320 }
2321
2322 //________________________________________________________________________
2323 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2324 {
2325   // Returns if id is part of cluster.
2326
2327   AliVCaloCells *cells = fEsdCells;
2328   if (!cells)
2329     cells = fAodCells;
2330   if (!cells)
2331     return 0;
2332
2333   Int_t ncells = c->GetNCells();
2334   for(Int_t j=0; j<ncells; ++j) {
2335     Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2336     if (cid == id)
2337       return 1;
2338   }
2339   return 0;
2340 }
2341
2342 //________________________________________________________________________
2343 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
2344 {
2345   // Print recursively daughter information.
2346
2347   if (!p || !arr)
2348     return;
2349
2350   const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2351   if (!amc)
2352     return;
2353   for (Int_t i=0; i<level; ++i) printf(" ");
2354   amc->Print();
2355
2356   Int_t n = amc->GetNDaughters();
2357   for (Int_t i=0; i<n; ++i) {
2358     Int_t d = amc->GetDaughter(i);
2359     const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2360     PrintDaughters(dmc,arr,level+1);
2361   }
2362 }
2363
2364 //________________________________________________________________________
2365 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2366 {
2367   // Print recursively daughter information.
2368
2369   if (!p || !arr)
2370     return;
2371
2372   for (Int_t i=0; i<level; ++i) printf(" ");
2373   Int_t d1 = p->GetFirstDaughter();
2374   Int_t d2 = p->GetLastDaughter();
2375   printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2376          p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2377   if (d1<0)
2378     return;
2379   if (d2<0)
2380     d2=d1;
2381   for (Int_t i=d1;i<=d2;++i) {
2382     const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
2383     PrintDaughters(dmc,arr,level+1);
2384   }
2385 }
2386
2387 //________________________________________________________________________
2388 void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2389 {
2390   // Print track ref array.
2391
2392   if (!p)
2393     return;
2394
2395   Int_t n = p->GetNumberOfTrackReferences();
2396   for (Int_t i=0; i<n; ++i) {
2397     AliTrackReference *ref = p->GetTrackReference(i);
2398     if (!ref)
2399       continue;
2400     ref->SetUserId(ref->DetectorId());
2401     ref->Print();
2402   }
2403 }
2404
2405 //________________________________________________________________________
2406 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2407 {
2408   // Process and create daughters.
2409
2410   if (!p || !arr)
2411     return;
2412
2413   AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2414   if (!amc)
2415     return;
2416
2417   //amc->Print();
2418
2419   Int_t nparts = arr->GetEntries();
2420   Int_t nents  = fMcParts->GetEntries();
2421
2422   AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2423   newp->fPt  = amc->Pt();
2424   newp->fEta = amc->Eta();
2425   newp->fPhi = amc->Phi();
2426   if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2427     TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2428     newp->fVR = vec.Perp();
2429     newp->fVEta = vec.Eta();
2430     newp->fVPhi = vec.Phi();
2431   }
2432   newp->fPid  = amc->PdgCode();
2433   newp->fLab  = nents;
2434   Int_t moi = amc->GetMother();
2435   if (moi>=0&&moi<nparts) {
2436     const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2437     moi = mmc->GetUniqueID();
2438   }
2439   newp->fMo = moi;
2440   p->SetUniqueID(nents);
2441
2442   // TODO: Determine which detector was hit
2443   //newp->fDet = ???
2444
2445   Int_t n = amc->GetNDaughters();
2446   for (Int_t i=0; i<n; ++i) {
2447     Int_t d = amc->GetDaughter(i);
2448     if (d<=index || d>=nparts)
2449       continue;
2450     AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2451     ProcessDaughters(dmc,d,arr);
2452   }
2453 }
2454
2455 //________________________________________________________________________
2456 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2457 {
2458   // Process and create daughters.
2459
2460   if (!p || !arr)
2461     return;
2462
2463   Int_t d1 = p->GetFirstDaughter();
2464   Int_t d2 = p->GetLastDaughter();
2465   if (0) {
2466     printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2467            index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
2468     PrintTrackRefs(p);
2469   }
2470   Int_t nents  = fMcParts->GetEntries();
2471
2472   AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2473   newp->fPt  = p->Pt();
2474   newp->fEta = p->Eta();
2475   newp->fPhi = p->Phi();
2476   if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2477     TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2478     newp->fVR = vec.Perp();
2479     newp->fVEta = vec.Eta();
2480     newp->fVPhi = vec.Phi();
2481   }
2482   newp->fPid  = p->PdgCode();
2483   newp->fLab  = nents;
2484   Int_t moi = p->GetMother();
2485   if (moi>=0) {
2486     const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2487     moi = mmc->GetUniqueID();
2488   }
2489   newp->fMo = moi;
2490   p->SetUniqueID(nents);
2491
2492   Int_t nref = p->GetNumberOfTrackReferences();
2493   if (nref>0) {
2494     AliTrackReference *ref = p->GetTrackReference(nref-1);
2495     if (ref) {
2496       newp->fDet = ref->DetectorId();
2497     }
2498   }
2499
2500   if (d1<0)
2501     return;
2502   if (d2<0)
2503     d2=d1;
2504   for (Int_t i=d1;i<=d2;++i) {
2505     AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2506     if (dmc->P()<0.01)
2507       continue;
2508     ProcessDaughters(dmc,i,arr);
2509   }
2510 }
2511
2512 //__________________________________________________________________________________________________
2513 void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2514 {
2515   // Calculate momentum.
2516
2517   TVector3 pos;
2518   pos.SetPtEtaPhi(fR,fEta,fPhi);
2519
2520   if(vertex){ //calculate direction relative to vertex
2521     pos -= vertex;
2522   }
2523
2524   Double_t r = pos.Mag();
2525   p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2526 }
2527
2528 //__________________________________________________________________________________________________
2529 void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2530 {
2531   // Calculate momentum.
2532
2533   Double_t v[3] = {0,0,0};
2534   if (vertex) {
2535     v[0] = vertex->fVx;
2536     v[1] = vertex->fVy;
2537     v[2] = vertex->fVz;
2538   }
2539   GetMom(p, v);
2540 }