PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.cxx

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