PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0PbPb.cxx

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