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 =  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         geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
 631       if (!geom)
 632         continue;
 633       geom->Print();
 634       fGeom->SetMisalMatrix(geom,i);
 635     }
 636     fIsGeoMatsSet = kTRUE;
 637   }
 638
 639   if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 640     if (fEsdEv) {
 641       const AliESDRun *erun = fEsdEv->GetESDRun();
 642       AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 643                                                erun->GetCurrentDip(),
 644                                                AliMagF::kConvLHC,
 645                                                kFALSE,
 646                                                erun->GetBeamEnergy(),
 647                                                erun->GetBeamType());
 648       TGeoGlobalMagField::Instance()->SetField(field);
 649     } else {
 650       Double_t pol = -1; //polarity
 651       Double_t be = -1;  //beam energy
 652       AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 653       Int_t runno = fAodEv->GetRunNumber();
 654       if (runno>=136851 && runno<138275) {
 655         pol = -1;
 656         be = 2760;
 657         btype = AliMagF::kBeamTypeAA;
 658       } else if (runno>=138275 && runno<=139517) {
 659         pol = +1;
 660         be = 2760;
 661         btype = AliMagF::kBeamTypeAA;
 662       } else {
 663         AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 664       }
 665       TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 666     }
 667   }
 668
 669   Int_t cut = 1;
 670   fHCuts->Fill(cut++);
 671
 672   TString trgclasses;
 673   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 674   if (h) {
 675     trgclasses = fEsdEv->GetFiredTriggerClasses();
 676   } else {
 677     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 678     if (h2)
 679       trgclasses = h2->GetFiredTriggerClasses();
 680   }
 681   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 682     const char *name = fTrClassNamesArr->At(i)->GetName();
 683     TRegexp regexp(name);
 684     if (trgclasses.Contains(regexp))
 685       fHTclsBeforeCuts->Fill(1+i);
 686   }
 687
 688   if (fDoPSel && offtrigger==0)
 689     return;
 690
 691   fHCuts->Fill(cut++);
 692
 693   const AliCentrality *centP = InputEvent()->GetCentrality();
 694   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 695   fHCent->Fill(cent);
 696   if (cent<fCentFrom||cent>fCentTo)
 697     return;
 698
 699   fHCuts->Fill(cut++);
 700
 701   if (fUseQualFlag) {
 702     if (centP->GetQuality()>0)
 703       return;
 704   }
 705
 706   fHCentQual->Fill(cent);
 707   fHCuts->Fill(cut++);
 708
 709   if (fEsdEv) {
 710     am->LoadBranch("PrimaryVertex.");
 711     am->LoadBranch("SPDVertex.");
 712     am->LoadBranch("TPCVertex.");
 713   } else {
 714     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 715     am->LoadBranch("vertices");
 716     if (!fAodEv) return;
 717   }
 718
 719   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 720   if (!vertex)
 721     return;
 722
 723   fHVertexZ->Fill(vertex->GetZ());
 724
 725   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 726     return;
 727
 728   fHCuts->Fill(cut++);
 729   fHVertexZ2->Fill(vertex->GetZ());
 730
 731   // count number of accepted events
 732   ++fNEvs;
 733
 734   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 735     const char *name = fTrClassNamesArr->At(i)->GetName();
 736     TRegexp regexp(name);
 737     if (trgclasses.Contains(regexp))
 738       fHTclsAfterCuts->Fill(1+i);
 739   }
 740
 741   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 742   fDigits      = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 743   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
 744   fEsdCells    = 0; // will be set if ESD input used
 745   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
 746                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 747   fAodCells    = 0; // will be set if AOD input used
 748
 749   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 750   Bool_t overwrite    = 0;
 751   Bool_t clusattached = 0;
 752   Bool_t recalibrated = 0;
 753   if (1 && !fClusName.IsNull()) {
 754     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 755     TObjArray *ts = am->GetTasks();
 756     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 757     if (cltask && cltask->GetClusters()) {
 758       fRecPoints   = cltask->GetClusters();
 759       fDigits      = cltask->GetDigits();
 760       clusattached = cltask->GetAttachClusters();
 761       overwrite    = cltask->GetOverwrite();
 762       if (cltask->GetCalibData()!=0)
 763         recalibrated = kTRUE;
 764     }
 765   }
 766   if (1 && !fClusName.IsNull()) {
 767     TList *l = 0;
 768     if (AODEvent())
 769       l = AODEvent()->GetList();
 770     else if (fAodEv)
 771       l = fAodEv->GetList();
 772     if (l) {
 773       fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 774     }
 775   }
 776
 777   if (fEsdEv) { // ESD input mode
 778     if (1 && (!fRecPoints||clusattached)) {
 779       if (!clusattached && !overwrite)
 780         am->LoadBranch("CaloClusters");
 781       TList *l = fEsdEv->GetList();
 782       if (clusattached) {
 783         fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 784       } else {
 785         fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 786       }
 787     }
 788     if (1) {
 789       if (!recalibrated)
 790         am->LoadBranch("EMCALCells.");
 791       fEsdCells = fEsdEv->GetEMCALCells();
 792     }
 793   } else if (fAodEv) { // AOD input mode
 794     if (1 && (!fAodClusters || clusattached)) {
 795       if (!clusattached)
 796         am->LoadBranch("caloClusters");
 797       TList *l = fAodEv->GetList();
 798       if (l) {
 799         fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 800       }
 801     }
 802     if (1) {
 803       if (!recalibrated)
 804         am->LoadBranch("emcalCells");
 805       fAodCells = fAodEv->GetEMCALCells();
 806     }
 807   } else {
 808     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 809   }
 810
 811   if (1) {
 812     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 813     AliDebug(2,Form("fDigits      set: %p", fDigits));
 814     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 815     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 816     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 817     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 818   }
 819
 820   if (fDoAfterburner)
 821     ClusterAfterburner();
 822
 823   CalcMcInfo();
 824   CalcCaloTriggers();
 825   CalcPrimTracks();
 826   CalcTracks();
 827   CalcClusterProps();
 828
 829   FillCellHists();
 830   if (!fTrainMode) {
 831     FillClusHists();
 832     FillPionHists();
 833     FillOtherHists();
 834   }
 835   FillMcHists();
 836   FillNtuple();
 837
 838   if (fTrainMode) {
 839     fSelTracks->Clear();
 840     fSelPrimTracks->Clear();
 841     if (fMcParts)
 842       fMcParts->Clear();
 843     if (fTriggers)
 844       fTriggers->Clear();
 845     if (fClusters)
 846       fClusters->Clear();
 847   }
 848
 849   PostData(1, fOutput);
 850 }
 851
 852 //________________________________________________________________________
 853 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 854 {
 855   // Terminate called at the end of analysis.
 856
 857   if (fNtuple) {
 858     TFile *f = OpenFile(1);
 859     TDirectory::TContext context(f);
 860     if (f)
 861       fNtuple->Write();
 862   }
 863
 864   AliInfo(Form("%s: Accepted %lld events          ", GetName(), fNEvs));
 865 }
 866
 867 //________________________________________________________________________
 868 void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
 869 {
 870   // Calculate triggers
 871
 872   if (fAodEv)
 873     return; // information not available in AOD
 874
 875   if (!fTriggers)
 876     return;
 877
 878   fTriggers->Clear();
 879
 880   if (fTrigName.Length()<=0)
 881     return;
 882
 883   TClonesArray *arr = dynamic_cast<TClonesArray*>(fEsdEv->FindListObject(fTrigName));
 884   if (!arr) {
 885     AliError(Form("Could not get array with name %s", fTrigName.Data()));
 886     return;
 887   }
 888
 889   Int_t nNumberOfCaloClusters = arr->GetEntries();
 890   for(Int_t j = 0, ntrigs = 0; j < nNumberOfCaloClusters; ++j) {
 891     AliVCluster *cl = dynamic_cast<AliVCluster*>(arr->At(j));
 892     if (!cl)
 893       continue;
 894     if (!cl->IsEMCAL())
 895       continue;
 896     if (cl->E()<1)
 897       continue;
 898     AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
 899     Float_t pos[3] = {0,0,0};
 900     cl->GetPosition(pos);
 901     TVector3 vpos(pos);
 902     trignew->fE       = cl->E();
 903     trignew->fEta     = vpos.Eta();
 904     trignew->fPhi     = vpos.Phi();
 905     Short_t  id    = -1;
 906     GetMaxCellEnergy(cl, id);
 907     trignew->fIdMax    = id;
 908   }
 909 }
 910
 911 //________________________________________________________________________
 912 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 913 {
 914   // Calculate cluster properties
 915
 916   if (!fClusters)
 917     return;
 918
 919   fClusters->Clear();
 920
 921   AliVCaloCells *cells = fEsdCells;
 922   if (!cells)
 923     cells = fAodCells;
 924   if (!cells)
 925     return;
 926
 927   TObjArray *clusters = fEsdClusters;
 928   if (!clusters)
 929     clusters = fAodClusters;
 930   if (!clusters)
 931     return;
 932
 933   const Int_t ncells = cells->GetNumberOfCells();
 934   const Int_t nclus  = clusters->GetEntries();
 935   const Int_t ntrks  = fSelTracks->GetEntries();
 936
 937   std::map<Short_t,Short_t> map;
 938   for (Short_t pos=0; pos<ncells; ++pos) {
 939     Short_t id = cells->GetCellNumber(pos);
 940     if (id>24*48*fGeom->GetNumberOfSuperModules()) {
 941       AliError(Form("Id of cell found to be too large: %d", id));
 942       continue;
 943     }
 944     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
 945     if (digit->GetId() != id) {
 946       AliError(Form("Ids should be equal: %d %d", id, digit->GetId()));
 947       return;
 948     }
 949     map[id]    = pos;
 950   }
 951
 952   TObjArray filtMcParts;
 953   if (fMcParts) {
 954     Int_t nmc = fMcParts->GetEntries();
 955     for (Int_t i=0; i<nmc; ++i) {
 956       AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
 957       if (pa->OnEmcal())
 958         filtMcParts.Add(pa);
 959     }
 960   }
 961
 962   Double_t vertex[3] = {0,0,0};
 963   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 964
 965   for(Int_t i=0, ncl=0; i<nclus; ++i) {
 966     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 967
 968     if (!clus)
 969       continue;
 970     if (!clus->IsEMCAL())
 971       continue;
 972     if (clus->E()<fMinE)
 973       continue;
 974
 975     Float_t clsPos[3] = {0,0,0};
 976     clus->GetPosition(clsPos);
 977     TVector3 clsVec(clsPos);
 978     TLorentzVector clusterVec;
 979     clus->GetMomentum(clusterVec,vertex);
 980     Double_t clsEta = clusterVec.Eta();
 981
 982     AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
 983     cl->fE        = clus->E();
 984     cl->fR        = clsVec.Perp();
 985     cl->fEta      = clsVec.Eta();
 986     cl->fPhi      = clsVec.Phi();
 987     cl->fN        = clus->GetNCells();
 988     cl->fN1       = GetNCells(clus,0.100);
 989     cl->fN3       = GetNCells(clus,0.300);
 990     Short_t id    = -1;
 991     Double_t emax = GetMaxCellEnergy(clus, id);
 992     cl->fIdMax    = id;
 993     cl->fSM       = fGeom->GetSuperModuleNumber(id);
 994     cl->fEmax     = emax;
 995     Short_t id2   = -1;
 996     cl->fE2max    = GetSecondMaxCellEnergy(clus,id2);
 997     cl->fTmax     = cells->GetCellTime(id);
 998     if (clus->GetDistanceToBadChannel()<10000)
 999       cl->fDbc    = clus->GetDistanceToBadChannel();
1000     if (!TMath::IsNaN(clus->GetDispersion()))
1001       cl->fDisp   = clus->GetDispersion();
1002     if (!TMath::IsNaN(clus->GetM20()))
1003       cl->fM20    = clus->GetM20();
1004     if (!TMath::IsNaN(clus->GetM02()))
1005       cl->fM02    = clus->GetM02();
1006     Double_t maxAxis = -1, minAxis = -1;
1007     GetSigma(clus,maxAxis,minAxis);
1008     clus->SetTOF(maxAxis);     // store sigma in TOF for later plotting
1009     cl->fSig      = maxAxis;
1010     Double_t sEtaEta = -1;
1011     Double_t sPhiPhi = -1;
1012     GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
1013     cl->fSigEtaEta = sEtaEta;
1014     cl->fSigPhiPhi = sPhiPhi;
1015     Double_t clusterEcc = -1;
1016     if (maxAxis > 0)
1017       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
1018     clus->SetChi2(clusterEcc); // store ecc in chi2 for later plotting
1019     cl->fEcc        = clusterEcc;
1020     cl->fTrIso      = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1021     cl->fTrIso1     = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1022     cl->fTrIso2     = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
1023     cl->fTrIsoD1    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
1024     cl->fTrIso1D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
1025     cl->fTrIso2D1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
1026     cl->fTrIsoD3    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
1027     cl->fTrIso1D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
1028     cl->fTrIso2D3   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
1029     cl->fTrIsoD4    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2);
1030     cl->fTrIso1D4   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 1);
1031     cl->fTrIso2D4   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 2);
1032     cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
1033     cl->fCeCore     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1034     cl->fCeIso      = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1035     cl->fCeIso1     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
1036     cl->fCeIso3     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
1037     cl->fCeIso4     = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.40);
1038     cl->fCeIso3x3   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 3);
1039     cl->fCeIso4x4   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
1040     cl->fCeIso5x5   = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
1041     cl->fCeIso3x22  = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
1042     cl->fIsShared   = IsShared(clus);
1043     cl->fTrigId     = -1;
1044     cl->fTrigE      = 0;
1045
1046     if (fTriggers) {
1047       Int_t ntrig = fTriggers->GetEntries();
1048       for (Int_t j = 0; j<ntrig; ++j) {
1049         AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1050         if (!sta)
1051           continue;
1052         Short_t idmax = sta->fIdMax;
1053         Bool_t inc = IsIdPartOfCluster(clus, idmax);
1054         if (inc) {
1055           cl->fTrigId     = j;
1056           cl->fTrigE      = sta->fE;
1057           break;
1058         }
1059       }
1060     }
1061
1062     // track matching
1063     if (fDoTrMatGeom) {
1064       Double_t mind2 = 1e10;
1065       for(Int_t j = 0; j<ntrks; ++j) {
1066         AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1067         if (!track)
1068           continue;
1069         if (TMath::Abs(clsEta-track->Eta())>0.5)
1070           continue;
1071         TVector3 vec(clsPos);
1072         Float_t tmpEta=-1, tmpPhi=-1, tmpR2=1e10;
1073         Double_t dedx = -1;
1074         if (fEsdEv) {
1075           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1076           dedx = esdTrack->GetTPCsignal();
1077         }
1078         if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1079           continue;
1080         AliExternalTrackParam tParam;
1081         tParam.CopyFromVTrack(track);
1082         if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpEta, tmpPhi))
1083           continue;
1084         tmpR2 = tmpEta*tmpEta + tmpPhi*tmpPhi;
1085         Double_t d2 = tmpR2;
1086         if (mind2>d2) {
1087           mind2=d2;
1088           cl->fTrDz   = tmpEta;
1089           cl->fTrDr   = tmpPhi;
1090           cl->fTrEp   = clus->E()/track->P();
1091           cl->fTrDedx = dedx;
1092           cl->fIsTrackM = 1;
1093         }
1094       }
1095       if (cl->fIsTrackM) {
1096         fHMatchDr->Fill(cl->fTrDr);
1097         fHMatchDz->Fill(cl->fTrDz);
1098         fHMatchEp->Fill(cl->fTrEp);
1099       }
1100     }
1101
1102     //mc matching
1103     if (fMcParts) {
1104       Int_t nmc = filtMcParts.GetEntries();
1105       Double_t diffR2 = 1e9;
1106       AliStaPart *msta = 0;
1107       for (Int_t j=0; j<nmc; ++j) {
1108         AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1109         Double_t t1=clsVec.Eta()-pa->fVEta;
1110         Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1111         Double_t tmp = t1*t1+t2*t2;
1112         if (tmp<diffR2) {
1113           diffR2 = tmp;
1114           msta   = pa;
1115         }
1116       }
1117       if (diffR2<10 && msta!=0) {
1118         cl->fMcLabel = msta->fLab;
1119       }
1120     }
1121
1122     cl->fEmbE = 0;
1123     if (fDigits && fEmbedMode) {
1124       for(Int_t j=0; j<cl->fN; ++j) {
1125         Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1126         Short_t pos = -1;
1127         std::map<Short_t,Short_t>::iterator it = map.find(cid);
1128         if (it!=map.end())
1129           pos = it->second;
1130         if (pos<0)
1131           continue;
1132         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1133         if (!digit)
1134           continue;
1135         if (digit->GetType()<-1) {
1136           cl->fEmbE += digit->GetChi2();
1137         }
1138       }
1139     }
1140   }
1141 }
1142
1143 //________________________________________________________________________
1144 void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1145 {
1146   // Calculate track properties for primary tracks.
1147
1148   fSelPrimTracks->Clear();
1149
1150   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1151   TClonesArray *tracks = 0;
1152   Bool_t doConstrain = kFALSE;
1153   TString trkname(fPrimTracksName);
1154   if (fEsdEv) {
1155     if (trkname.IsNull()) {
1156       trkname = "Tracks";
1157       am->LoadBranch("Tracks");
1158       fSelPrimTracks->SetOwner(kTRUE);
1159       doConstrain = kTRUE;
1160     }
1161     TList *l = fEsdEv->GetList();
1162     tracks = dynamic_cast<TClonesArray*>(l->FindObject(trkname));
1163   } else {
1164     trkname = "tracks";
1165     am->LoadBranch("tracks");
1166     TList *l = fAodEv->GetList();
1167     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1168   }
1169
1170   if (!tracks) {
1171     AliError(Form("Pointer to tracks %s == 0", trkname.Data() ));
1172     return;
1173   }
1174
1175   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1176   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1177   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
1178
1179   if (fEsdEv) {
1180     am->LoadBranch("PrimaryVertex.");
1181     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1182     am->LoadBranch("SPDVertex.");
1183     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1184     am->LoadBranch("Tracks");
1185     const Int_t Ntracks = tracks->GetEntries();
1186     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1187       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1188       if (!track) {
1189         AliWarning(Form("Could not receive track %d\n", iTracks));
1190         continue;
1191       }
1192       if (fPrimTrCuts && !fPrimTrCuts->IsSelected(track))
1193         continue;
1194       Double_t eta = track->Eta();
1195       if (eta<-1||eta>1)
1196         continue;
1197       if (track->Phi()<phimin||track->Phi()>phimax)
1198         continue;
1199       if (!doConstrain) {
1200         fSelPrimTracks->Add(track);
1201         continue;
1202       }
1203       AliESDtrack copyt(*track);
1204       Double_t bfield[3];
1205       copyt.GetBxByBz(bfield);
1206       AliExternalTrackParam tParam;
1207       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1208       if (!relate)
1209         continue;
1210       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1211       Double_t p[3]      = { 0. };
1212       copyt.GetPxPyPz(p);
1213       Double_t pos[3]    = { 0. };
1214       copyt.GetXYZ(pos);
1215       Double_t covTr[21] = { 0. };
1216       copyt.GetCovarianceXYZPxPyPz(covTr);
1217       Double_t pid[10]   = { 0. };
1218       copyt.GetESDpid(pid);
1219       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1220                                             copyt.GetLabel(),
1221                                             p,
1222                                             kTRUE,
1223                                             pos,
1224                                             kFALSE,
1225                                             covTr,
1226                                             (Short_t)copyt.GetSign(),
1227                                             copyt.GetITSClusterMap(),
1228                                             pid,
1229                                             0,/*fPrimaryVertex,*/
1230                                             kTRUE, // check if this is right
1231                                             vtx->UsesTrack(copyt.GetID()));
1232       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1233       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1234       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1235       if(ndf>0)
1236         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1237       else
1238         aTrack->SetChi2perNDF(-1);
1239       aTrack->SetFlags(copyt.GetStatus());
1240       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
1241       fSelPrimTracks->Add(aTrack);
1242     }
1243   } else {
1244     Int_t ntracks = tracks->GetEntries();
1245     for (Int_t i=0; i<ntracks; ++i) {
1246       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1247       if (!track)
1248         continue;
1249       Double_t eta = track->Eta();
1250       if (eta<-1||eta>1)
1251         continue;
1252       if (track->Phi()<phimin||track->Phi()>phimax)
1253         continue;
1254       if(track->GetTPCNcls()<fMinNClusPerTr)
1255         continue;
1256       //todo: Learn how to set/filter AODs for prim/sec tracks
1257       fSelPrimTracks->Add(track);
1258     }
1259   }
1260 }
1261
1262 //________________________________________________________________________
1263 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
1264 {
1265   // Calculate track properties (including secondaries).
1266
1267   fSelTracks->Clear();
1268
1269   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1270   TClonesArray *tracks = 0;
1271   if (fEsdEv) {
1272     am->LoadBranch("Tracks");
1273     TList *l = fEsdEv->GetList();
1274     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1275   } else {
1276     am->LoadBranch("tracks");
1277     TList *l = fAodEv->GetList();
1278     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1279   }
1280
1281   if (!tracks)
1282     return;
1283
1284   if (fEsdEv) {
1285     const Int_t Ntracks = tracks->GetEntries();
1286     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1287       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1288       if (!track) {
1289         AliWarning(Form("Could not receive track %d\n", iTracks));
1290         continue;
1291       }
1292       if (fTrCuts && !fTrCuts->IsSelected(track))
1293         continue;
1294       Double_t eta = track->Eta();
1295       if (eta<-1||eta>1)
1296         continue;
1297       fSelTracks->Add(track);
1298     }
1299   } else {
1300     Int_t ntracks = tracks->GetEntries();
1301     for (Int_t i=0; i<ntracks; ++i) {
1302       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1303       if (!track)
1304         continue;
1305       Double_t eta = track->Eta();
1306       if (eta<-1||eta>1)
1307         continue;
1308       if(track->GetTPCNcls()<fMinNClusPerTr)
1309         continue;
1310
1311       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
1312         AliExternalTrackParam tParam;
1313         tParam.CopyFromVTrack(track);
1314         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 430, 0.139, 1, kTRUE)) {
1315           Double_t trkPos[3];
1316           tParam.GetXYZ(trkPos);
1317           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1318         }
1319       }
1320       fSelTracks->Add(track);
1321     }
1322   }
1323 }
1324
1325 //________________________________________________________________________
1326 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
1327 {
1328   // Run custer reconstruction afterburner.
1329
1330   AliVCaloCells *cells = fEsdCells;
1331   if (!cells)
1332     cells = fAodCells;
1333
1334   if (!cells)
1335     return;
1336
1337   Int_t ncells = cells->GetNumberOfCells();
1338   if (ncells<=0)
1339     return;
1340
1341   Double_t cellMeanE = 0, cellSigE = 0;
1342   for (Int_t i = 0; i<ncells; ++i) {
1343     Double_t cellE = cells->GetAmplitude(i);
1344     cellMeanE += cellE;
1345     cellSigE += cellE*cellE;
1346   }
1347   cellMeanE /= ncells;
1348   cellSigE /= ncells;
1349   cellSigE -= cellMeanE*cellMeanE;
1350   if (cellSigE<0)
1351     cellSigE = 0;
1352   cellSigE = TMath::Sqrt(cellSigE / ncells);
1353
1354   Double_t subE = cellMeanE - 7*cellSigE;
1355   if (subE<0)
1356     return;
1357
1358   for (Int_t i = 0; i<ncells; ++i) {
1359     Short_t id=-1;
1360     Int_t mclabel = -1;
1361     Double_t amp=0,time=0, efrac = 0;
1362     if (!cells->GetCell(i, id, amp, time,mclabel,efrac))
1363       continue;
1364     amp -= cellMeanE;
1365     if (amp<0.001)
1366       amp = 0;
1367     cells->SetCell(i, id, amp, time);
1368   }
1369
1370   TObjArray *clusters = fEsdClusters;
1371   if (!clusters)
1372     clusters = fAodClusters;
1373   if (!clusters)
1374     return;
1375
1376   Int_t nclus = clusters->GetEntries();
1377   for (Int_t i = 0; i<nclus; ++i) {
1378     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1379     if (!clus->IsEMCAL())
1380       continue;
1381     Int_t nc = clus->GetNCells();
1382     Double_t clusE = 0;
1383     UShort_t ids[100] = {0};
1384     Double_t fra[100] = {0};
1385     for (Int_t j = 0; j<nc; ++j) {
1386       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1387       Double_t cen = cells->GetCellAmplitude(id);
1388       clusE += cen;
1389       if (cen>0) {
1390         ids[nc] = id;
1391         ++nc;
1392       }
1393     }
1394     if (clusE<=0) {
1395       clusters->RemoveAt(i);
1396       continue;
1397     }
1398
1399     for (Int_t j = 0; j<nc; ++j) {
1400       Short_t id = ids[j];
1401       Double_t cen = cells->GetCellAmplitude(id);
1402       fra[j] = cen/clusE;
1403     }
1404     clus->SetE(clusE);
1405     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1406     if (aodclus) {
1407       aodclus->Clear("");
1408       aodclus->SetNCells(nc);
1409       aodclus->SetCellsAmplitudeFraction(fra);
1410       aodclus->SetCellsAbsId(ids);
1411     }
1412   }
1413   clusters->Compress();
1414 }
1415
1416 //________________________________________________________________________
1417 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1418 {
1419   // Fill histograms related to cell properties.
1420
1421   AliVCaloCells *cells = fEsdCells;
1422   if (!cells)
1423     cells = fAodCells;
1424
1425   if (!cells)
1426     return;
1427
1428   Int_t cellModCount[12] = {0};
1429   Double_t cellMaxE = 0;
1430   Double_t cellMeanE = 0;
1431   Int_t ncells = cells->GetNumberOfCells();
1432   if (ncells==0)
1433     return;
1434
1435   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1436
1437   for (Int_t i = 0; i<ncells; ++i) {
1438     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1439     Double_t cellE = cells->GetAmplitude(i);
1440     fHCellE->Fill(cellE);
1441     if (cellE>cellMaxE)
1442       cellMaxE = cellE;
1443     cellMeanE += cellE;
1444
1445     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1446     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1447     if (!ret) {
1448       AliError(Form("Could not get cell index for %d", absID));
1449       continue;
1450     }
1451     ++cellModCount[iSM];
1452     Int_t iPhi=-1, iEta=-1;
1453     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1454     fHColuRow[iSM]->Fill(iEta,iPhi,1);
1455     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1456     fHCellFreqNoCut[iSM]->Fill(absID);
1457     if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1458     if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
1459     if (fHCellCheckE && fHCellCheckE[absID])
1460       fHCellCheckE[absID]->Fill(cellE);
1461     fHCellFreqE[iSM]->Fill(absID, cellE);
1462   }
1463   fHCellH->Fill(cellMaxE);
1464   cellMeanE /= ncells;
1465   fHCellM->Fill(cellMeanE);
1466   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1467   for (Int_t i=0; i<nsm; ++i)
1468     fHCellMult[i]->Fill(cellModCount[i]);
1469 }
1470
1471 //________________________________________________________________________
1472 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1473 {
1474   // Fill histograms related to cluster properties.
1475
1476   TObjArray *clusters = fEsdClusters;
1477   if (!clusters)
1478     clusters = fAodClusters;
1479   if (!clusters)
1480     return;
1481
1482   Double_t vertex[3] = {0};
1483   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1484
1485   Int_t nclus = clusters->GetEntries();
1486   for(Int_t i = 0; i<nclus; ++i) {
1487     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1488     if (!clus)
1489       continue;
1490     if (!clus->IsEMCAL())
1491       continue;
1492     TLorentzVector clusterVec;
1493     clus->GetMomentum(clusterVec,vertex);
1494     Double_t maxAxis    = clus->GetTOF(); //sigma
1495     Double_t clusterEcc = clus->Chi2();   //eccentricity
1496     fHClustEccentricity->Fill(clusterEcc);
1497     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1498     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1499     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1500     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1501     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1502     fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
1503   }
1504 }
1505
1506 //________________________________________________________________________
1507 void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1508 {
1509   // Get Mc truth particle information.
1510
1511   if (!fMcParts)
1512     return;
1513
1514   fMcParts->Clear();
1515
1516   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1517   Double_t etamin = -0.7;
1518   Double_t etamax = +0.7;
1519   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1520   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1521
1522   if (fAodEv) {
1523     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1524     am->LoadBranch(AliAODMCParticle::StdBranchName());
1525     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1526     if (!tca)
1527       return;
1528
1529     Int_t nents = tca->GetEntries();
1530     for(int it=0; it<nents; ++it) {
1531       AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
1532       part->Print();
1533
1534       // pion or eta meson or direct photon
1535       if(part->GetPdgCode() == 111) {
1536       } else if(part->GetPdgCode() == 221) {
1537       } else if(part->GetPdgCode() == 22 ) {
1538       } else
1539         continue;
1540
1541       // primary particle
1542       Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1543       if(dR > 1.0)
1544         continue;
1545
1546       // kinematic cuts
1547       Double_t pt = part->Pt() ;
1548       if (pt<0.5)
1549         continue;
1550       Double_t eta = part->Eta();
1551       if (eta<etamin||eta>etamax)
1552         continue;
1553       Double_t phi  = part->Phi();
1554       if (phi<phimin||phi>phimax)
1555         continue;
1556
1557       ProcessDaughters(part, it, tca);
1558     }
1559     return;
1560   }
1561
1562   AliMCEvent *mcEvent = MCEvent();
1563   if (!mcEvent)
1564     return;
1565
1566   const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1567   if (!evtVtx)
1568     return;
1569
1570   mcEvent->PreReadAll();
1571
1572   Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1573   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1574     AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1575     if (!mcP)
1576       continue;
1577
1578     // pion or eta meson or direct photon
1579     if(mcP->PdgCode() == 111) {
1580     } else if(mcP->PdgCode() == 221) {
1581     } else if(mcP->PdgCode() == 22 ) {
1582     } else
1583       continue;
1584
1585     // primary particle
1586     Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1587                               (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1588     if(dR > 1.0)
1589       continue;
1590
1591     // kinematic cuts
1592     Double_t pt = mcP->Pt() ;
1593     if (pt<0.5)
1594       continue;
1595     Double_t eta = mcP->Eta();
1596     if (eta<etamin||eta>etamax)
1597       continue;
1598     Double_t phi  = mcP->Phi();
1599     if (phi<phimin||phi>phimax)
1600       continue;
1601
1602     ProcessDaughters(mcP, iTrack, mcEvent);
1603   }
1604 }
1605
1606 //________________________________________________________________________
1607 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1608 {
1609   // Fill ntuple.
1610
1611   if (!fNtuple)
1612     return;
1613
1614   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1615   if (fAodEv) {
1616     fHeader->fRun            = fAodEv->GetRunNumber();
1617     fHeader->fOrbit          = fAodEv->GetHeader()->GetOrbitNumber();
1618     fHeader->fPeriod         = fAodEv->GetHeader()->GetPeriodNumber();
1619     fHeader->fBx             = fAodEv->GetHeader()->GetBunchCrossNumber();
1620     fHeader->fL0             = fAodEv->GetHeader()->GetL0TriggerInputs();
1621     fHeader->fL1             = fAodEv->GetHeader()->GetL1TriggerInputs();
1622     fHeader->fL2             = fAodEv->GetHeader()->GetL2TriggerInputs();
1623     fHeader->fTrClassMask    = fAodEv->GetHeader()->GetTriggerMask();
1624     fHeader->fTrCluster      = fAodEv->GetHeader()->GetTriggerCluster();
1625     fHeader->fOffTriggers    = fAodEv->GetHeader()->GetOfflineTrigger();
1626     fHeader->fFiredTriggers  = fAodEv->GetHeader()->GetFiredTriggerClasses();
1627   } else {
1628     fHeader->fRun            = fEsdEv->GetRunNumber();
1629     fHeader->fOrbit          = fEsdEv->GetHeader()->GetOrbitNumber();
1630     fHeader->fPeriod         = fEsdEv->GetESDRun()->GetPeriodNumber();
1631     fHeader->fBx             = fEsdEv->GetHeader()->GetBunchCrossNumber();
1632     fHeader->fL0             = fEsdEv->GetHeader()->GetL0TriggerInputs();
1633     fHeader->fL1             = fEsdEv->GetHeader()->GetL1TriggerInputs();
1634     fHeader->fL2             = fEsdEv->GetHeader()->GetL2TriggerInputs();
1635     fHeader->fTrClassMask    = fEsdEv->GetHeader()->GetTriggerMask();
1636     fHeader->fTrCluster      = fEsdEv->GetHeader()->GetTriggerCluster();
1637     fHeader->fOffTriggers    = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1638     fHeader->fFiredTriggers  = fEsdEv->GetFiredTriggerClasses();
1639     Float_t v0CorrR = 0;
1640     fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1641     const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1642     if (mult)
1643       fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1644     fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
1645     AliTriggerAnalysis trAn; /// Trigger Analysis
1646     Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1647     Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
1648     fHeader->fV0And        = v0A && v0B;
1649     fHeader->fIsHT         = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1650     am->LoadBranch("SPDPileupVertices");
1651     am->LoadBranch("TrkPileupVertices");
1652     fHeader->fIsPileup     = fEsdEv->IsPileupFromSPD(3,0.8);
1653     fHeader->fIsPileup2    = fEsdEv->IsPileupFromSPD(3,0.4);
1654     fHeader->fIsPileup4    = fEsdEv->IsPileupFromSPD(3,0.2);
1655     fHeader->fIsPileup8    = fEsdEv->IsPileupFromSPD(3,0.1);
1656     fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1657     fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
1658   }
1659
1660   AliCentrality *cent = InputEvent()->GetCentrality();
1661   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
1662   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
1663   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
1664   fHeader->fCqual     = cent->GetQuality();
1665
1666   AliEventplane *ep = InputEvent()->GetEventplane();
1667   if (ep) {
1668     if (ep->GetQVector())
1669       fHeader->fPsi     = ep->GetQVector()->Phi()/2. ;
1670     else
1671       fHeader->fPsi = -1;
1672     if (ep->GetQsub1()&&ep->GetQsub2())
1673       fHeader->fPsiRes  = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1674     else
1675       fHeader->fPsiRes = 0;
1676   }
1677
1678   Double_t val = 0;
1679   TString trgclasses(fHeader->fFiredTriggers);
1680   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1681     const char *name = fTrClassNamesArr->At(j)->GetName();
1682     TRegexp regexp(name);
1683     if (trgclasses.Contains(regexp))
1684       val += TMath::Power(2,j);
1685   }
1686   fHeader->fTcls = (UInt_t)val;
1687
1688   fHeader->fNSelTr     = fSelTracks->GetEntries();
1689   fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1690   fHeader->fNSelPrimTr1   = 0;
1691   fHeader->fNSelPrimTr2   = 0;
1692   for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1693     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1694     if(track->Pt()>1)
1695       ++fHeader->fNSelPrimTr1;
1696     if(track->Pt()>2)
1697       ++fHeader->fNSelPrimTr2;
1698   }
1699
1700   fHeader->fNCells   = 0;
1701   fHeader->fNCells0  = 0;
1702   fHeader->fNCells01 = 0;
1703   fHeader->fNCells03 = 0;
1704   fHeader->fNCells1  = 0;
1705   fHeader->fNCells2  = 0;
1706   fHeader->fNCells5  = 0;
1707   fHeader->fMaxCellE = 0;
1708
1709   AliVCaloCells *cells = fEsdCells;
1710   if (!cells)
1711     cells = fAodCells;
1712
1713   if (cells) {
1714     Int_t ncells = cells->GetNumberOfCells();
1715     for(Int_t j=0; j<ncells; ++j) {
1716       Double_t cellen = cells->GetAmplitude(j);
1717       if (cellen>0.045)
1718         ++fHeader->fNCells0;
1719       if (cellen>0.1)
1720         ++fHeader->fNCells01;
1721       if (cellen>0.3)
1722         ++fHeader->fNCells03;
1723       if (cellen>1)
1724         ++fHeader->fNCells1;
1725       if (cellen>2)
1726         ++fHeader->fNCells2;
1727       if (cellen>5)
1728         ++fHeader->fNCells5;
1729       if (cellen>fHeader->fMaxCellE)
1730         fHeader->fMaxCellE = cellen;
1731     }
1732     fHeader->fNCells = ncells;
1733   }
1734
1735   fHeader->fNClus      = 0;
1736   fHeader->fNClus1     = 0;
1737   fHeader->fNClus2     = 0;
1738   fHeader->fNClus5     = 0;
1739   fHeader->fMaxClusE   = 0;
1740
1741   TObjArray *clusters = fEsdClusters;
1742   if (!clusters)
1743     clusters = fAodClusters;
1744
1745   if (clusters) {
1746     Int_t nclus = clusters->GetEntries();
1747     for(Int_t j=0; j<nclus; ++j) {
1748       AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1749       if (!clus->IsEMCAL())
1750         continue;
1751       Double_t clusen = clus->E();
1752       if (clusen>1)
1753         ++fHeader->fNClus1;
1754       if (clusen>2)
1755         ++fHeader->fNClus2;
1756       if (clusen>5)
1757         ++fHeader->fNClus5;
1758       if (clusen>fHeader->fMaxClusE)
1759         fHeader->fMaxClusE = clusen;
1760     }
1761     fHeader->fNClus = nclus;
1762   }
1763
1764   fHeader->fMaxTrE     = 0;
1765   if (fTriggers) {
1766     Int_t ntrig = fTriggers->GetEntries();
1767     for (Int_t j = 0; j<ntrig; ++j) {
1768       AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1769       if (!sta)
1770         continue;
1771       if (sta->fE>fHeader->fMaxTrE)
1772         fHeader->fMaxTrE = sta->fE;
1773     }
1774   }
1775
1776   // count cells above 100 MeV on super modules
1777   fHeader->fNcSM0 = GetNCells(0, 0.100);
1778   fHeader->fNcSM1 = GetNCells(1, 0.100);
1779   fHeader->fNcSM2 = GetNCells(2, 0.100);
1780   fHeader->fNcSM3 = GetNCells(3, 0.100);
1781   fHeader->fNcSM4 = GetNCells(4, 0.100);
1782   fHeader->fNcSM5 = GetNCells(5, 0.100);
1783   fHeader->fNcSM6 = GetNCells(6, 0.100);
1784   fHeader->fNcSM7 = GetNCells(7, 0.100);
1785   fHeader->fNcSM8 = GetNCells(8, 0.100);
1786   fHeader->fNcSM9 = GetNCells(9, 0.100);
1787
1788   if (fAodEv) {
1789     am->LoadBranch("vertices");
1790     AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1791     FillVertex(fPrimVert, pv);
1792     AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1793     FillVertex(fSpdVert, sv);
1794   } else {
1795     am->LoadBranch("PrimaryVertex.");
1796     const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1797     FillVertex(fPrimVert, pv);
1798     am->LoadBranch("SPDVertex.");
1799     const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1800     FillVertex(fSpdVert, sv);
1801     am->LoadBranch("TPCVertex.");
1802     const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1803     FillVertex(fTpcVert, tv);
1804   }
1805
1806   fNtuple->Fill();
1807 }
1808
1809 //________________________________________________________________________
1810 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1811 {
1812   // Fill histograms related to pions.
1813
1814   Double_t vertex[3] = {0};
1815   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1816
1817   TObjArray *clusters = fEsdClusters;
1818   if (!clusters)
1819     clusters = fAodClusters;
1820
1821   if (clusters) {
1822     TLorentzVector clusterVec1;
1823     TLorentzVector clusterVec2;
1824     TLorentzVector pionVec;
1825
1826     Int_t nclus = clusters->GetEntries();
1827     for (Int_t i = 0; i<nclus; ++i) {
1828       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1829       if (!clus1)
1830         continue;
1831       if (!clus1->IsEMCAL())
1832         continue;
1833       if (clus1->E()<fMinE)
1834         continue;
1835       if (clus1->GetNCells()<fNminCells)
1836         continue;
1837       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1838         continue;
1839       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1840         continue;
1841       clus1->GetMomentum(clusterVec1,vertex);
1842       for (Int_t j = i+1; j<nclus; ++j) {
1843         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1844         if (!clus2)
1845           continue;
1846         if (!clus2->IsEMCAL())
1847           continue;
1848         if (clus2->E()<fMinE)
1849           continue;
1850         if (clus2->GetNCells()<fNminCells)
1851           continue;
1852         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1853           continue;
1854         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1855           continue;
1856         clus2->GetMomentum(clusterVec2,vertex);
1857         pionVec = clusterVec1 + clusterVec2;
1858         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1859         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1860         if (pionZgg < fAsymMax) {
1861           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1862           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1863           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1864           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1865           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1866           fHPionInvMasses[bin]->Fill(pionVec.M());
1867         }
1868       }
1869     }
1870   }
1871 }
1872
1873 //________________________________________________________________________
1874 void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1875 {
1876   // Fill additional MC information histograms.
1877
1878   if (!fMcParts)
1879     return;
1880
1881   // check if aod or esd mc mode and the fill histos
1882 }
1883
1884 //________________________________________________________________________
1885 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1886 {
1887   // Fill other histograms.
1888
1889   for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1890     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1891     if(!track)
1892       continue;
1893     fHPrimTrackPt->Fill(track->Pt());
1894     fHPrimTrackEta->Fill(track->Eta());
1895     fHPrimTrackPhi->Fill(track->Phi());
1896   }
1897 }
1898
1899 //________________________________________________________________________
1900 void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1901 {
1902   // Fill track histograms.
1903
1904   if (fSelPrimTracks) {
1905     for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1906       AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1907       if(!track)
1908         continue;
1909       fHPrimTrackPt->Fill(track->Pt());
1910       fHPrimTrackEta->Fill(track->Eta());
1911       fHPrimTrackPhi->Fill(track->Phi());
1912     }
1913   }
1914 }
1915
1916 //__________________________________________________________________________________________________
1917 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1918 {
1919   // Fill vertex from ESD vertex info.
1920
1921   v->fVx   = esdv->GetX();
1922   v->fVy   = esdv->GetY();
1923   v->fVz   = esdv->GetZ();
1924   v->fVc   = esdv->GetNContributors();
1925   v->fDisp = esdv->GetDispersion();
1926   v->fZres = esdv->GetZRes();
1927   v->fChi2 = esdv->GetChi2();
1928   v->fSt   = esdv->GetStatus();
1929   v->fIs3D = esdv->IsFromVertexer3D();
1930   v->fIsZ  = esdv->IsFromVertexerZ();
1931 }
1932
1933 //__________________________________________________________________________________________________
1934 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1935 {
1936   // Fill vertex from AOD vertex info.
1937
1938   v->fVx   = aodv->GetX();
1939   v->fVy   = aodv->GetY();
1940   v->fVz   = aodv->GetZ();
1941   v->fVc   = aodv->GetNContributors();
1942   v->fChi2 = aodv->GetChi2();
1943 }
1944
1945 //________________________________________________________________________
1946 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1947 {
1948   // Compute isolation based on cell content.
1949
1950   AliVCaloCells *cells = fEsdCells;
1951   if (!cells)
1952     cells = fAodCells;
1953   if (!cells)
1954     return 0;
1955
1956   Double_t cellIsolation = 0;
1957   Double_t rad2 = radius*radius;
1958   Int_t ncells = cells->GetNumberOfCells();
1959   for (Int_t i = 0; i<ncells; ++i) {
1960     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1961     Float_t eta=-1, phi=-1;
1962     fGeom->EtaPhiFromIndex(absID,eta,phi);
1963     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1964     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1965     if(dist>rad2)
1966       continue;
1967     Double_t cellE = cells->GetAmplitude(i);
1968     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1969     Double_t cellEt = cellE*sin(theta);
1970     cellIsolation += cellEt;
1971   }
1972   return cellIsolation;
1973 }
1974
1975 //________________________________________________________________________
1976 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1977 {
1978   // Compute isolation based on cell content, in a NxM rectangle.
1979
1980   AliVCaloCells *cells = fEsdCells;
1981   if (!cells)
1982     cells = fAodCells;
1983   if (!cells)
1984     return 0;
1985
1986   Double_t cellIsolation = 0;
1987   Int_t ncells = cells->GetNumberOfCells();
1988   for (Int_t i = 0; i<ncells; ++i) {
1989     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1990     Float_t eta=-1, phi=-1;
1991     fGeom->EtaPhiFromIndex(absID,eta,phi);
1992     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1993     Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1994     if(TMath::Abs(etadiff)/0.014>N)
1995       continue;
1996     if(TMath::Abs(phidiff)/0.014>M)
1997       continue;
1998     Double_t cellE = cells->GetAmplitude(i);
1999     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
2000     Double_t cellEt = cellE*sin(theta);
2001     cellIsolation += cellEt;
2002   }
2003   return cellIsolation;
2004 }
2005
2006 //________________________________________________________________________
2007 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
2008 {
2009   // Get maximum energy of attached cell.
2010
2011   Double_t ret = 0;
2012   Int_t ncells = cluster->GetNCells();
2013   if (fEsdCells) {
2014     for (Int_t i=0; i<ncells; i++) {
2015       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2016       ret += e;
2017       }
2018   } else {
2019     for (Int_t i=0; i<ncells; i++) {
2020       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2021       ret += e;
2022     }
2023   }
2024   return ret;
2025 }
2026
2027 //________________________________________________________________________
2028 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
2029 {
2030   // Get maximum energy of attached cell.
2031
2032   id = -1;
2033   Double_t maxe = 0;
2034   Int_t ncells = cluster->GetNCells();
2035   if (fEsdCells) {
2036     for (Int_t i=0; i<ncells; i++) {
2037       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2038       if (e>maxe) {
2039         maxe = e;
2040         id   = cluster->GetCellAbsId(i);
2041       }
2042     }
2043   } else {
2044     for (Int_t i=0; i<ncells; i++) {
2045       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2046       if (e>maxe)
2047         maxe = e;
2048         id   = cluster->GetCellAbsId(i);
2049     }
2050   }
2051   return maxe;
2052 }
2053
2054 //________________________________________________________________________
2055 Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
2056 {
2057   // Get second maximum cell.
2058
2059   AliVCaloCells *cells = fEsdCells;
2060   if (!cells)
2061     cells = fAodCells;
2062   if (!cells)
2063     return -1;
2064
2065   Double_t secondEmax=0, firstEmax=0;
2066   Double_t cellen;
2067   for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2068     Int_t absId = clus->GetCellAbsId(iCell);
2069     cellen = cells->GetCellAmplitude(absId);
2070     if(cellen > firstEmax)
2071       firstEmax = cellen;
2072   }
2073   for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2074     Int_t absId = clus->GetCellAbsId(iCell);
2075     cellen = cells->GetCellAmplitude(absId);
2076     if(cellen < firstEmax && cellen > secondEmax) {
2077       secondEmax = cellen;
2078       id = absId;
2079     }
2080   }
2081   return secondEmax;
2082 }
2083
2084 //________________________________________________________________________
2085 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
2086 {
2087   // Calculate the (E) weighted variance along the longer (eigen) axis.
2088
2089   sigmaMax = 0;          // cluster variance along its longer axis
2090   sigmaMin = 0;          // cluster variance along its shorter axis
2091   Double_t Ec  = c->E(); // cluster energy
2092   if(Ec<=0)
2093     return;
2094   Double_t Xc  = 0 ;     // cluster first moment along X
2095   Double_t Yc  = 0 ;     // cluster first moment along Y
2096   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
2097   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
2098   Double_t Syy = 0 ;     // cluster covariance^2
2099
2100   AliVCaloCells *cells = fEsdCells;
2101   if (!cells)
2102     cells = fAodCells;
2103
2104   if (!cells)
2105     return;
2106
2107   Int_t ncells = c->GetNCells();
2108   if (ncells==1)
2109     return;
2110
2111   for(Int_t j=0; j<ncells; ++j) {
2112     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2113     Double_t cellen = cells->GetCellAmplitude(id);
2114     TVector3 pos;
2115     fGeom->GetGlobal(id,pos);
2116     Xc  += cellen*pos.X();
2117     Yc  += cellen*pos.Y();
2118     Sxx += cellen*pos.X()*pos.X();
2119     Syy += cellen*pos.Y()*pos.Y();
2120     Sxy += cellen*pos.X()*pos.Y();
2121   }
2122   Xc  /= Ec;
2123   Yc  /= Ec;
2124   Sxx /= Ec;
2125   Syy /= Ec;
2126   Sxy /= Ec;
2127   Sxx -= Xc*Xc;
2128   Syy -= Yc*Yc;
2129   Sxy -= Xc*Yc;
2130   Sxx = TMath::Abs(Sxx);
2131   Syy = TMath::Abs(Syy);
2132   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2133   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2134   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2135   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
2136 }
2137
2138 //________________________________________________________________________
2139 void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2140 {
2141   // Calculate the (E) weighted variance along the pseudorapidity.
2142
2143   sEtaEta = 0;
2144   sPhiPhi = 0;
2145
2146   Double_t Ec  = c->E(); // cluster energy
2147   if(Ec<=0)
2148     return;
2149
2150   const Int_t ncells = c->GetNCells();
2151
2152   Double_t EtaC    = 0;  // cluster first moment along eta
2153   Double_t PhiC    = 0;  // cluster first moment along phi
2154   Double_t Setaeta = 0;  // cluster second central moment along eta
2155   Double_t Sphiphi = 0;  // cluster second central moment along phi
2156   Double_t w[ncells];    // weight max(0,4.5*log(E_i/Ec))
2157   Double_t sumw = 0;
2158   Int_t id[ncells];
2159
2160   AliVCaloCells *cells = fEsdCells;
2161   if (!cells)
2162     cells = fAodCells;
2163
2164   if (!cells)
2165     return;
2166
2167   if (ncells==1)
2168     return;
2169
2170   for(Int_t j=0; j<ncells; ++j) {
2171     id[j] = TMath::Abs(c->GetCellAbsId(j));
2172     Double_t cellen = cells->GetCellAmplitude(id[j]);
2173     w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2174     TVector3 pos;
2175     fGeom->GetGlobal(id[j],pos);
2176     EtaC += w[j]*pos.Eta();
2177     PhiC += w[j]*pos.Phi();
2178     sumw += w[j];
2179   }
2180   EtaC /= sumw;
2181   PhiC /= sumw;
2182
2183   for(Int_t j=0; j<ncells; ++j) {
2184     TVector3 pos;
2185     fGeom->GetGlobal(id[j],pos);
2186     Setaeta =  w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2187     Sphiphi =  w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2188   }
2189   Setaeta /= sumw;
2190   sEtaEta = TMath::Sqrt(Setaeta);
2191   Sphiphi /= sumw;
2192   sPhiPhi = TMath::Sqrt(Sphiphi);
2193 }
2194
2195 //________________________________________________________________________
2196 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
2197 {
2198   // Calculate number of attached cells above emin.
2199
2200   AliVCaloCells *cells = fEsdCells;
2201   if (!cells)
2202     cells = fAodCells;
2203   if (!cells)
2204     return 0;
2205
2206   Int_t n = 0;
2207   Int_t ncells = c->GetNCells();
2208   for(Int_t j=0; j<ncells; ++j) {
2209     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2210     Double_t cellen = cells->GetCellAmplitude(id);
2211     if (cellen>=emin)
2212       ++n;
2213   }
2214   return n;
2215 }
2216
2217 //________________________________________________________________________
2218 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2219 {
2220   // Calculate number of cells per SM above emin.
2221
2222   AliVCaloCells *cells = fEsdCells;
2223   if (!cells)
2224     cells = fAodCells;
2225   if (!cells)
2226     return 0;
2227
2228   Int_t n = 0;
2229   Int_t ncells = cells->GetNumberOfCells();
2230   for(Int_t j=0; j<ncells; ++j) {
2231     Int_t id = TMath::Abs(cells->GetCellNumber(j));
2232     Double_t cellen = cells->GetCellAmplitude(id);
2233     if (cellen<emin)
2234       continue;
2235     Int_t fsm = fGeom->GetSuperModuleNumber(id);
2236     if (fsm != sm)
2237       continue;
2238     ++n;
2239   }
2240   return n;
2241 }
2242
2243 //________________________________________________________________________
2244 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
2245 {
2246   // Compute isolation based on tracks.
2247
2248   Double_t trkIsolation = 0;
2249   Double_t rad2 = radius*radius;
2250   Int_t ntrks = fSelPrimTracks->GetEntries();
2251   for(Int_t j = 0; j<ntrks; ++j) {
2252     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
2253     if (!track)
2254       continue;
2255     if (track->Pt()<pt)
2256       continue;
2257     Float_t eta = track->Eta();
2258     Float_t phi = track->Phi();
2259     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2260     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2261     if(dist>rad2)
2262       continue;
2263     trkIsolation += track->Pt();
2264   }
2265   return trkIsolation;
2266 }
2267
2268 //________________________________________________________________________
2269 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2270 {
2271   // Compute isolation based on tracks.
2272
2273   Double_t trkIsolation = 0;
2274   Int_t ntrks = fSelPrimTracks->GetEntries();
2275   for(Int_t j = 0; j<ntrks; ++j) {
2276     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
2277     if (!track)
2278       continue;
2279     if (track->Pt()<pt)
2280       continue;
2281     Float_t eta = track->Eta();
2282     Float_t phi = track->Phi();
2283     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2284     Double_t etadiff = (eta-cEta);
2285     if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2286       continue;
2287     trkIsolation += track->Pt();
2288   }
2289   return trkIsolation;
2290 }
2291
2292 //________________________________________________________________________
2293 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2294 {
2295   // Returns if cluster shared across super modules.
2296
2297   if (!c)
2298     return 0;
2299
2300   Int_t n = -1;
2301   Int_t ncells = c->GetNCells();
2302   for(Int_t j=0; j<ncells; ++j) {
2303     Int_t id = TMath::Abs(c->GetCellAbsId(j));
2304     Int_t got = id / (24*48);
2305     if (n==-1) {
2306       n = got;
2307       continue;
2308     }
2309     if (got!=n)
2310       return 1;
2311   }
2312   return 0;
2313 }
2314
2315 //________________________________________________________________________
2316 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2317 {
2318   // Returns if id is part of cluster.
2319
2320   AliVCaloCells *cells = fEsdCells;
2321   if (!cells)
2322     cells = fAodCells;
2323   if (!cells)
2324     return 0;
2325
2326   Int_t ncells = c->GetNCells();
2327   for(Int_t j=0; j<ncells; ++j) {
2328     Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2329     if (cid == id)
2330       return 1;
2331   }
2332   return 0;
2333 }
2334
2335 //________________________________________________________________________
2336 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
2337 {
2338   // Print recursively daughter information.
2339
2340   if (!p || !arr)
2341     return;
2342
2343   const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2344   if (!amc)
2345     return;
2346   for (Int_t i=0; i<level; ++i) printf(" ");
2347   amc->Print();
2348
2349   Int_t n = amc->GetNDaughters();
2350   for (Int_t i=0; i<n; ++i) {
2351     Int_t d = amc->GetDaughter(i);
2352     const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2353     PrintDaughters(dmc,arr,level+1);
2354   }
2355 }
2356
2357 //________________________________________________________________________
2358 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2359 {
2360   // Print recursively daughter information.
2361
2362   if (!p || !arr)
2363     return;
2364
2365   for (Int_t i=0; i<level; ++i) printf(" ");
2366   Int_t d1 = p->GetFirstDaughter();
2367   Int_t d2 = p->GetLastDaughter();
2368   printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2369          p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2370   if (d1<0)
2371     return;
2372   if (d2<0)
2373     d2=d1;
2374   for (Int_t i=d1;i<=d2;++i) {
2375     const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
2376     PrintDaughters(dmc,arr,level+1);
2377   }
2378 }
2379
2380 //________________________________________________________________________
2381 void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2382 {
2383   // Print track ref array.
2384
2385   if (!p)
2386     return;
2387
2388   Int_t n = p->GetNumberOfTrackReferences();
2389   for (Int_t i=0; i<n; ++i) {
2390     AliTrackReference *ref = p->GetTrackReference(i);
2391     if (!ref)
2392       continue;
2393     ref->SetUserId(ref->DetectorId());
2394     ref->Print();
2395   }
2396 }
2397
2398 //________________________________________________________________________
2399 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2400 {
2401   // Process and create daughters.
2402
2403   if (!p || !arr)
2404     return;
2405
2406   AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2407   if (!amc)
2408     return;
2409
2410   //amc->Print();
2411
2412   Int_t nparts = arr->GetEntries();
2413   Int_t nents  = fMcParts->GetEntries();
2414
2415   AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2416   newp->fPt  = amc->Pt();
2417   newp->fEta = amc->Eta();
2418   newp->fPhi = amc->Phi();
2419   if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2420     TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2421     newp->fVR = vec.Perp();
2422     newp->fVEta = vec.Eta();
2423     newp->fVPhi = vec.Phi();
2424   }
2425   newp->fPid  = amc->PdgCode();
2426   newp->fLab  = nents;
2427   Int_t moi = amc->GetMother();
2428   if (moi>=0&&moi<nparts) {
2429     const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2430     moi = mmc->GetUniqueID();
2431   }
2432   newp->fMo = moi;
2433   p->SetUniqueID(nents);
2434
2435   // TODO: Determine which detector was hit
2436   //newp->fDet = ???
2437
2438   Int_t n = amc->GetNDaughters();
2439   for (Int_t i=0; i<n; ++i) {
2440     Int_t d = amc->GetDaughter(i);
2441     if (d<=index || d>=nparts)
2442       continue;
2443     AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2444     ProcessDaughters(dmc,d,arr);
2445   }
2446 }
2447
2448 //________________________________________________________________________
2449 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2450 {
2451   // Process and create daughters.
2452
2453   if (!p || !arr)
2454     return;
2455
2456   Int_t d1 = p->GetFirstDaughter();
2457   Int_t d2 = p->GetLastDaughter();
2458   if (0) {
2459     printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2460            index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
2461     PrintTrackRefs(p);
2462   }
2463   Int_t nents  = fMcParts->GetEntries();
2464
2465   AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2466   newp->fPt  = p->Pt();
2467   newp->fEta = p->Eta();
2468   newp->fPhi = p->Phi();
2469   if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2470     TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2471     newp->fVR = vec.Perp();
2472     newp->fVEta = vec.Eta();
2473     newp->fVPhi = vec.Phi();
2474   }
2475   newp->fPid  = p->PdgCode();
2476   newp->fLab  = nents;
2477   Int_t moi = p->GetMother();
2478   if (moi>=0) {
2479     const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2480     moi = mmc->GetUniqueID();
2481   }
2482   newp->fMo = moi;
2483   p->SetUniqueID(nents);
2484
2485   Int_t nref = p->GetNumberOfTrackReferences();
2486   if (nref>0) {
2487     AliTrackReference *ref = p->GetTrackReference(nref-1);
2488     if (ref) {
2489       newp->fDet = ref->DetectorId();
2490     }
2491   }
2492
2493   if (d1<0)
2494     return;
2495   if (d2<0)
2496     d2=d1;
2497   for (Int_t i=d1;i<=d2;++i) {
2498     AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2499     if (dmc->P()<0.01)
2500       continue;
2501     ProcessDaughters(dmc,i,arr);
2502   }
2503 }
2504
2505 //__________________________________________________________________________________________________
2506 void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2507 {
2508   // Calculate momentum.
2509
2510   TVector3 pos;
2511   pos.SetPtEtaPhi(fR,fEta,fPhi);
2512
2513   if(vertex){ //calculate direction relative to vertex
2514     pos -= vertex;
2515   }
2516
2517   Double_t r = pos.Mag();
2518   p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2519 }
2520
2521 //__________________________________________________________________________________________________
2522 void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2523 {
2524   // Calculate momentum.
2525
2526   Double_t v[3] = {0,0,0};
2527   if (vertex) {
2528     v[0] = vertex->fVx;
2529     v[1] = vertex->fVy;
2530     v[2] = vertex->fVz;
2531   }
2532   GetMom(p, v);
2533 }