PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx

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