PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx

   1 // $Id$
   2
   3 #include "AliAnalysisTaskEMCALPi0PbPb.h"
   4 #include <TAxis.h>
   5 #include <TChain.h>
   6 #include <TClonesArray.h>
   7 #include <TFile.h>
   8 #include <TGeoGlobalMagField.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TNtuple.h>
  14 #include <TProfile.h>
  15 #include <TString.h>
  16 #include <TVector2.h>
  17 #include "AliAODEvent.h"
  18 #include "AliAODVertex.h"
  19 #include "AliAnalysisManager.h"
  20 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  21 #include "AliCDBManager.h"
  22 #include "AliCentrality.h"
  23 #include "AliEMCALGeoUtils.h"
  24 #include "AliEMCALGeometry.h"
  25 #include "AliEMCALRecoUtils.h"
  26 #include "AliESDCaloTrigger.h"
  27 #include "AliESDEvent.h"
  28 #include "AliESDVertex.h"
  29 #include "AliESDtrack.h"
  30 #include "AliESDtrackCuts.h"
  31 #include "AliEventplane.h"
  32 #include "AliGeomManager.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliLog.h"
  35 #include "AliMagF.h"
  36 #include "AliTrackerBase.h"
  37
  38 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  42   : AliAnalysisTaskSE(name),
  43     fCentVar("V0M"),
  44     fCentFrom(0),
  45     fCentTo(100),
  46     fVtxZMin(-10),
  47     fVtxZMax(+10),
  48     fUseQualFlag(1),
  49     fClusName(),
  50     fDoNtuple(0),
  51     fDoAfterburner(0),
  52     fAsymMax(1),
  53     fNminCells(1),
  54     fMinE(0.100),
  55     fMinErat(0),
  56     fMinEcc(0),
  57     fGeoName("EMCAL_FIRSTYEARV1"),
  58     fMinNClusPerTr(50),
  59     fIsoDist(0.2),
  60     fTrClassNames(""),
  61     fTrCuts(0),
  62     fPrimTrCuts(0),
  63     fDoTrMatGeom(0),
  64     fTrainMode(0),
  65     fMarkCells(),
  66     fMinL0Time(-1),
  67     fMaxL0Time(1024),
  68     fIsGeoMatsSet(0),
  69     fNEvs(0),
  70     fGeom(0),
  71     fReco(0),
  72     fOutput(0),
  73     fTrClassNamesArr(0),
  74     fEsdEv(0),
  75     fAodEv(0),
  76     fRecPoints(0),
  77     fEsdClusters(0),
  78     fEsdCells(0),
  79     fAodClusters(0),
  80     fAodCells(0),
  81     fPtRanges(0),
  82     fSelTracks(0),
  83     fSelPrimTracks(0),
  84     fNAmpInTrigger(0),
  85     fAmpInTrigger(0),
  86     fNtuple(0),
  87     fHeader(0),
  88     fPrimVert(0),
  89     fSpdVert(0),
  90     fTpcVert(0),
  91     fClusters(0),
  92     fTriggers(0),
  93     fHCuts(0x0),
  94     fHVertexZ(0x0),
  95     fHVertexZ2(0x0),
  96     fHCent(0x0),
  97     fHCentQual(0x0),
  98     fHTclsBeforeCuts(0x0),
  99     fHTclsAfterCuts(0x0),
 100     fHColuRow(0x0),
 101     fHColuRowE(0x0),
 102     fHCellMult(0x0),
 103     fHCellE(0x0),
 104     fHCellH(0x0),
 105     fHCellM(0x0),
 106     fHCellM2(0x0),
 107     fHCellFreqNoCut(0x0),
 108     fHCellFreqCut100M(0x0),
 109     fHCellFreqCut300M(0x0),
 110     fHCellFreqE(0x0),
 111     fHCellCheckE(0x0),
 112     fHClustEccentricity(0),
 113     fHClustEtaPhi(0x0),
 114     fHClustEnergyPt(0x0),
 115     fHClustEnergySigma(0x0),
 116     fHClustSigmaSigma(0x0),
 117     fHClustNCellEnergyRatio(0x0),
 118     fHMatchDr(0x0),
 119     fHMatchDz(0x0),
 120     fHMatchEp(0x0),
 121     fHPionEtaPhi(0x0),
 122     fHPionMggPt(0x0),
 123     fHPionMggAsym(0x0),
 124     fHPionMggDgg(0x0)
 125 {
 126   // Constructor.
 127
 128   if (!name)
 129     return;
 130   SetName(name);
 131   DefineInput(0, TChain::Class());
 132   DefineOutput(1, TList::Class());
 133   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks  EMCALTrigger."
 134                "AOD:header,vertices,emcalCells,tracks";
 135 }
 136
 137 //________________________________________________________________________
 138 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
 139 {
 140   // Destructor.
 141
 142   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 143     delete fOutput; fOutput = 0;
 144   }
 145   delete fPtRanges; fPtRanges = 0;
 146   delete fGeom; fGeom = 0;
 147   delete fReco; fReco = 0;
 148   delete fTrClassNamesArr;
 149   delete fSelTracks;
 150   delete fSelPrimTracks;
 151   delete [] fAmpInTrigger;
 152   delete [] fHColuRow;
 153   delete [] fHColuRowE;
 154   delete [] fHCellMult;
 155   delete [] fHCellFreqNoCut;
 156   delete [] fHCellFreqCut100M;
 157   delete [] fHCellFreqCut300M;
 158 }
 159
 160 //________________________________________________________________________
 161 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 162 {
 163   // Create user objects here.
 164
 165   cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
 166   cout << " fCentVar:       " << fCentVar << endl;
 167   cout << " fCentFrom:      " << fCentFrom << endl;
 168   cout << " fCentTo:        " << fCentTo << endl;
 169   cout << " fVtxZMin:       " << fVtxZMin << endl;
 170   cout << " fVtxZMax:       " << fVtxZMax << endl;
 171   cout << " fUseQualFlag:   " << fUseQualFlag << endl;
 172   cout << " fClusName:      \"" << fClusName << "\"" << endl;
 173   cout << " fDoNtuple:      " << fDoNtuple << endl;
 174   cout << " fDoAfterburner: " << fDoAfterburner << endl;
 175   cout << " fAsymMax:       " << fAsymMax << endl;
 176   cout << " fNminCells:     " << fNminCells << endl;
 177   cout << " fMinE:          " << fMinE << endl;
 178   cout << " fMinErat:       " << fMinErat << endl;
 179   cout << " fMinEcc:        " << fMinEcc << endl;
 180   cout << " fGeoName:       \"" << fGeoName << "\"" << endl;
 181   cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
 182   cout << " fTrClassNames:  \"" << fTrClassNames << "\"" << endl;
 183   cout << " fDoNtuple:      " << fDoNtuple << endl;
 184   cout << " fTrCuts:        " << fTrCuts << endl;
 185   cout << " fPrimTrCuts:    " << fPrimTrCuts << endl;
 186   cout << " fDoTrMatGeom:   " << fDoTrMatGeom << endl;
 187   cout << " fTrainMode:     " << fTrainMode << endl;
 188   cout << " fMarkCells:     " << fMarkCells << endl;
 189   cout << " fMinL0Time:     " << fMinL0Time << endl;
 190   cout << " fMaxL0Time:     " << fMaxL0Time << endl;
 191
 192   fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
 193   fReco = new AliEMCALRecoUtils();
 194   fTrClassNamesArr = fTrClassNames.Tokenize(" ");
 195   fOutput = new TList();
 196   fOutput->SetOwner();
 197   fSelTracks = new TObjArray;
 198   fSelPrimTracks = new TObjArray;
 199
 200   if (fDoNtuple) {
 201     TFile *f = OpenFile(1);
 202     if (f) {
 203       f->SetCompressionLevel(2);
 204       fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
 205       fNtuple->SetDirectory(f);
 206       if (fTrainMode) {
 207         fNtuple->SetAutoFlush(-2*1024*1024);
 208         fNtuple->SetAutoSave(0);
 209       } else {
 210         fNtuple->SetAutoFlush(-32*1024*1024);
 211         fNtuple->SetAutoSave(0);
 212       }
 213
 214       fHeader = new AliStaHeader;
 215       fNtuple->Branch("header", &fHeader, 16*1024, 99);
 216       fPrimVert = new AliStaVertex;
 217       fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
 218       fSpdVert = new AliStaVertex;
 219       fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
 220       fTpcVert = new AliStaVertex;
 221       fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
 222       if (TClass::GetClass("AliStaCluster"))
 223         TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
 224       fClusters = new TClonesArray("AliStaCluster");
 225       fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
 226       if (TClass::GetClass("AliStaTrigger"))
 227         TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
 228       fTriggers = new TClonesArray("AliStaTrigger");
 229       fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
 230     }
 231   }
 232
 233   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 234   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
 235   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
 236
 237   // histograms
 238   TH1::SetDefaultSumw2(kTRUE);
 239   TH2::SetDefaultSumw2(kTRUE);
 240   fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
 241   fHCuts->GetXaxis()->SetBinLabel(1,"All");
 242   fHCuts->GetXaxis()->SetBinLabel(2,"PS");
 243   fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 244   fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
 245   fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 246   fOutput->Add(fHCuts);
 247   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 248   fHVertexZ->SetXTitle("z [cm]");
 249   fOutput->Add(fHVertexZ);
 250   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 251   fHVertexZ2->SetXTitle("z [cm]");
 252   fOutput->Add(fHVertexZ2);
 253   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 254   fHCent->SetXTitle(fCentVar.Data());
 255   fOutput->Add(fHCent);
 256   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 257   fHCentQual->SetXTitle(fCentVar.Data());
 258   fOutput->Add(fHCentQual);
 259   fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 260   fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
 261   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 262     const char *name = fTrClassNamesArr->At(i)->GetName();
 263     fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
 264     fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
 265   }
 266   fOutput->Add(fHTclsBeforeCuts);
 267   fOutput->Add(fHTclsAfterCuts);
 268
 269   // histograms for cells
 270   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 271   fHColuRow   = new TH2*[nsm];
 272   fHColuRowE  = new TH2*[nsm];
 273   fHCellMult  = new TH1*[nsm];
 274   for (Int_t i = 0; i<nsm; ++i) {
 275     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
 276     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 277     fHColuRow[i]->SetXTitle("col (i#eta)");
 278     fHColuRow[i]->SetYTitle("row (i#phi)");
 279     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
 280     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 281     fHColuRowE[i]->SetXTitle("col (i#eta)");
 282     fHColuRowE[i]->SetYTitle("row (i#phi)");
 283     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 284     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 285     fHCellMult[i]->SetXTitle("# of cells");
 286     fOutput->Add(fHColuRow[i]);
 287     fOutput->Add(fHColuRowE[i]);
 288     fOutput->Add(fHCellMult[i]);
 289   }
 290   fHCellE = new TH1F("hCellE","",250,0.,25.);
 291   fHCellE->SetXTitle("E_{cell} [GeV]");
 292   fOutput->Add(fHCellE);
 293   fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
 294   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 295   fOutput->Add(fHCellH);
 296   fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
 297   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 298   fOutput->Add(fHCellM);
 299   fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
 300   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 301   fOutput->Add(fHCellM2);
 302
 303   fHCellFreqNoCut   = new TH1*[nsm];
 304   fHCellFreqCut100M = new TH1*[nsm];
 305   fHCellFreqCut300M = new TH1*[nsm];
 306   fHCellFreqE       = new TH1*[nsm];
 307   for (Int_t i = 0; i<nsm; ++i){
 308     Double_t lbin = i*24*48-0.5;
 309     Double_t hbin = lbin+24*48;
 310     fHCellFreqNoCut[i]   = new TH1F(Form("hCellFreqNoCut_SM%d",i),
 311                                     Form("Frequency SM%d (no cut);id;#",i),  1152, lbin, hbin);
 312     fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
 313                                     Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
 314     fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
 315                                     Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
 316     fHCellFreqE[i]       = new TH1F(Form("hCellFreqE_SM%d",i),
 317                                     Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
 318     fOutput->Add(fHCellFreqNoCut[i]);
 319     fOutput->Add(fHCellFreqCut100M[i]);
 320     fOutput->Add(fHCellFreqCut300M[i]);
 321     fOutput->Add(fHCellFreqE[i]);
 322   }
 323   if (!fMarkCells.IsNull()) {
 324     fHCellCheckE = new TH1*[24*48*nsm];
 325     memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
 326     TObjArray *cells = fMarkCells.Tokenize(" ");
 327     Int_t n = cells->GetEntries();
 328     Int_t *tcs = new Int_t[n];
 329     for (Int_t i=0;i<n;++i) {
 330       TString name(cells->At(i)->GetName());
 331       tcs[i]=name.Atoi();
 332     }
 333     for (Int_t i = 0; i<n; ++i) {
 334       Int_t c=tcs[i];
 335       if (c<24*48*nsm) {
 336         fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
 337         fOutput->Add(fHCellCheckE[i]);
 338       }
 339     }
 340     delete cells;
 341     delete [] tcs;
 342   }
 343
 344   // histograms for clusters
 345   if (!fTrainMode) {
 346     fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 347     fHClustEccentricity->SetXTitle("#epsilon_{C}");
 348     fOutput->Add(fHClustEccentricity);
 349     fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
 350     fHClustEtaPhi->SetXTitle("#eta");
 351     fHClustEtaPhi->SetYTitle("#varphi");
 352     fOutput->Add(fHClustEtaPhi);
 353     fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 354     fHClustEnergyPt->SetXTitle("E [GeV]");
 355     fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 356     fOutput->Add(fHClustEnergyPt);
 357     fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 358     fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 359     fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 360     fOutput->Add(fHClustEnergySigma);
 361     fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 362     fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 363     fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 364     fOutput->Add(fHClustSigmaSigma);
 365     fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 366     fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 367     fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 368     fOutput->Add(fHClustNCellEnergyRatio);
 369   }
 370
 371   // histograms for track matching
 372   fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
 373   fOutput->Add(fHMatchDr);
 374   fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
 375   fOutput->Add(fHMatchDz);
 376   fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
 377   fOutput->Add(fHMatchEp);
 378
 379   // histograms for pion candidates
 380   if (!fTrainMode) {
 381     fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
 382     fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 383     fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 384     fOutput->Add(fHPionEtaPhi);
 385     fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 386     fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 387     fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 388     fOutput->Add(fHPionMggPt);
 389     fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 390     fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 391     fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 392     fOutput->Add(fHPionMggAsym);
 393     fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 394     fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 395     fHPionMggDgg->SetYTitle("opening angle [grad]");
 396     fOutput->Add(fHPionMggDgg);
 397     const Int_t nbins = 20;
 398     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};
 399     fPtRanges = new TAxis(nbins-1,xbins);
 400     for (Int_t i = 0; i<=nbins; ++i) {
 401       fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 402       fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 403       if (i==0)
 404         fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 405       else if (i==nbins)
 406         fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 407       else
 408         fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 409       fOutput->Add(fHPionInvMasses[i]);
 410     }
 411   }
 412
 413   PostData(1, fOutput);
 414 }
 415
 416 //________________________________________________________________________
 417 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 418 {
 419   // Called for each event.
 420
 421   if (!InputEvent())
 422     return;
 423
 424   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 425   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 426   if (fEsdEv) {
 427     am->LoadBranch("AliESDRun.");
 428     am->LoadBranch("AliESDHeader.");
 429   } else {
 430     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 431     if (!fAodEv) {
 432       AliFatal("Neither ESD nor AOD event found");
 433       return;
 434     }
 435     am->LoadBranch("header");
 436   }
 437
 438   if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
 439     AliWarning("Accessing geometry from OCDB, this is not very efficient!");
 440     AliCDBManager *cdb = AliCDBManager::Instance();
 441     if (!cdb->IsDefaultStorageSet())
 442       cdb->SetDefaultStorage("raw://");
 443     Int_t runno = InputEvent()->GetRunNumber();
 444     if (runno != cdb->GetRun())
 445       cdb->SetRun(runno);
 446     AliGeomManager::LoadGeometry();
 447   }
 448
 449   if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
 450     Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
 451     if (fEsdEv) {
 452       for (Int_t i=0; i<nsm; ++i)
 453         fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 454     } else {
 455       for (Int_t i=0; i<nsm; ++i)
 456         fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 457     }
 458     fIsGeoMatsSet = kTRUE;
 459   }
 460
 461   if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 462     if (fEsdEv) {
 463       const AliESDRun *erun = fEsdEv->GetESDRun();
 464       AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
 465                                                erun->GetCurrentDip(),
 466                                                AliMagF::kConvLHC,
 467                                                kFALSE,
 468                                                erun->GetBeamEnergy(),
 469                                                erun->GetBeamType());
 470       TGeoGlobalMagField::Instance()->SetField(field);
 471     } else {
 472       Double_t pol = -1; //polarity
 473       Double_t be = -1;  //beam energy
 474       AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
 475       Int_t runno = fAodEv->GetRunNumber();
 476       if (runno>=136851 && runno<138275) {
 477         pol = -1;
 478         be = 2760;
 479         btype = AliMagF::kBeamTypeAA;
 480       } else if (runno>=138275 && runno<=139517) {
 481         pol = +1;
 482         be = 2760;
 483         btype = AliMagF::kBeamTypeAA;
 484       } else {
 485         AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
 486       }
 487       TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
 488     }
 489   }
 490
 491   Int_t cut = 1;
 492   fHCuts->Fill(cut++);
 493
 494   TString trgclasses;
 495   AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
 496   if (h) {
 497     trgclasses = fEsdEv->GetFiredTriggerClasses();
 498   } else {
 499     AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
 500     if (h2)
 501       trgclasses = h2->GetFiredTriggerClasses();
 502   }
 503   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 504     const char *name = fTrClassNamesArr->At(i)->GetName();
 505     if (trgclasses.Contains(name))
 506       fHTclsBeforeCuts->Fill(1+i);
 507   }
 508
 509   UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 510   if (res==0)
 511     return;
 512
 513   fHCuts->Fill(cut++);
 514
 515   const AliCentrality *centP = InputEvent()->GetCentrality();
 516   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 517   fHCent->Fill(cent);
 518   if (cent<fCentFrom||cent>fCentTo)
 519     return;
 520
 521   fHCuts->Fill(cut++);
 522
 523   if (fUseQualFlag) {
 524     if (centP->GetQuality()>0)
 525       return;
 526   }
 527
 528   fHCentQual->Fill(cent);
 529   fHCuts->Fill(cut++);
 530
 531   if (fEsdEv) {
 532     am->LoadBranch("PrimaryVertex.");
 533     am->LoadBranch("SPDVertex.");
 534     am->LoadBranch("TPCVertex.");
 535   } else {
 536     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 537     am->LoadBranch("vertices");
 538     if (!fAodEv) return;
 539   }
 540
 541   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 542   if (!vertex)
 543     return;
 544
 545   fHVertexZ->Fill(vertex->GetZ());
 546
 547   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 548     return;
 549
 550   fHCuts->Fill(cut++);
 551   fHVertexZ2->Fill(vertex->GetZ());
 552
 553   // count number of accepted events
 554   ++fNEvs;
 555
 556   for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
 557     const char *name = fTrClassNamesArr->At(i)->GetName();
 558     if (trgclasses.Contains(name))
 559       fHTclsAfterCuts->Fill(1+i);
 560   }
 561
 562   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 563   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
 564   fEsdCells    = 0; // will be set if ESD input used
 565   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
 566                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 567   fAodCells    = 0; // will be set if AOD input used
 568
 569   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 570   Bool_t clusattached = 0;
 571   Bool_t recalibrated = 0;
 572   if (1 && !fClusName.IsNull()) {
 573     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 574     TObjArray *ts = am->GetTasks();
 575     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 576     if (cltask && cltask->GetClusters()) {
 577       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 578       clusattached = cltask->GetAttachClusters();
 579       if (cltask->GetCalibData()!=0)
 580         recalibrated = kTRUE;
 581     }
 582   }
 583   if (1 && AODEvent() && !fClusName.IsNull()) {
 584     TList *l = AODEvent()->GetList();
 585     TClonesArray *clus = 0;
 586     if (l) {
 587       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 588       fAodClusters = clus;
 589     }
 590   }
 591
 592   if (fEsdEv) { // ESD input mode
 593     if (1 && (!fRecPoints||clusattached)) {
 594       if (!clusattached)
 595         am->LoadBranch("CaloClusters");
 596       TList *l = fEsdEv->GetList();
 597       TClonesArray *clus = 0;
 598       if (l) {
 599         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 600         fEsdClusters = clus;
 601       }
 602     }
 603     if (1) {
 604       if (!recalibrated)
 605         am->LoadBranch("EMCALCells.");
 606       fEsdCells = fEsdEv->GetEMCALCells();
 607     }
 608   } else if (fAodEv) { // AOD input mode
 609     if (1 && (!fAodClusters || clusattached)) {
 610       if (!clusattached)
 611         am->LoadBranch("caloClusters");
 612       TList *l = fAodEv->GetList();
 613       TClonesArray *clus = 0;
 614       if (l) {
 615         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 616         fAodClusters = clus;
 617       }
 618     }
 619     if (1) {
 620       if (!recalibrated)
 621         am->LoadBranch("emcalCells");
 622       fAodCells = fAodEv->GetEMCALCells();
 623     }
 624   } else {
 625     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 626   }
 627
 628   if (1) {
 629     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 630     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 631     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 632     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 633     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 634   }
 635
 636   if (fDoAfterburner)
 637     ClusterAfterburner();
 638
 639   CalcCaloTriggers();
 640   CalcPrimTracks();
 641   CalcTracks();
 642   CalcClusterProps();
 643
 644   FillCellHists();
 645   if (!fTrainMode) {
 646     FillClusHists();
 647     FillPionHists();
 648     FillOtherHists();
 649   }
 650   FillNtuple();
 651
 652   if (fTrainMode) {
 653     fSelTracks->Clear();
 654     fSelPrimTracks->Clear();
 655     fTriggers->Clear();
 656   }
 657
 658   PostData(1, fOutput);
 659 }
 660
 661 //________________________________________________________________________
 662 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 663 {
 664   // Terminate called at the end of analysis.
 665
 666   if (fNtuple) {
 667     TFile *f = OpenFile(1);
 668     if (f)
 669       fNtuple->Write();
 670   }
 671
 672   AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
 673 }
 674
 675 //________________________________________________________________________
 676 void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
 677 {
 678   // Calculate triggers
 679
 680   fTriggers->Clear();
 681
 682   AliVCaloCells *cells = fEsdCells;
 683   if (!cells)
 684     cells = fAodCells;
 685   if (!cells)
 686     return;
 687
 688   Int_t ncells = cells->GetNumberOfCells();
 689   if (ncells<=0)
 690     return;
 691
 692   if (ncells>fNAmpInTrigger) {
 693     delete [] fAmpInTrigger;
 694     fAmpInTrigger = new Float_t[ncells];
 695     fNAmpInTrigger = ncells;
 696   }
 697   for (Int_t i=0;i<ncells;++i)
 698     fAmpInTrigger[i] = 0;
 699
 700   std::map<Short_t,Short_t> map;
 701   for (Short_t pos=0;pos<ncells;++pos) {
 702     Short_t id = cells->GetCellNumber(pos);
 703     map[id]=pos;
 704   }
 705
 706   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 707   am->LoadBranch("EMCALTrigger.");
 708
 709   AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
 710   if (!triggers)
 711     return;
 712   if (triggers->GetEntries()<=0)
 713     return;
 714
 715   triggers->Reset();
 716   Int_t ntrigs=0;
 717   while (triggers->Next()) {
 718     Int_t gCol=0, gRow=0, ntimes=0;
 719     triggers->GetPosition(gCol,gRow);
 720     triggers->GetNL0Times(ntimes);
 721     if (ntimes<1)
 722       continue;
 723     Float_t amp=0;
 724     triggers->GetAmplitude(amp);
 725     Int_t find = -1;
 726     fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
 727     if (find<0)
 728       continue;
 729     Int_t cidx[4] = {-1};
 730     Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
 731     if (!ret)
 732       continue;
 733     Int_t trgtimes[25];
 734     triggers->GetL0Times(trgtimes);
 735     Int_t mintime = trgtimes[0];
 736     Int_t maxtime = trgtimes[0];
 737     Bool_t trigInTimeWindow = 0;
 738     for (Int_t i=0;i<ntimes;++i) {
 739       if (trgtimes[i]<mintime)
 740         mintime = trgtimes[i];
 741       if (maxtime<trgtimes[i])
 742         maxtime = trgtimes[i];
 743       if ((fMinL0Time<trgtimes[i]) && (fMaxL0Time>trgtimes[i]))
 744         trigInTimeWindow = 1;
 745     }
 746
 747     Double_t tenergy = 0;
 748     Double_t tphi=0;
 749     Double_t teta=0;
 750     for (Int_t i=0;i<3;++i) {
 751       Short_t pos = -1;
 752       std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
 753       if (it!=map.end())
 754         pos = it->second;
 755       if (pos<0)
 756         continue;
 757       if (trigInTimeWindow)
 758         fAmpInTrigger[pos] = amp; // save for usage in CalcClusProperties
 759       Float_t eta=-1, phi=-1;
 760       fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
 761       Double_t en= cells->GetAmplitude(pos);
 762       tenergy+=en;
 763       teta+=eta*en;
 764       tphi+=phi*en;
 765     }
 766     teta/=tenergy;
 767     tphi/=tenergy;
 768     AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
 769     trignew->fE       = tenergy;
 770     trignew->fEta     = teta;
 771     trignew->fPhi     = tphi;
 772     trignew->fAmp     = amp;
 773     trignew->fMinTime = mintime;
 774     trignew->fMaxTime = maxtime;
 775   }
 776 }
 777
 778 //________________________________________________________________________
 779 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
 780 {
 781   // Calculate cluster properties
 782
 783   fClusters->Clear();
 784
 785   AliVCaloCells *cells = fEsdCells;
 786   if (!cells)
 787     cells = fAodCells;
 788   if (!cells)
 789     return;
 790
 791   TObjArray *clusters = fEsdClusters;
 792   if (!clusters)
 793     clusters = fAodClusters;
 794   if (!clusters)
 795     return;
 796
 797   Int_t ncells = cells->GetNumberOfCells();
 798   Int_t nclus  = clusters->GetEntries();
 799   Int_t ntrks  = fSelTracks->GetEntries();
 800   Bool_t btracks[6][ntrks];
 801   memset(btracks,0,sizeof(btracks));
 802
 803   std::map<Short_t,Short_t> map;
 804   for (Short_t pos=0;pos<ncells;++pos) {
 805     Short_t id = cells->GetCellNumber(pos);
 806     map[id]=pos;
 807   }
 808
 809   for(Int_t i=0, ncl=0; i<nclus; ++i) {
 810     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 811     if (!clus)
 812       continue;
 813     if (!clus->IsEMCAL())
 814       continue;
 815     if (clus->E()<fMinE)
 816       continue;
 817
 818     Float_t clsPos[3] = {0};
 819     clus->GetPosition(clsPos);
 820     TVector3 clsVec(clsPos);
 821     Double_t vertex[3] = {0};
 822     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 823     TLorentzVector clusterVec;
 824     clus->GetMomentum(clusterVec,vertex);
 825     Double_t clsEta = clusterVec.Eta();
 826
 827     AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
 828     cl->fE        = clus->E();
 829     cl->fR        = clsVec.Perp();
 830     cl->fEta      = clsVec.Eta();
 831     cl->fPhi      = clsVec.Phi();
 832     cl->fN        = clus->GetNCells();
 833     cl->fN1       = GetNCells(clus,0.100);
 834     cl->fN3       = GetNCells(clus,0.300);
 835     Short_t id    = -1;
 836     Double_t emax = GetMaxCellEnergy(clus, id);
 837     cl->fIdMax    = id;
 838     cl->fEmax     = emax;
 839     if (clus->GetDistanceToBadChannel()<10000)
 840       cl->fDbc    = clus->GetDistanceToBadChannel();
 841     if (!TMath::IsNaN(clus->GetDispersion()))
 842       cl->fDisp   = clus->GetDispersion();
 843     if (!TMath::IsNaN(clus->GetM20()))
 844       cl->fM20      = clus->GetM20();
 845     if (!TMath::IsNaN(clus->GetM02()))
 846       cl->fM02      = clus->GetM02();
 847     Double_t maxAxis = 0, minAxis = 0;
 848     GetSigma(clus,maxAxis,minAxis);
 849     clus->SetTOF(maxAxis);     // store sigma in TOF
 850     cl->fSig      = maxAxis;
 851     Double_t clusterEcc = 0;
 852     if (maxAxis > 0)
 853       clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 854     clus->SetChi2(clusterEcc); // store ecc in chi2
 855     cl->fEcc      = clusterEcc;
 856     cl->fTrIso    = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
 857     cl->fTrIso1   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
 858     cl->fTrIso2   = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
 859     cl->fCeCore   = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
 860     cl->fCeIso    = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
 861     Double_t trigpen = 0;
 862     Double_t trignen = 0;
 863     for(Int_t j=0; j<cl->fN; ++j) {
 864       Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
 865       Short_t pos = -1;
 866       std::map<Short_t,Short_t>::iterator it = map.find(cid);
 867       if (it!=map.end())
 868         pos = it->second;
 869       if (pos<0)
 870         continue;
 871       if (fAmpInTrigger[pos]>0)
 872         trigpen += cells->GetAmplitude(pos);
 873       else if (fAmpInTrigger[pos]<0)
 874         trignen += cells->GetAmplitude(pos);
 875     }
 876     if (trigpen>0) {
 877       cl->fIsTrigM = 1;
 878       cl->fTrigE   = trigpen;
 879     }
 880     if (trignen>0) {
 881       cl->fIsTrigM   = 1;
 882       cl->fTrigMaskE = trignen;
 883     }
 884
 885     // track matching
 886     Double_t mind2 = 1e10;
 887     for(Int_t j = 0; j<ntrks; ++j) {
 888       AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
 889       if (!track)
 890         continue;
 891
 892       if (TMath::Abs(clsEta-track->Eta())>0.5)
 893         continue;
 894
 895       TVector3 vec(clsPos);
 896       Int_t index =  (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
 897       if (btracks[index-4][j]) {
 898         continue;
 899       }
 900
 901       Float_t tmpR=-1, tmpZ=-1;
 902       if (!fDoTrMatGeom) {
 903         AliExternalTrackParam *tParam = 0;
 904         if (fEsdEv) {
 905           AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
 906           tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
 907         } else
 908           tParam = new AliExternalTrackParam(track);
 909
 910         Double_t bfield[3] = {0};
 911         track->GetBxByBz(bfield);
 912         Double_t alpha = (index+0.5)*20*TMath::DegToRad();
 913         vec.RotateZ(-alpha);   //Rotate the cluster to the local extrapolation coordinate system
 914         tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
 915         Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
 916         if (!ret) {
 917           btracks[index-4][j]=1;
 918           delete tParam;
 919           continue;
 920         }
 921         Double_t trkPos[3] = {0};
 922         tParam->GetXYZ(trkPos); //Get the extrapolated global position
 923         tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
 924         tmpZ = clsPos[2]-trkPos[2];
 925         delete tParam;
 926       } else {
 927         if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
 928           continue;
 929         AliExternalTrackParam tParam(track);
 930         if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
 931           continue;
 932       }
 933
 934       Double_t d2 = tmpR;
 935       if (mind2>d2) {
 936         mind2=d2;
 937         cl->fTrDz   = tmpZ;
 938         cl->fTrDr   = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
 939         cl->fTrEp   = clus->E()/track->P();
 940         cl->fIsTrackM = 1;
 941       }
 942     }
 943
 944     if (cl->fIsTrackM) {
 945       fHMatchDr->Fill(cl->fTrDr);
 946       fHMatchDz->Fill(cl->fTrDz);
 947       fHMatchEp->Fill(cl->fTrEp);
 948     }
 949   }
 950 }
 951
 952 //________________________________________________________________________
 953 void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
 954 {
 955   // Calculate track properties.
 956
 957   fSelPrimTracks->Clear();
 958
 959   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 960   TClonesArray *tracks = 0;
 961   if (fEsdEv) {
 962     am->LoadBranch("Tracks");
 963     TList *l = fEsdEv->GetList();
 964     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
 965   } else {
 966     am->LoadBranch("tracks");
 967     TList *l = fAodEv->GetList();
 968     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
 969   }
 970
 971   if (!tracks)
 972     return;
 973
 974   AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
 975   Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
 976   Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
 977
 978   if (fEsdEv) {
 979     fSelPrimTracks->SetOwner(kTRUE);
 980     am->LoadBranch("PrimaryVertex.");
 981     const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
 982     am->LoadBranch("SPDVertex.");
 983     const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
 984     am->LoadBranch("Tracks");
 985     const Int_t Ntracks = tracks->GetEntries();
 986     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 987       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
 988       if (!track) {
 989         AliWarning(Form("Could not receive track %d\n", iTracks));
 990         continue;
 991       }
 992       if (fTrCuts && !fTrCuts->IsSelected(track))
 993         continue;
 994       Double_t eta = track->Eta();
 995       if (eta<-1||eta>1)
 996         continue;
 997       if (track->Phi()<phimin||track->Phi()>phimax)
 998         continue;
 999
1000       AliESDtrack copyt(*track);
1001       Double_t bfield[3];
1002       copyt.GetBxByBz(bfield);
1003       AliExternalTrackParam tParam;
1004       Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1005       if (!relate)
1006         continue;
1007       copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1008
1009       Double_t p[3]      = { 0. };
1010       copyt.GetPxPyPz(p);
1011       Double_t pos[3]    = { 0. };
1012       copyt.GetXYZ(pos);
1013       Double_t covTr[21] = { 0. };
1014       copyt.GetCovarianceXYZPxPyPz(covTr);
1015       Double_t pid[10]   = { 0. };
1016       copyt.GetESDpid(pid);
1017       AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1018                                             copyt.GetLabel(),
1019                                             p,
1020                                             kTRUE,
1021                                             pos,
1022                                             kFALSE,
1023                                             covTr,
1024                                             (Short_t)copyt.GetSign(),
1025                                             copyt.GetITSClusterMap(),
1026                                             pid,
1027                                             0,/*fPrimaryVertex,*/
1028                                             kTRUE, // check if this is right
1029                                             vtx->UsesTrack(copyt.GetID()));
1030       aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1031       aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1032       Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1033       if(ndf>0)
1034         aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1035       else
1036         aTrack->SetChi2perNDF(-1);
1037       aTrack->SetFlags(copyt.GetStatus());
1038       aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
1039       fSelPrimTracks->Add(aTrack);
1040     }
1041   } else {
1042     Int_t ntracks = tracks->GetEntries();
1043     for (Int_t i=0; i<ntracks; ++i) {
1044       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1045       if (!track)
1046         continue;
1047       Double_t eta = track->Eta();
1048       if (eta<-1||eta>1)
1049         continue;
1050       if (track->Phi()<phimin||track->Phi()>phimax)
1051         continue;
1052       if(track->GetTPCNcls()<fMinNClusPerTr)
1053         continue;
1054       //todo: Learn how to set/filter AODs for prim/sec tracks
1055       fSelPrimTracks->Add(track);
1056     }
1057   }
1058 }
1059
1060 //________________________________________________________________________
1061 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
1062 {
1063   // Calculate track properties (including secondaries).
1064
1065   fSelTracks->Clear();
1066
1067   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1068   TClonesArray *tracks = 0;
1069   if (fEsdEv) {
1070     am->LoadBranch("Tracks");
1071     TList *l = fEsdEv->GetList();
1072     tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1073   } else {
1074     am->LoadBranch("tracks");
1075     TList *l = fAodEv->GetList();
1076     tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1077   }
1078
1079   if (!tracks)
1080     return;
1081
1082   if (fEsdEv) {
1083     const Int_t Ntracks = tracks->GetEntries();
1084     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1085       AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1086       if (!track) {
1087         AliWarning(Form("Could not receive track %d\n", iTracks));
1088         continue;
1089       }
1090       if (fTrCuts && !fTrCuts->IsSelected(track))
1091         continue;
1092       Double_t eta = track->Eta();
1093       if (eta<-1||eta>1)
1094         continue;
1095       fSelTracks->Add(track);
1096     }
1097   } else {
1098     Int_t ntracks = tracks->GetEntries();
1099     for (Int_t i=0; i<ntracks; ++i) {
1100       AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1101       if (!track)
1102         continue;
1103       Double_t eta = track->Eta();
1104       if (eta<-1||eta>1)
1105         continue;
1106       if(track->GetTPCNcls()<fMinNClusPerTr)
1107         continue;
1108
1109       if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
1110         AliExternalTrackParam tParam(track);
1111         if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1112           Double_t trkPos[3];
1113           tParam.GetXYZ(trkPos);
1114           track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1115         }
1116       }
1117       fSelTracks->Add(track);
1118     }
1119   }
1120 }
1121
1122 //________________________________________________________________________
1123 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
1124 {
1125   // Run custer reconstruction afterburner.
1126
1127   AliVCaloCells *cells = fEsdCells;
1128   if (!cells)
1129     cells = fAodCells;
1130
1131   if (!cells)
1132     return;
1133
1134   Int_t ncells = cells->GetNumberOfCells();
1135   if (ncells<=0)
1136     return;
1137
1138   Double_t cellMeanE = 0, cellSigE = 0;
1139   for (Int_t i = 0; i<ncells; ++i) {
1140     Double_t cellE = cells->GetAmplitude(i);
1141     cellMeanE += cellE;
1142     cellSigE += cellE*cellE;
1143   }
1144   cellMeanE /= ncells;
1145   cellSigE /= ncells;
1146   cellSigE -= cellMeanE*cellMeanE;
1147   if (cellSigE<0)
1148     cellSigE = 0;
1149   cellSigE = TMath::Sqrt(cellSigE / ncells);
1150
1151   Double_t subE = cellMeanE - 7*cellSigE;
1152   if (subE<0)
1153     return;
1154
1155   for (Int_t i = 0; i<ncells; ++i) {
1156     Short_t id=-1;
1157     Double_t amp=0,time=0;
1158     if (!cells->GetCell(i, id, amp, time))
1159       continue;
1160     amp -= cellMeanE;
1161     if (amp<0.001)
1162       amp = 0;
1163     cells->SetCell(i, id, amp, time);
1164   }
1165
1166   TObjArray *clusters = fEsdClusters;
1167   if (!clusters)
1168     clusters = fAodClusters;
1169   if (!clusters)
1170     return;
1171
1172   Int_t nclus = clusters->GetEntries();
1173   for (Int_t i = 0; i<nclus; ++i) {
1174     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1175     if (!clus->IsEMCAL())
1176       continue;
1177     Int_t nc = clus->GetNCells();
1178     Double_t clusE = 0;
1179     UShort_t ids[100] = {0};
1180     Double_t fra[100] = {0};
1181     for (Int_t j = 0; j<nc; ++j) {
1182       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1183       Double_t cen = cells->GetCellAmplitude(id);
1184       clusE += cen;
1185       if (cen>0) {
1186         ids[nc] = id;
1187         ++nc;
1188       }
1189     }
1190     if (clusE<=0) {
1191       clusters->RemoveAt(i);
1192       continue;
1193     }
1194
1195     for (Int_t j = 0; j<nc; ++j) {
1196       Short_t id = ids[j];
1197       Double_t cen = cells->GetCellAmplitude(id);
1198       fra[j] = cen/clusE;
1199     }
1200     clus->SetE(clusE);
1201     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1202     if (aodclus) {
1203       aodclus->Clear("");
1204       aodclus->SetNCells(nc);
1205       aodclus->SetCellsAmplitudeFraction(fra);
1206       aodclus->SetCellsAbsId(ids);
1207     }
1208   }
1209   clusters->Compress();
1210 }
1211
1212 //________________________________________________________________________
1213 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1214 {
1215   // Fill histograms related to cell properties.
1216
1217   AliVCaloCells *cells = fEsdCells;
1218   if (!cells)
1219     cells = fAodCells;
1220
1221   if (!cells)
1222     return;
1223
1224   Int_t cellModCount[12] = {0};
1225   Double_t cellMaxE = 0;
1226   Double_t cellMeanE = 0;
1227   Int_t ncells = cells->GetNumberOfCells();
1228   if (ncells==0)
1229     return;
1230
1231   Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1232
1233   for (Int_t i = 0; i<ncells; ++i) {
1234     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1235     Double_t cellE = cells->GetAmplitude(i);
1236     fHCellE->Fill(cellE);
1237     if (cellE>cellMaxE)
1238       cellMaxE = cellE;
1239     cellMeanE += cellE;
1240
1241     Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1242     Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1243     if (!ret) {
1244       AliError(Form("Could not get cell index for %d", absID));
1245       continue;
1246     }
1247     ++cellModCount[iSM];
1248     Int_t iPhi=-1, iEta=-1;
1249     fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1250     fHColuRow[iSM]->Fill(iEta,iPhi,1);
1251     fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1252     fHCellFreqNoCut[iSM]->Fill(absID);
1253     if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1254     if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
1255     if (fHCellCheckE && fHCellCheckE[absID])
1256       fHCellCheckE[absID]->Fill(cellE);
1257     fHCellFreqE[iSM]->Fill(absID, cellE);
1258   }
1259   fHCellH->Fill(cellMaxE);
1260   cellMeanE /= ncells;
1261   fHCellM->Fill(cellMeanE);
1262   fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1263   for (Int_t i=0; i<nsm; ++i)
1264     fHCellMult[i]->Fill(cellModCount[i]);
1265 }
1266
1267 //________________________________________________________________________
1268 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1269 {
1270   // Fill histograms related to cluster properties.
1271
1272   TObjArray *clusters = fEsdClusters;
1273   if (!clusters)
1274     clusters = fAodClusters;
1275   if (!clusters)
1276     return;
1277
1278   Double_t vertex[3] = {0};
1279   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1280
1281   Int_t nclus = clusters->GetEntries();
1282   for(Int_t i = 0; i<nclus; ++i) {
1283     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1284     if (!clus)
1285       continue;
1286     if (!clus->IsEMCAL())
1287       continue;
1288     TLorentzVector clusterVec;
1289     clus->GetMomentum(clusterVec,vertex);
1290     Double_t maxAxis    = clus->GetTOF(); //sigma
1291     Double_t clusterEcc = clus->Chi2();   //eccentricity
1292     fHClustEccentricity->Fill(clusterEcc);
1293     fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1294     fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1295     fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1296     fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1297     fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1298   }
1299 }
1300
1301 //________________________________________________________________________
1302 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1303 {
1304   // Fill ntuple.
1305
1306   if (!fNtuple)
1307     return;
1308
1309   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1310   if (fAodEv) {
1311     fHeader->fRun            = fAodEv->GetRunNumber();
1312     fHeader->fOrbit          = fAodEv->GetHeader()->GetOrbitNumber();
1313     fHeader->fPeriod         = fAodEv->GetHeader()->GetPeriodNumber();
1314     fHeader->fBx             = fAodEv->GetHeader()->GetBunchCrossNumber();
1315     fHeader->fL0             = fAodEv->GetHeader()->GetL0TriggerInputs();
1316     fHeader->fL1             = fAodEv->GetHeader()->GetL1TriggerInputs();
1317     fHeader->fL2             = fAodEv->GetHeader()->GetL2TriggerInputs();
1318     fHeader->fTrClassMask    = fAodEv->GetHeader()->GetTriggerMask();
1319     fHeader->fTrCluster      = fAodEv->GetHeader()->GetTriggerCluster();
1320     fHeader->fOffTriggers    = fAodEv->GetHeader()->GetOfflineTrigger();
1321     fHeader->fFiredTriggers  = fAodEv->GetHeader()->GetFiredTriggerClasses();
1322   } else {
1323     fHeader->fRun            = fEsdEv->GetRunNumber();
1324     fHeader->fOrbit          = fEsdEv->GetHeader()->GetOrbitNumber();
1325     fHeader->fPeriod         = fEsdEv->GetESDRun()->GetPeriodNumber();
1326     fHeader->fBx             = fEsdEv->GetHeader()->GetBunchCrossNumber();
1327     fHeader->fL0             = fEsdEv->GetHeader()->GetL0TriggerInputs();
1328     fHeader->fL1             = fEsdEv->GetHeader()->GetL1TriggerInputs();
1329     fHeader->fL2             = fEsdEv->GetHeader()->GetL2TriggerInputs();
1330     fHeader->fTrClassMask    = fEsdEv->GetHeader()->GetTriggerMask();
1331     fHeader->fTrCluster      = fEsdEv->GetHeader()->GetTriggerCluster();
1332     fHeader->fOffTriggers    = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1333     fHeader->fFiredTriggers  = fEsdEv->GetFiredTriggerClasses();
1334   }
1335   AliCentrality *cent = InputEvent()->GetCentrality();
1336   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
1337   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
1338   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
1339   fHeader->fCqual     = cent->GetQuality();
1340
1341   AliEventplane *ep = InputEvent()->GetEventplane();
1342   if (ep) {
1343     if (ep->GetQVector())
1344       fHeader->fPsi     = ep->GetQVector()->Phi()/2. ;
1345     fHeader->fPsi = -1;
1346     if (ep->GetQsub1()&&ep->GetQsub2())
1347       fHeader->fPsiRes  = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1348     else
1349       fHeader->fPsiRes = 0;
1350   }
1351
1352   Double_t val = 0;
1353   TString trgclasses(fHeader->fFiredTriggers);
1354   for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1355     const char *name = fTrClassNamesArr->At(j)->GetName();
1356     if (trgclasses.Contains(name))
1357       val += TMath::Power(2,j);
1358   }
1359   fHeader->fTcls = (UInt_t)val;
1360
1361   if (fAodEv) {
1362     am->LoadBranch("vertices");
1363     AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1364     FillVertex(fPrimVert, pv);
1365     AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1366     FillVertex(fSpdVert, sv);
1367   } else {
1368     am->LoadBranch("PrimaryVertex.");
1369     const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1370     FillVertex(fPrimVert, pv);
1371     am->LoadBranch("SPDVertex.");
1372     const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1373     FillVertex(fSpdVert, sv);
1374     am->LoadBranch("TPCVertex.");
1375     const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1376     FillVertex(fTpcVert, tv);
1377   }
1378
1379   fNtuple->Fill();
1380 }
1381
1382 //________________________________________________________________________
1383 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1384 {
1385   // Fill histograms related to pions.
1386
1387   Double_t vertex[3] = {0};
1388   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1389
1390   TObjArray *clusters = fEsdClusters;
1391   if (!clusters)
1392     clusters = fAodClusters;
1393
1394   if (clusters) {
1395     TLorentzVector clusterVec1;
1396     TLorentzVector clusterVec2;
1397     TLorentzVector pionVec;
1398
1399     Int_t nclus = clusters->GetEntries();
1400     for (Int_t i = 0; i<nclus; ++i) {
1401       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1402       if (!clus1)
1403         continue;
1404       if (!clus1->IsEMCAL())
1405         continue;
1406       if (clus1->E()<fMinE)
1407         continue;
1408       if (clus1->GetNCells()<fNminCells)
1409         continue;
1410       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1411         continue;
1412       if (clus1->Chi2()<fMinEcc) // eccentricity cut
1413         continue;
1414       clus1->GetMomentum(clusterVec1,vertex);
1415       for (Int_t j = i+1; j<nclus; ++j) {
1416         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1417         if (!clus2)
1418           continue;
1419         if (!clus2->IsEMCAL())
1420           continue;
1421         if (clus2->E()<fMinE)
1422           continue;
1423         if (clus2->GetNCells()<fNminCells)
1424           continue;
1425         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1426           continue;
1427         if (clus2->Chi2()<fMinEcc) // eccentricity cut
1428           continue;
1429         clus2->GetMomentum(clusterVec2,vertex);
1430         pionVec = clusterVec1 + clusterVec2;
1431         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1432         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1433         if (pionZgg < fAsymMax) {
1434           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1435           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1436           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1437           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1438           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1439           fHPionInvMasses[bin]->Fill(pionVec.M());
1440         }
1441       }
1442     }
1443   }
1444 }
1445
1446 //________________________________________________________________________
1447 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1448 {
1449   // Fill other histograms.
1450 }
1451
1452 //__________________________________________________________________________________________________
1453 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1454 {
1455   // Fill vertex from ESD vertex info.
1456
1457   v->fVx   = esdv->GetX();
1458   v->fVy   = esdv->GetY();
1459   v->fVz   = esdv->GetZ();
1460   v->fVc   = esdv->GetNContributors();
1461   v->fDisp = esdv->GetDispersion();
1462   v->fZres = esdv->GetZRes();
1463   v->fChi2 = esdv->GetChi2();
1464   v->fSt   = esdv->GetStatus();
1465   v->fIs3D = esdv->IsFromVertexer3D();
1466   v->fIsZ  = esdv->IsFromVertexerZ();
1467 }
1468
1469 //__________________________________________________________________________________________________
1470 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1471 {
1472   // Fill vertex from AOD vertex info.
1473
1474   v->fVx   = aodv->GetX();
1475   v->fVy   = aodv->GetY();
1476   v->fVz   = aodv->GetZ();
1477   v->fVc   = aodv->GetNContributors();
1478   v->fChi2 = aodv->GetChi2();
1479 }
1480
1481 //________________________________________________________________________
1482 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1483 {
1484   // Compute isolation based on cell content.
1485
1486   AliVCaloCells *cells = fEsdCells;
1487   if (!cells)
1488     cells = fAodCells;
1489   if (!cells)
1490     return 0;
1491
1492   Double_t cellIsolation = 0;
1493   Double_t rad2 = radius*radius;
1494   Int_t ncells = cells->GetNumberOfCells();
1495   for (Int_t i = 0; i<ncells; ++i) {
1496     Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
1497     Double_t cellE = cells->GetAmplitude(i);
1498     Float_t eta, phi;
1499     fGeom->EtaPhiFromIndex(absID,eta,phi);
1500     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1501     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1502     if(dist>rad2)
1503       continue;
1504     cellIsolation += cellE;
1505   }
1506   return cellIsolation;
1507 }
1508
1509 //________________________________________________________________________
1510 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster, Short_t &id) const
1511 {
1512   // Get maximum energy of attached cell.
1513
1514   id = -1;
1515   Double_t maxe = 0;
1516   Int_t ncells = cluster->GetNCells();
1517   if (fEsdCells) {
1518     for (Int_t i=0; i<ncells; i++) {
1519       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1520       if (e>maxe) {
1521         maxe = e;
1522         id   = cluster->GetCellAbsId(i);
1523       }
1524     }
1525   } else {
1526     for (Int_t i=0; i<ncells; i++) {
1527       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1528       if (e>maxe)
1529         maxe = e;
1530         id   = cluster->GetCellAbsId(i);
1531     }
1532   }
1533   return maxe;
1534 }
1535
1536 //________________________________________________________________________
1537 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
1538 {
1539   // Calculate the (E) weighted variance along the longer (eigen) axis.
1540
1541   sigmaMax = 0;          // cluster variance along its longer axis
1542   sigmaMin = 0;          // cluster variance along its shorter axis
1543   Double_t Ec  = c->E(); // cluster energy
1544   if(Ec<=0)
1545     return;
1546   Double_t Xc  = 0 ;     // cluster first moment along X
1547   Double_t Yc  = 0 ;     // cluster first moment along Y
1548   Double_t Sxx = 0 ;     // cluster second central moment along X (variance_X^2)
1549   Double_t Sxy = 0 ;     // cluster second central moment along Y (variance_Y^2)
1550   Double_t Syy = 0 ;     // cluster covariance^2
1551
1552   AliVCaloCells *cells = fEsdCells;
1553   if (!cells)
1554     cells = fAodCells;
1555
1556   if (!cells)
1557     return;
1558
1559   Int_t ncells = c->GetNCells();
1560   if (ncells==1)
1561     return;
1562
1563   for(Int_t j=0; j<ncells; ++j) {
1564     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1565     Double_t cellen = cells->GetCellAmplitude(id);
1566     TVector3 pos;
1567     fGeom->GetGlobal(id,pos);
1568     Xc  += cellen*pos.X();
1569     Yc  += cellen*pos.Y();
1570     Sxx += cellen*pos.X()*pos.X();
1571     Syy += cellen*pos.Y()*pos.Y();
1572     Sxy += cellen*pos.X()*pos.Y();
1573   }
1574   Xc  /= Ec;
1575   Yc  /= Ec;
1576   Sxx /= Ec;
1577   Syy /= Ec;
1578   Sxy /= Ec;
1579   Sxx -= Xc*Xc;
1580   Syy -= Yc*Yc;
1581   Sxy -= Xc*Yc;
1582   Sxx = TMath::Abs(Sxx);
1583   Syy = TMath::Abs(Syy);
1584   sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1585   sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
1586   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
1587   sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
1588 }
1589
1590 //________________________________________________________________________
1591 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
1592 {
1593   // Calculate number of attached cells above emin.
1594
1595   AliVCaloCells *cells = fEsdCells;
1596   if (!cells)
1597     cells = fAodCells;
1598   if (!cells)
1599     return 0;
1600
1601   Int_t n = 0;
1602   Int_t ncells = c->GetNCells();
1603   for(Int_t j=0; j<ncells; ++j) {
1604     Int_t id = TMath::Abs(c->GetCellAbsId(j));
1605     Double_t cellen = cells->GetCellAmplitude(id);
1606     if (cellen>=emin)
1607       ++n;
1608   }
1609   return n;
1610 }
1611
1612 //________________________________________________________________________
1613 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1614 {
1615   // Compute isolation based on tracks.
1616
1617   Double_t trkIsolation = 0;
1618   Double_t rad2 = radius*radius;
1619   Int_t ntrks = fSelPrimTracks->GetEntries();
1620   for(Int_t j = 0; j<ntrks; ++j) {
1621     AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1622     if (!track)
1623       continue;
1624     if (track->Pt()<pt)
1625       continue;
1626     Float_t eta = track->Eta();
1627     Float_t phi = track->Phi();
1628     Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
1629     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1630     if(dist>rad2)
1631       continue;
1632     trkIsolation += track->Pt();
1633   }
1634   return trkIsolation;
1635 }
1636
1637