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 <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12 #include <TNtuple.h>
  13 #include "AliAODEvent.h"
  14 #include "AliAODVertex.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  17 #include "AliCentrality.h"
  18 #include "AliEMCALGeoUtils.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDVertex.h"
  21 #include "AliLog.h"
  22
  23 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
  27   : AliAnalysisTaskSE(name),
  28     fCentVar("V0M"),
  29     fCentFrom(0),
  30     fCentTo(100),
  31     fVtxZMin(-10),
  32     fVtxZMax(+10),
  33     fUseQualFlag(1),
  34     fClusName(),
  35     fDoNtuple(0),
  36     fDoAfterburner(0),
  37     fAsymMax(1),
  38     fNminCells(1),
  39     fMinE(0.100),
  40     fMinErat(0),
  41     fMinEcc(0),
  42     fNEvs(0),
  43     fGeom(0),
  44     fOutput(0),
  45     fEsdEv(0),
  46     fAodEv(0),
  47     fRecPoints(0),
  48     fEsdClusters(0),
  49     fEsdCells(0),
  50     fAodClusters(0),
  51     fAodCells(0),
  52     fPtRanges(0),
  53     fNtuple(0),
  54     fHCuts(0x0),
  55     fHVertexZ(0x0),
  56     fHVertexZ2(0x0),
  57     fHCent(0x0),
  58     fHCentQual(0x0),
  59     fHColuRow(0x0),
  60     fHColuRowE(0x0),
  61     fHCellMult(0x0),
  62     fHCellE(0x0),
  63     fHCellH(0x0),
  64     fHCellM(0x0),
  65     fHCellM2(0x0),
  66     fHClustEccentricity(0),
  67     fHClustEtaPhi(0x0),
  68     fHClustEnergyPt(0x0),
  69     fHClustEnergySigma(0x0),
  70     fHClustSigmaSigma(0x0),
  71     fHClustNCellEnergyRatio(0x0),
  72     fHPionEtaPhi(0x0),
  73     fHPionMggPt(0x0),
  74     fHPionMggAsym(0x0),
  75     fHPionMggDgg(0x0)
  76 {
  77   // Constructor.
  78
  79   if (!name)
  80     return;
  81   SetName(name);
  82   DefineInput(0, TChain::Class());
  83   DefineOutput(1, TList::Class());
  84   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells. "
  85                "AOD:header,vertices,emcalCells";
  86 }
  87
  88 //________________________________________________________________________
  89 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
  90 {
  91   // Destructor.
  92
  93   delete fOutput; fOutput = 0;
  94   delete fPtRanges; fPtRanges = 0;
  95   delete fGeom; fGeom = 0;
  96   delete [] fHColuRow;
  97   delete [] fHColuRowE;
  98   delete [] fHCellMult;
  99 }
 100
 101 //________________________________________________________________________
 102 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
 103 {
 104   // Create user objects here.
 105
 106   fOutput = new TList();
 107   fOutput->SetOwner();
 108
 109   if (fDoNtuple) {
 110     TFile *f = OpenFile(1);
 111     if (f)
 112       fNtuple = new TNtuple(Form("nt%.0fto%.0f",fCentFrom,fCentTo),"nt",
 113                             "run:evt:cent:pt:e:emax:n:db:disp:mn:ms:chi:cpv:ecc:sig:eta:phi");
 114   }
 115
 116   // histograms
 117   TH1::SetDefaultSumw2(kTRUE);
 118   TH2::SetDefaultSumw2(kTRUE);
 119   fHCuts = new TH1F("hEventCuts","",4,0.5,4.5);
 120   fHCuts->GetXaxis()->SetBinLabel(1,"All (PS)");
 121   fHCuts->GetXaxis()->SetBinLabel(2,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
 122   fHCuts->GetXaxis()->SetBinLabel(3,"QFlag");
 123   fHCuts->GetXaxis()->SetBinLabel(4,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
 124   fOutput->Add(fHCuts);
 125   fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
 126   fHVertexZ->SetXTitle("z [cm]");
 127   fOutput->Add(fHVertexZ);
 128   fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
 129   fHVertexZ2->SetXTitle("z [cm]");
 130   fOutput->Add(fHVertexZ2);
 131   fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
 132   fHCent->SetXTitle(fCentVar.Data());
 133   fOutput->Add(fHCent);
 134   fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
 135   fHCentQual->SetXTitle(fCentVar.Data());
 136   fOutput->Add(fHCentQual);
 137
 138   // histograms for cells
 139   fHColuRow   = new TH2F*[4];
 140   fHColuRowE  = new TH2F*[4];
 141   fHCellMult  = new TH1F*[4];
 142   for (Int_t i = 0; i<4; ++i) {
 143     fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",49,0,49,25,0,25);
 144     fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
 145     fHColuRow[i]->SetXTitle("col (i#eta)");
 146     fHColuRow[i]->SetYTitle("row (i#phi)");
 147     fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",49,0,49,25,0,25);
 148     fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
 149     fHColuRowE[i]->SetXTitle("col (i#eta)");
 150     fHColuRowE[i]->SetYTitle("row (i#phi)");
 151     fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
 152     fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
 153     fHCellMult[i]->SetXTitle("# of cells");
 154     fOutput->Add(fHColuRow[i]);
 155     fOutput->Add(fHColuRowE[i]);
 156     fOutput->Add(fHCellMult[i]);
 157   }
 158   fHCellE = new TH1F("hCellE","",150,0.,15.);
 159   fHCellE->SetXTitle("E_{cell} [GeV]");
 160   fOutput->Add(fHCellE);
 161   fHCellH = new TH1F ("fHCellHighestE","",150,0.,15.);
 162   fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
 163   fOutput->Add(fHCellH);
 164   fHCellM = new TH1F ("fHCellMeanEperHitCell","",250,0.,2.5);
 165   fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
 166   fOutput->Add(fHCellM);
 167   fHCellM2 = new TH1F ("fHCellMeanEperAllCells","",250,0.,1);
 168   fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
 169   fOutput->Add(fHCellM2);
 170
 171   // histograms for clusters
 172   fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
 173   fHClustEccentricity->SetXTitle("#epsilon_{C}");
 174   fOutput->Add(fHClustEccentricity);
 175   fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,500,1.2,2.2);
 176   fHClustEtaPhi->SetXTitle("#eta");
 177   fHClustEtaPhi->SetYTitle("#varphi");
 178   fOutput->Add(fHClustEtaPhi);
 179   fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
 180   fHClustEnergyPt->SetXTitle("E [GeV]");
 181   fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
 182   fOutput->Add(fHClustEnergyPt);
 183   fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
 184   fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
 185   fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
 186   fOutput->Add(fHClustEnergySigma);
 187   fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
 188   fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
 189   fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
 190   fOutput->Add(fHClustSigmaSigma);
 191   fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
 192   fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
 193   fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
 194   fOutput->Add(fHClustNCellEnergyRatio);
 195
 196   // histograms for pion candidates
 197   fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100,1.2,2.2);
 198   fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
 199   fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
 200   fOutput->Add(fHPionEtaPhi);
 201   fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
 202   fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 203   fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
 204   fOutput->Add(fHPionMggPt);
 205   fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
 206   fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 207   fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
 208   fOutput->Add(fHPionMggAsym);
 209   fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
 210   fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 211   fHPionMggDgg->SetYTitle("opening angle [grad]");
 212   fOutput->Add(fHPionMggDgg);
 213   const Int_t nbins = 20;
 214   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};
 215   fPtRanges = new TAxis(nbins-1,xbins);
 216   for (Int_t i = 0; i<=nbins; ++i) {
 217     fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
 218     fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
 219     if (i==0)
 220       fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
 221     else if (i==nbins)
 222       fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
 223     else
 224       fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
 225     fOutput->Add(fHPionInvMasses[i]);
 226   }
 227
 228   PostData(1, fOutput);
 229 }
 230
 231 //________________________________________________________________________
 232 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
 233 {
 234   // Called for each event.
 235
 236   if (!InputEvent())
 237     return;
 238
 239   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 240   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 241   if (fEsdEv) {
 242     am->LoadBranch("AliESDRun.");
 243     am->LoadBranch("AliESDHeader.");
 244   } else {
 245     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 246     am->LoadBranch("header");
 247   }
 248
 249   if (!fGeom) { // set misalignment matrices (stored in first event)
 250     fGeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEARV1","EMCAL");
 251     if (fEsdEv) {
 252       for (Int_t i=0; i<4; ++i)
 253         fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
 254     } else {
 255       for (Int_t i=0; i<4; ++i)
 256         fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
 257     }
 258     fGeom->GetEMCGeometry();
 259   }
 260
 261   Int_t cut = 1;
 262   fHCuts->Fill(cut++);
 263
 264   const AliCentrality *centP = InputEvent()->GetCentrality();
 265   Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
 266   fHCent->Fill(cent);
 267   if (cent<fCentFrom||cent>fCentTo)
 268     return;
 269
 270   fHCuts->Fill(cut++);
 271
 272   if (fUseQualFlag) {
 273     if (centP->GetQuality()>0)
 274       return;
 275   }
 276
 277   fHCentQual->Fill(cent);
 278   fHCuts->Fill(cut++);
 279
 280   // count number of accepted events
 281   ++fNEvs;
 282
 283   if (fEsdEv) {
 284     am->LoadBranch("PrimaryVertex.");
 285     am->LoadBranch("SPDVertex.");
 286     am->LoadBranch("TPCVertex.");
 287   } else {
 288     fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 289     am->LoadBranch("vertices");
 290     if (!fAodEv) return;
 291   }
 292
 293   const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
 294   if (!vertex)
 295     return;
 296
 297   fHVertexZ->Fill(vertex->GetZ());
 298
 299   if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
 300     return;
 301
 302   fHCuts->Fill(cut++);
 303   fHVertexZ2->Fill(vertex->GetZ());
 304
 305   fRecPoints   = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
 306   fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
 307   fEsdCells    = 0; // will be set if ESD input used
 308   fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
 309                     //             or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
 310   fAodCells    = 0; // will be set if AOD input used
 311
 312   // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
 313   Bool_t clusattached = 0;
 314   Bool_t recalibrated = 0;
 315   if (1 && !fClusName.IsNull()) {
 316     AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
 317     TObjArray *ts = am->GetTasks();
 318     cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
 319     if (cltask && cltask->GetClusters()) {
 320       fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
 321       clusattached = cltask->GetAttachClusters();
 322       if (cltask->GetCalibData()!=0)
 323         recalibrated = kTRUE;
 324     }
 325   }
 326   if (1 && AODEvent() && !fClusName.IsNull()) {
 327     TList *l = AODEvent()->GetList();
 328     TClonesArray *clus = 0;
 329     if (l) {
 330       clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
 331       fAodClusters = clus;
 332     }
 333   }
 334
 335   if (fEsdEv) { // ESD input mode
 336     if (1 && (!fRecPoints||clusattached)) {
 337       if (!clusattached)
 338         am->LoadBranch("CaloClusters");
 339       TList *l = fEsdEv->GetList();
 340       TClonesArray *clus = 0;
 341       if (l) {
 342         clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 343         fEsdClusters = clus;
 344       }
 345     }
 346     if (1) {
 347       if (!recalibrated)
 348         am->LoadBranch("EMCALCells.");
 349       fEsdCells = fEsdEv->GetEMCALCells();
 350     }
 351   } else if (fAodEv) { // AOD input mode
 352     if (1 && (!fAodClusters || clusattached)) {
 353       if (!clusattached)
 354         am->LoadBranch("caloClusters");
 355       TList *l = fAodEv->GetList();
 356       TClonesArray *clus = 0;
 357       if (l) {
 358         clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
 359         fAodClusters = clus;
 360       }
 361     }
 362     if (1) {
 363       if (!recalibrated)
 364         am->LoadBranch("emcalCells");
 365       fAodCells = fAodEv->GetEMCALCells();
 366     }
 367   } else {
 368     AliFatal("Impossible to not have either pointer to ESD or AOD event");
 369   }
 370
 371   if (1) {
 372     AliDebug(2,Form("fRecPoints   set: %p", fRecPoints));
 373     AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
 374     AliDebug(2,Form("fEsdCells    set: %p", fEsdCells));
 375     AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
 376     AliDebug(2,Form("fAodCells    set: %p", fAodCells));
 377   }
 378
 379   if (fDoAfterburner)
 380     ClusterAfterburner();
 381
 382   FillCellHists();
 383   FillClusHists();
 384   FillPionHists();
 385
 386   PostData(1, fOutput);
 387 }
 388
 389 //________________________________________________________________________
 390 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
 391 {
 392   // Terminate called at the end of analysis.
 393
 394   if (fNtuple) {
 395     TFile *f = OpenFile(1);
 396     if (f)
 397       fNtuple->Write();
 398   }
 399
 400   AliInfo(Form("\n%s: Accepted %lld events", GetName(), fNEvs));
 401 }
 402
 403 //________________________________________________________________________
 404 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
 405 {
 406   // Fill histograms related to cell properties.
 407
 408   AliVCaloCells *cells = fEsdCells;
 409   if (!cells)
 410     cells = fAodCells;
 411
 412   if (cells) {
 413     Int_t cellModCount[4] = {0,0,0,0};
 414     Double_t cellMaxE = 0;
 415     Double_t cellMeanE = 0;
 416     Int_t ncells = cells->GetNumberOfCells();
 417     if (ncells==0)
 418       return;
 419
 420     for (Int_t i = 0; i<ncells; ++i ) {
 421       Int_t absID    = TMath::Abs(cells->GetCellNumber(i));
 422       Double_t cellE = cells->GetAmplitude(i);
 423       fHCellE->Fill(cellE);
 424       if (cellE>cellMaxE)
 425         cellMaxE = cellE;
 426       cellMeanE += cellE;
 427
 428       Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
 429       Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
 430       if (!ret) {
 431         AliError(Form("Could not get cell index for %d", absID));
 432         continue;
 433       }
 434       ++cellModCount[iSM];
 435       Int_t iPhi=-1, iEta=-1;
 436       fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
 437       fHColuRow[iSM]->Fill(iEta,iPhi,1);
 438       fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
 439     }
 440     fHCellH->Fill(cellMaxE);
 441     cellMeanE /= ncells;
 442     fHCellM->Fill(cellMeanE);
 443     fHCellM2->Fill(cellMeanE*ncells/24/48/4); //hard-coded but there is a way to figure out from geometry
 444     for (Int_t i=0; i<4; ++i)
 445       fHCellMult[i]->Fill(cellModCount[i]);
 446   }
 447 }
 448
 449 //________________________________________________________________________
 450 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
 451 {
 452   // Fill histograms related to cluster properties.
 453   Double_t clusterEcc = 0;
 454
 455   Double_t vertex[3] = {0,0,0};
 456   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 457
 458   TObjArray *clusters = fEsdClusters;
 459   if (!clusters)
 460     clusters = fAodClusters;
 461
 462   if (clusters) {
 463     TLorentzVector clusterVec;
 464     Int_t nclus = clusters->GetEntries();
 465     for(Int_t i = 0; i<nclus; ++i) {
 466       AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 467       if (!clus)
 468         continue;
 469       if (!clus->IsEMCAL())
 470         continue;
 471       clus->GetMomentum(clusterVec,vertex);
 472       Double_t maxAxis = 0, minAxis = 0;
 473       GetSigma(clus,maxAxis,minAxis);
 474       if (maxAxis > 0)
 475         clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
 476       clus->SetChi2(clusterEcc); // store ecc in chi2
 477       fHClustEccentricity->Fill(clusterEcc);
 478       fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
 479       fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
 480       fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
 481       fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
 482       fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
 483       if (fNtuple) {
 484         Float_t vals[17];
 485         vals[0]  = InputEvent()->GetRunNumber();
 486         vals[1]  = (((UInt_t)InputEvent()->GetOrbitNumber()  << 12) | (UInt_t)InputEvent()->GetBunchCrossNumber());
 487         if (vals[1]<1)
 488           vals[1] = fNEvs;
 489         vals[2]  = InputEvent()->GetCentrality()->GetCentralityPercentileUnchecked(fCentVar);
 490         vals[3]  = clusterVec.Pt();
 491         vals[4]  = clusterVec.E();
 492         vals[5]  = GetMaxCellEnergy(clus);
 493         vals[6]  = clus->GetNCells();
 494         vals[7]  = clus->GetDistanceToBadChannel();
 495         vals[8]  = clus->GetDispersion();
 496         vals[9]  = clus->GetM20();
 497         vals[10] = clus->GetM02();
 498         vals[11] = clus->Chi2();
 499         vals[12] = clus->GetEmcCpvDistance();
 500         vals[13] = clusterEcc;
 501         vals[14] = maxAxis;
 502         vals[15] = clusterVec.Eta();
 503         vals[16] = clusterVec.Phi();
 504         fNtuple->Fill(vals);
 505       }
 506     }
 507   }
 508 }
 509
 510 //________________________________________________________________________
 511 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
 512 {
 513   // Fill histograms related to pions.
 514
 515   Double_t vertex[3] = {0,0,0};
 516   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 517
 518   TObjArray *clusters = fEsdClusters;
 519   if (!clusters)
 520     clusters = fAodClusters;
 521
 522   if (clusters) {
 523     TLorentzVector clusterVec1;
 524     TLorentzVector clusterVec2;
 525     TLorentzVector pionVec;
 526
 527     Int_t nclus = clusters->GetEntries();
 528     for (Int_t i = 0; i<nclus; ++i) {
 529       AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
 530       if (!clus1)
 531         continue;
 532       if (!clus1->IsEMCAL())
 533         continue;
 534       if (clus1->E()<fMinE)
 535         continue;
 536       if (clus1->GetNCells()<fNminCells)
 537         continue;
 538       if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
 539         continue;
 540       if (clus1->Chi2()<fMinEcc) // eccentricity cut
 541         continue;
 542       clus1->GetMomentum(clusterVec1,vertex);
 543       for (Int_t j = i+1; j<nclus; ++j) {
 544         AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
 545         if (!clus2)
 546           continue;
 547         if (!clus2->IsEMCAL())
 548           continue;
 549         if (clus2->E()<fMinE)
 550           continue;
 551         if (clus2->GetNCells()<fNminCells)
 552           continue;
 553         if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
 554           continue;
 555         if (clus2->Chi2()<fMinEcc) // eccentricity cut
 556           continue;
 557         clus2->GetMomentum(clusterVec2,vertex);
 558         pionVec = clusterVec1 + clusterVec2;
 559         Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
 560         Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
 561         if (pionZgg < fAsymMax) {
 562           fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
 563           fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
 564           fHPionMggAsym->Fill(pionVec.M(),pionZgg);
 565           fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
 566           Int_t bin = fPtRanges->FindBin(pionVec.Pt());
 567           fHPionInvMasses[bin]->Fill(pionVec.M());
 568         }
 569       }
 570     }
 571   }
 572 }
 573
 574 //________________________________________________________________________
 575 void  AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
 576 {
 577   //
 578
 579   AliVCaloCells *cells = fEsdCells;
 580   if (!cells)
 581     cells = fAodCells;
 582
 583   if (!cells)
 584     return;
 585
 586   Int_t ncells = cells->GetNumberOfCells();
 587   if (ncells<=0)
 588     return;
 589
 590   Double_t cellMeanE = 0, cellSigE = 0;
 591   for (Int_t i = 0; i<ncells; ++i ) {
 592     Double_t cellE = cells->GetAmplitude(i);
 593     cellMeanE += cellE;
 594     cellSigE += cellE*cellE;
 595   }
 596   cellMeanE /= ncells;
 597   cellSigE /= ncells;
 598   cellSigE -= cellMeanE*cellMeanE;
 599   if (cellSigE<0)
 600     cellSigE = 0;
 601   cellSigE = TMath::Sqrt(cellSigE / ncells);
 602
 603   Double_t subE = cellMeanE - 7*cellSigE;
 604   if (subE<0)
 605     return;
 606
 607   for (Int_t i = 0; i<ncells; ++i) {
 608     Short_t id=-1;
 609     Double_t amp=0,time=0;
 610     if (!cells->GetCell(i, id, amp, time))
 611       continue;
 612     amp -= cellMeanE;
 613     if (amp<0.001)
 614       amp = 0;
 615     cells->SetCell(i, id, amp, time);
 616   }
 617
 618   TObjArray *clusters = fEsdClusters;
 619   if (!clusters)
 620     clusters = fAodClusters;
 621
 622   if (!clusters)
 623     return;
 624
 625   Int_t nclus = clusters->GetEntries();
 626   for (Int_t i = 0; i<nclus; ++i) {
 627     AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
 628     Int_t nc = clus->GetNCells();
 629     Double_t clusE = 0;
 630     UShort_t ids[100] = {0};
 631     Double_t fra[100] = {0};
 632     for (Int_t j = 0; j<nc; ++j) {
 633       Short_t id = TMath::Abs(clus->GetCellAbsId(j));
 634       Double_t cen = cells->GetCellAmplitude(id);
 635       clusE += cen;
 636       if (cen>0) {
 637         ids[nc] = id;
 638         ++nc;
 639       }
 640     }
 641     if (clusE<=0) {
 642       clusters->RemoveAt(i);
 643       continue;
 644     }
 645
 646     for (Int_t j = 0; j<nc; ++j) {
 647       Short_t id = ids[j];
 648       Double_t cen = cells->GetCellAmplitude(id);
 649       fra[j] = cen/clusE;
 650     }
 651     clus->SetE(clusE);
 652     AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
 653     if (aodclus) {
 654       aodclus->Clear("");
 655       aodclus->SetNCells(nc);
 656       aodclus->SetCellsAmplitudeFraction(fra);
 657       aodclus->SetCellsAbsId(ids);
 658     }
 659   }
 660   clusters->Compress();
 661 }
 662
 663 //________________________________________________________________________
 664 Double_t  AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster)
 665 {
 666   // Get maximum energy of attached cell.
 667
 668   Double_t maxe = 0;
 669   Int_t ncells = cluster->GetNCells();
 670   if (fEsdCells) {
 671     for (Int_t i=0; i<ncells; i++) {
 672       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 673       if (e>maxe)
 674         maxe = e;
 675     }
 676   } else {
 677     for (Int_t i=0; i<ncells; i++) {
 678       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 679       if (e>maxe)
 680         maxe = e;
 681     }
 682   }
 683   return maxe;
 684 }
 685
 686 //________________________________________________________________________
 687 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *cluster, Double_t& sigmaMax, Double_t &sigmaMin)
 688 {
 689   // Calculate the (E) weighted variance along the longer (eigen) axis.
 690
 691   sigmaMax = 0;               // cluster variance along its longer axis
 692   sigmaMin = 0;               // cluster variance along its shorter axis
 693   Double_t Ec = cluster->E(); // cluster energy
 694   Double_t Xc = 0 ;           // cluster first moment along X
 695   Double_t Yc = 0 ;           // cluster first moment along Y
 696   Double_t Sxx = 0 ;          // cluster second central moment along X (variance_X^2)
 697   Double_t Sxy = 0 ;          // cluster second central moment along Y (variance_Y^2)
 698   Double_t Syy = 0 ;          // cluster covariance^2
 699
 700   AliVCaloCells *cells = fEsdCells;
 701   if (!cells)
 702     cells = fAodCells;
 703
 704   if (!cells)
 705     return;
 706
 707   Int_t ncells = cluster->GetNCells();
 708   if (ncells==1)
 709     return;
 710
 711   TVector3 pos;
 712   for(Int_t j=0; j<ncells; ++j) {
 713     Int_t id = TMath::Abs(cluster->GetCellAbsId(j));
 714     fGeom->GetGlobal(id,pos);
 715     Double_t cellen = cells->GetCellAmplitude(id);
 716     Xc  += cellen*pos.X();
 717     Yc  += cellen*pos.Y();
 718     Sxx += cellen*pos.X()*pos.X();
 719     Syy += cellen*pos.Y()*pos.Y();
 720     Sxy += cellen*pos.X()*pos.Y();
 721   }
 722   Xc  /= Ec;
 723   Yc  /= Ec;
 724   Sxx /= Ec;
 725   Syy /= Ec;
 726   Sxy /= Ec;
 727   Sxx -= Xc*Xc;
 728   Syy -= Yc*Yc;
 729   Sxy -= Xc*Yc;
 730   sigmaMax = (Sxx + Syy + TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 731   sigmaMax = TMath::Sqrt(sigmaMax);
 732   sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy))/2.0;
 733   sigmaMin = TMath::Sqrt(sigmaMin);
 734 }
 735