PWGHF/hfe/AliAnalysisTaskHFECal.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // Class for heavy-flavour electron  with EMCal triggered events
  17 // Author: Shingo Sakai
  18
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH2F.h"
  23 #include "TF1.h"
  24 #include "TMath.h"
  25 #include "TCanvas.h"
  26 #include "THnSparse.h"
  27 #include "TLorentzVector.h"
  28 #include "TString.h"
  29 #include "TFile.h"
  30 #include "TGraphErrors.h"
  31
  32 #include "TDatabasePDG.h"
  33
  34 #include "AliAnalysisTask.h"
  35 #include "AliAnalysisManager.h"
  36
  37 #include "AliESDEvent.h"
  38 #include "AliESDHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODHandler.h"
  41
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCParticle.h"
  45
  46 #include "AliAnalysisTaskHFECal.h"
  47 #include "TGeoGlobalMagField.h"
  48 #include "AliLog.h"
  49 #include "AliAnalysisTaskSE.h"
  50 #include "TRefArray.h"
  51 #include "TVector.h"
  52 #include "AliESDInputHandler.h"
  53 #include "AliESDpid.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliPhysicsSelection.h"
  56 #include "AliESDCaloCluster.h"
  57 #include "AliAODCaloCluster.h"
  58 #include "AliEMCALRecoUtils.h"
  59 #include "AliEMCALGeometry.h"
  60 #include "AliGeomManager.h"
  61 #include "stdio.h"
  62 #include "TGeoManager.h"
  63 #include "iostream"
  64 #include "fstream"
  65
  66 #include "AliEMCALTrack.h"
  67 #include "AliMagF.h"
  68
  69 #include "AliKFParticle.h"
  70 #include "AliKFVertex.h"
  71
  72 #include "AliPID.h"
  73 #include "AliPIDResponse.h"
  74 #include "AliHFEcontainer.h"
  75 #include "AliHFEcuts.h"
  76 #include "AliHFEpid.h"
  77 #include "AliHFEpidBase.h"
  78 #include "AliHFEpidQAmanager.h"
  79 #include "AliHFEtools.h"
  80 #include "AliCFContainer.h"
  81 #include "AliCFManager.h"
  82
  83 #include "AliStack.h"
  84
  85 #include "AliCentrality.h"
  86
  87 using namespace std;
  88
  89 ClassImp(AliAnalysisTaskHFECal)
  90 //________________________________________________________________________
  91 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
  92   : AliAnalysisTaskSE(name)
  93   ,fESD(0)
  94   ,fMC(0)
  95   ,stack(0)
  96   ,fGeom(0)
  97   ,fOutputList(0)
  98   ,fqahist(0)
  99   ,fTrackCuts(0)
 100   ,fCuts(0)
 101   ,fIdentifiedAsOutInz(kFALSE)
 102   ,fPassTheEventCut(kFALSE)
 103   ,fRejectKinkMother(kFALSE)
 104   ,fmcData(kFALSE)
 105   ,fVz(0.0)
 106   ,fCFM(0)
 107   ,fPID(0)
 108   ,fPIDqa(0)
 109   ,fOpeningAngleCut(0.1)
 110   ,fMimpTassCut(0.5)
 111   ,fMimNsigassCut(-3)
 112   ,fInvmassCut(0)        // no mass
 113   ,fSetMassConstraint(kTRUE)
 114   ,fSetMassWidthCut(kFALSE)
 115   ,fSetMassNonlinear(kFALSE)
 116   ,fSetKFpart(kTRUE)
 117   ,fNoEvents(0)
 118   ,fEMCAccE(0)
 119   ,hEMCAccE(0)
 120   ,fTrkpt(0)
 121   ,fTrkEovPBef(0)
 122   ,fTrkEovPAft(0)
 123   ,fdEdxBef(0)
 124   ,fdEdxAft(0)
 125   ,fIncpT(0)
 126   ,fIncpTM20(0)
 127   ,fInvmassLS(0)
 128   ,fInvmassULS(0)
 129   ,fInvmassLSmc(0)
 130   ,fInvmassULSmc(0)
 131   ,fInvmassLSreco(0)
 132   ,fInvmassULSreco(0)
 133   ,fInvmassLSmc0(0)
 134   ,fInvmassLSmc1(0)
 135   ,fInvmassLSmc2(0)
 136   ,fInvmassLSmc3(0)
 137   ,fInvmassULSmc0(0)
 138   ,fInvmassULSmc1(0)
 139   ,fInvmassULSmc2(0)
 140   ,fInvmassULSmc3(0)
 141   ,fOpeningAngleLS(0)
 142   ,fOpeningAngleULS(0)
 143   ,fPhotoElecPt(0)
 144   ,fPhoElecPt(0)
 145   ,fPhoElecPtM20(0)
 146   ,fSameElecPt(0)
 147   ,fSameElecPtM20(0)
 148   ,fSemiElecPtM20(0)
 149   ,fTrackPtBefTrkCuts(0)
 150   ,fTrackPtAftTrkCuts(0)
 151   ,fTPCnsigma(0)
 152   ,fCent(0)
 153   ,fEleInfo(0)
 154   ,fElenSigma(0)
 155   /*
 156   ,fClsEBftTrigCut(0)
 157   ,fClsEAftTrigCut(0)
 158   ,fClsEAftTrigCut1(0)
 159   ,fClsEAftTrigCut2(0)
 160   ,fClsEAftTrigCut3(0)
 161   ,fClsEAftTrigCut4(0)
 162   ,fClsETime(0)
 163   ,fClsETime1(0)
 164   ,fTrigTimes(0)
 165   ,fCellCheck(0)
 166   */
 167   ,fInputHFEMC(0)
 168   ,fInputAlle(0)
 169   ,fIncpTMChfe(0)
 170   ,fIncpTMChfeAll(0)
 171   ,fIncpTMCM20hfe(0)
 172   ,fIncpTMCM20hfeAll(0)
 173   ,fIncpTMCM20hfeCheck(0)
 174   ,fInputHFEMC_weight(0)
 175   ,fIncpTMCM20hfeCheck_weight(0)
 176   ,fIncpTMCpho(0)
 177   ,fIncpTMCM20pho(0)
 178   ,fPhoElecPtMC(0)
 179   ,fPhoElecPtMCM20(0)
 180   ,fSameElecPtMC(0)
 181   ,fSameElecPtMCM20(0)
 182   ,fIncpTMCM20pho_pi0e(0)
 183   ,fPhoElecPtMCM20_pi0e(0)
 184   ,fSameElecPtMCM20_pi0e(0)
 185   ,fIncpTMCM20pho_eta(0)
 186   ,fPhoElecPtMCM20_eta(0)
 187   ,fSameElecPtMCM20_eta(0)
 188   ,fIncpTMCpho_pi0e_TPC(0)
 189   ,fPhoElecPtMC_pi0e_TPC(0)
 190   ,fSameElecPtMC_pi0e_TPC(0)
 191   ,fIncpTMCpho_eta_TPC(0)
 192   ,fPhoElecPtMC_eta_TPC(0)
 193   ,fSameElecPtMC_eta_TPC(0)
 194   ,CheckNclust(0)
 195   ,CheckNits(0)
 196   ,CheckDCA(0)
 197   ,Hpi0pTcheck(0)
 198   ,HETApTcheck(0)
 199   ,HphopTcheck(0)
 200   ,HDpTcheck(0)
 201   ,HBpTcheck(0)
 202   ,fpTCheck(0)
 203   ,fMomDtoE(0)
 204   ,fLabelCheck(0)
 205   ,fgeoFake(0)
 206   ,fFakeTrk0(0)
 207   ,fFakeTrk1(0)
 208   ,ftimingEle(0)
 209   ,fIncMaxE(0)
 210   ,fIncReco(0)
 211   ,fPhoReco(0)
 212   ,fSamReco(0)
 213   ,fIncRecoMaxE(0)
 214   ,fPhoRecoMaxE(0)
 215   ,fSamRecoMaxE(0)
 216   ,fPhoVertexReco_HFE(0)
 217   ,fPhoVertexReco_EMCal(0)
 218   ,fPhoVertexReco_Invmass(0)
 219   ,fPhoVertexReco_step0(0)
 220   ,fPhoVertexReco_step1(0)
 221   ,fMatchV0_0(0)
 222   ,fMatchV0_1(0)
 223   ,fMatchMC_0(0)
 224   ,fMatchMC_1(0)
 225   //,fnSigEtaCorr(NULL)
 226 {
 227   //Named constructor
 228
 229   fPID = new AliHFEpid("hfePid");
 230   fTrackCuts = new AliESDtrackCuts();
 231
 232   for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
 233
 234   // Define input and output slots here
 235   // Input slot #0 works with a TChain
 236   DefineInput(0, TChain::Class());
 237   // Output slot #0 id reserved by the base class for AOD
 238   // Output slot #1 writes into a TH1 container
 239   // DefineOutput(1, TH1I::Class());
 240   DefineOutput(1, TList::Class());
 241   //  DefineOutput(3, TTree::Class());
 242 }
 243
 244 //________________________________________________________________________
 245 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
 246   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
 247   ,fESD(0)
 248   ,fMC(0)
 249   ,stack(0)
 250   ,fGeom(0)
 251   ,fOutputList(0)
 252   ,fqahist(0)
 253   ,fTrackCuts(0)
 254   ,fCuts(0)
 255   ,fIdentifiedAsOutInz(kFALSE)
 256   ,fPassTheEventCut(kFALSE)
 257   ,fRejectKinkMother(kFALSE)
 258   ,fmcData(kFALSE)
 259   ,fVz(0.0)
 260   ,fCFM(0)
 261   ,fPID(0)
 262   ,fPIDqa(0)
 263   ,fOpeningAngleCut(0.1)
 264   ,fMimpTassCut(0.5)
 265   ,fMimNsigassCut(-3)
 266   ,fInvmassCut(0)        // no mass
 267   ,fSetMassConstraint(kTRUE)
 268   ,fSetMassWidthCut(kFALSE)
 269   ,fSetMassNonlinear(kFALSE)
 270   ,fSetKFpart(kTRUE)
 271   ,fNoEvents(0)
 272   ,fEMCAccE(0)
 273   ,hEMCAccE(0)
 274   ,fTrkpt(0)
 275   ,fTrkEovPBef(0)
 276   ,fTrkEovPAft(0)
 277   ,fdEdxBef(0)
 278   ,fdEdxAft(0)
 279   ,fIncpT(0)
 280   ,fIncpTM20(0)
 281   ,fInvmassLS(0)
 282   ,fInvmassULS(0)
 283   ,fInvmassLSmc(0)
 284   ,fInvmassULSmc(0)
 285   ,fInvmassLSreco(0)
 286   ,fInvmassULSreco(0)
 287   ,fInvmassLSmc0(0)
 288   ,fInvmassLSmc1(0)
 289   ,fInvmassLSmc2(0)
 290   ,fInvmassLSmc3(0)
 291   ,fInvmassULSmc0(0)
 292   ,fInvmassULSmc1(0)
 293   ,fInvmassULSmc2(0)
 294   ,fInvmassULSmc3(0)
 295   ,fOpeningAngleLS(0)
 296   ,fOpeningAngleULS(0)
 297   ,fPhotoElecPt(0)
 298   ,fPhoElecPt(0)
 299   ,fPhoElecPtM20(0)
 300   ,fSameElecPt(0)
 301   ,fSameElecPtM20(0)
 302   ,fSemiElecPtM20(0)
 303   ,fTrackPtBefTrkCuts(0)
 304   ,fTrackPtAftTrkCuts(0)
 305   ,fTPCnsigma(0)
 306   ,fCent(0)
 307   ,fEleInfo(0)
 308   ,fElenSigma(0)
 309   /*
 310   ,fClsEBftTrigCut(0)
 311   ,fClsEAftTrigCut(0)
 312   ,fClsEAftTrigCut1(0)
 313   ,fClsEAftTrigCut2(0)
 314   ,fClsEAftTrigCut3(0)
 315   ,fClsEAftTrigCut4(0)
 316   ,fClsETime(0)
 317   ,fClsETime1(0)
 318   ,fTrigTimes(0)
 319   ,fCellCheck(0)
 320   */
 321   ,fInputHFEMC(0)
 322   ,fInputAlle(0)
 323   ,fIncpTMChfe(0)
 324   ,fIncpTMChfeAll(0)
 325   ,fIncpTMCM20hfe(0)
 326   ,fIncpTMCM20hfeAll(0)
 327   ,fIncpTMCM20hfeCheck(0)
 328   ,fInputHFEMC_weight(0)
 329   ,fIncpTMCM20hfeCheck_weight(0)
 330   ,fIncpTMCpho(0)
 331   ,fIncpTMCM20pho(0)
 332   ,fPhoElecPtMC(0)
 333   ,fPhoElecPtMCM20(0)
 334   ,fSameElecPtMC(0)
 335   ,fSameElecPtMCM20(0)
 336   ,fIncpTMCM20pho_pi0e(0)
 337   ,fPhoElecPtMCM20_pi0e(0)
 338   ,fSameElecPtMCM20_pi0e(0)
 339   ,fIncpTMCM20pho_eta(0)
 340   ,fPhoElecPtMCM20_eta(0)
 341   ,fSameElecPtMCM20_eta(0)
 342   ,fIncpTMCpho_pi0e_TPC(0)
 343   ,fPhoElecPtMC_pi0e_TPC(0)
 344   ,fSameElecPtMC_pi0e_TPC(0)
 345   ,fIncpTMCpho_eta_TPC(0)
 346   ,fPhoElecPtMC_eta_TPC(0)
 347   ,fSameElecPtMC_eta_TPC(0)
 348   ,CheckNclust(0)
 349   ,CheckNits(0)
 350   ,CheckDCA(0)
 351   ,Hpi0pTcheck(0)
 352   ,HETApTcheck(0)
 353   ,HphopTcheck(0)
 354   ,HDpTcheck(0)
 355   ,HBpTcheck(0)
 356   ,fpTCheck(0)
 357   ,fMomDtoE(0)
 358   ,fLabelCheck(0)
 359   ,fgeoFake(0)
 360   ,fFakeTrk0(0)
 361   ,fFakeTrk1(0)
 362   ,ftimingEle(0)
 363   ,fIncMaxE(0)
 364   ,fIncReco(0)
 365   ,fPhoReco(0)
 366   ,fSamReco(0)
 367   ,fIncRecoMaxE(0)
 368   ,fPhoRecoMaxE(0)
 369   ,fSamRecoMaxE(0)
 370   ,fPhoVertexReco_HFE(0)
 371   ,fPhoVertexReco_EMCal(0)
 372   ,fPhoVertexReco_Invmass(0)
 373   ,fPhoVertexReco_step0(0)
 374   ,fPhoVertexReco_step1(0)
 375   ,fMatchV0_0(0)
 376   ,fMatchV0_1(0)
 377   ,fMatchMC_0(0)
 378   ,fMatchMC_1(0)
 379   //,fnSigEtaCorr(NULL)
 380 {
 381         //Default constructor
 382         fPID = new AliHFEpid("hfePid");
 383
 384         fTrackCuts = new AliESDtrackCuts();
 385
 386        for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
 387
 388         // Constructor
 389         // Define input and output slots here
 390         // Input slot #0 works with a TChain
 391         DefineInput(0, TChain::Class());
 392         // Output slot #0 id reserved by the base class for AOD
 393         // Output slot #1 writes into a TH1 container
 394         // DefineOutput(1, TH1I::Class());
 395         DefineOutput(1, TList::Class());
 396         //DefineOutput(3, TTree::Class());
 397 }
 398 //_________________________________________
 399
 400 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
 401 {
 402   //Destructor
 403
 404   delete fOutputList;
 405   delete fGeom;
 406   delete fPID;
 407   delete fCuts;
 408   delete fCFM;
 409   delete fPIDqa;
 410   delete fTrackCuts;
 411 }
 412 //_________________________________________
 413
 414 void AliAnalysisTaskHFECal::UserExec(Option_t*)
 415 {
 416   //Main loop
 417   //Called for each event
 418
 419   // create pointer to event
 420   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 421   if (!fESD) {
 422     printf("ERROR: fESD not available\n");
 423     return;
 424   }
 425
 426   if(!fCuts){
 427     AliError("HFE cuts not available");
 428     return;
 429   }
 430
 431   if(!fPID->IsInitialized()){
 432     // Initialize PID with the given run number
 433     AliWarning("PID not initialised, get from Run no");
 434     fPID->InitializePID(fESD->GetRunNumber());
 435   }
 436
 437   if(fmcData)fMC = MCEvent();
 438   //AliStack* stack = NULL;
 439   if(fmcData && fMC)stack = fMC->Stack();
 440
 441   Float_t cent = -1.;
 442   AliCentrality *centrality = fESD->GetCentrality();
 443   cent = centrality->GetCentralityPercentile("V0M");
 444
 445   //---- fill MC track info
 446   if(fmcData && fMC)
 447     {
 448     Int_t nParticles = stack->GetNtrack();
 449     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 450       TParticle* particle = stack->Particle(iParticle);
 451       int fPDG = particle->GetPdgCode();
 452       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
 453       double pTMC = particle->Pt();
 454       double proR = particle->R();
 455       double etaMC = particle->Eta();
 456       if(fabs(etaMC)>0.6)continue;
 457       Bool_t mcInDtoE= kFALSE;
 458       Bool_t mcInBtoE= kFALSE;
 459
 460       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 461       //if(!MChijing)printf("not MC hijing");
 462       int iHijing = 1;
 463       if(!MChijing)iHijing = 0;
 464       double mcphoinfo[5];
 465       mcphoinfo[0] = cent;
 466       mcphoinfo[1] = pTMC;
 467       mcphoinfo[2] = iHijing;
 468       if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
 469       if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
 470       if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
 471       if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
 472
 473       if(particle->GetFirstMother()>-1 && fPDG==22)
 474         {
 475          int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 476          if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
 477         }
 478
 479       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
 480         {
 481             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 482             double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();
 483             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
 484             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
 485             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
 486                {
 487                 fInputHFEMC->Fill(cent,pTMC);
 488                 double mcinfo[5];
 489                 mcinfo[0] = cent;
 490                 mcinfo[1] = pTMC;
 491                 mcinfo[2] = 0.0;
 492                 mcinfo[3] = iHijing;
 493                 mcinfo[4] = pTMCparent;
 494                 fInputHFEMC_weight->Fill(mcinfo);
 495                }
 496          }
 497
 498
 499          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
 500
 501       }
 502     }
 503
 504   fNoEvents->Fill(0);
 505
 506   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 507   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 508
 509   Double_t pVtxZ = -999;
 510   pVtxZ = pVtx->GetZ();
 511
 512   if(TMath::Abs(pVtxZ)>10) return;
 513   fNoEvents->Fill(1);
 514
 515   if(fNOtrks<1) return;
 516   fNoEvents->Fill(2);
 517
 518   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 519   if(!pidResponse){
 520     AliDebug(1, "Using default PID Response");
 521     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 522   }
 523
 524   fPID->SetPIDResponse(pidResponse);
 525
 526   fCFM->SetRecEventInfo(fESD);
 527
 528   fCent->Fill(cent);
 529
 530   //if(cent>90.) return;
 531
 532  // Calorimeter info.
 533
 534    FindTriggerClusters();
 535
 536   // make EMCAL array
 537   double maxE = 0.0;
 538   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
 539      {
 540       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
 541       if(clust && clust->IsEMCAL())
 542         {
 543          double clustE = clust->E();
 544          float  emcx[3]; // cluster pos
 545          clust->GetPosition(emcx);
 546          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 547          double emcphi = clustpos.Phi();
 548          double emceta = clustpos.Eta();
 549          double calInfo[5];
 550          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
 551          //fEMCAccE->Fill(calInfo);
 552          hEMCAccE->Fill(cent,clustE);
 553          if(clustE>maxE)maxE = clustE;
 554         }
 555    }
 556
 557   // Track loop
 558   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 559     AliESDtrack* track = fESD->GetTrack(iTracks);
 560     if (!track) {
 561       printf("ERROR: Could not receive track %d\n", iTracks);
 562       continue;
 563     }
 564
 565     //--- Get MC informtion
 566
 567     int parentlabel = 99999;
 568     int parentPID = 99999;
 569     int grand_parentlabel = 99999;
 570     int grand_parentPID = 99999;
 571     Bool_t mcPho = kFALSE;
 572     Bool_t mcDtoE= kFALSE;
 573     Bool_t mcBtoE= kFALSE;
 574     Bool_t mcOrgPi0 = kFALSE;
 575     Bool_t mcOrgEta = kFALSE;
 576     double mcele = -1.;
 577     double mcpT = 0.0;
 578     double mcMompT = 0.0;
 579     //double mcGrandMompT = 0.0;
 580     double mcWeight = -10.0;
 581     double conv_proR = -1.0;
 582
 583     int iHijing = 1;
 584     int mcLabel = -1;
 585
 586     if(fmcData && fMC && stack)
 587       {
 588        Int_t label = TMath::Abs(track->GetLabel());
 589        mcLabel = track->GetLabel();
 590
 591        if(mcLabel>-1)
 592        {
 593
 594                Bool_t MChijing = fMC->IsFromBGEvent(label);
 595                if(!MChijing)iHijing = 0;
 596
 597                TParticle* particle = stack->Particle(label);
 598                int mcpid = particle->GetPdgCode();
 599                mcpT = particle->Pt();
 600                conv_proR = particle->R();
 601                //printf("MCpid = %d",mcpid);
 602                if(particle->GetFirstMother()>-1)
 603                {
 604                        //int parentlabel = particle->GetFirstMother();
 605                        //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 606                        //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
 607                        FindMother(particle, parentlabel, parentPID);
 608                        mcMompT = stack->Particle(parentlabel)->Pt();
 609                        if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
 610                        if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
 611                        if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
 612
 613                        // make D->e pT correlation
 614                        if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
 615
 616                        //cout << "check PID = " << parentPID << endl;
 617                        //cout << "check pho = " << mcPho << endl;
 618                        //cout << "check D or B = " << mcDtoE << endl;
 619                        // pi->e (Dalitz)
 620                        if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
 621                        {
 622                                //cout << "find pi0->e " <<  endl;
 623                                mcOrgPi0 = kTRUE;
 624                                mcWeight = GetMCweight(mcMompT);
 625                        }
 626                        // eta->e (Dalitz)
 627                        if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
 628                        {
 629                                //cout << "find Eta->e " <<  endl;
 630                                mcOrgEta = kTRUE;
 631                                mcWeight = GetMCweightEta(mcMompT);
 632                        }
 633
 634                        // access grand parent
 635                        TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
 636                        //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
 637                        if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
 638                        {
 639                                //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
 640                                //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
 641                                FindMother(particle_parent, grand_parentlabel, grand_parentPID);
 642                                double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
 643                                if(grand_parentPID==111 && mcGrandpT>0.0)
 644                                {
 645                                        // check eta->pi0 decay !
 646                                        int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
 647                                        TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
 648                                        FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
 649                                        if(grand2_parentPID==221)
 650                                        {
 651                                                //cout << "find Eta->e " <<  endl;
 652                                                double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
 653                                                mcOrgEta = kTRUE;
 654                                                mcWeight = GetMCweight(mcGrandpT2);
 655                                                mcMompT = mcGrandpT2;
 656                                        }
 657                                        else
 658                                        {
 659                                                //cout << "find pi0->e " <<  endl;
 660                                                mcOrgPi0 = kTRUE;
 661                                                mcWeight = GetMCweight(mcGrandpT);
 662                                                mcMompT = mcGrandpT;
 663                                        }
 664                                }
 665
 666                                if(grand_parentPID==221 && mcGrandpT>0.0)
 667                                {
 668                                        //cout << "find Eta->e " <<  endl;
 669                                        mcOrgEta = kTRUE;
 670                                        mcOrgPi0 = kFALSE;
 671                                        mcWeight = GetMCweightEta(mcGrandpT);
 672                                        mcMompT = mcGrandpT;
 673                                }
 674                        }
 675                }
 676
 677                //cout << "===================="<<endl;
 678                //cout << "mcDtoE : " << mcDtoE << endl;
 679                //cout << "mcBtoE : " << mcBtoE << endl;
 680                //cout << "mcPho : " << mcPho << endl;
 681
 682                if(fabs(mcpid)==11)mcele= 0.;
 683                //cout << "check e: " << mcele << endl;
 684                if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
 685                //cout << "check D->e: " << mcele << endl;
 686                if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
 687                //cout << "check B->e: " << mcele << endl;
 688                if(fabs(mcpid)==11 && mcPho)mcele= 3.;
 689                //cout << "check Pho->e: " << mcele << endl;
 690
 691                //cout << "check PID " << endl;
 692                if(fabs(mcpid)!=11)
 693                  {
 694                   //cout << "!= 11" << endl;
 695                   //cout << mcpid << endl;
 696                  }
 697                if(mcele==-1)
 698                  {
 699                   //cout << "mcele==-1" << endl;
 700                   //cout << mcele << endl;
 701                   //cout << mcpid << endl;
 702                  }
 703
 704        } // end of mcLabel>-1
 705
 706       } // <------ end of MC info.
 707
 708     //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl;
 709     //printf("weight = %f\n",mcWeight);
 710
 711     if(TMath::Abs(track->Eta())>0.6) continue;
 712     //if(TMath::Abs(track->Pt()<2.5)) continue;
 713     if(TMath::Abs(track->Pt()<0.1)) continue;
 714
 715     int nITS = track->GetNcls(0);
 716
 717     fTrackPtBefTrkCuts->Fill(track->Pt());
 718
 719     /*
 720     UChar_t itsPixel = track->GetITSClusterMap();
 721     cout << "nITS = " << nITS << endl;
 722     if(itsPixel & BIT(0))cout << "1st layer hit" << endl;
 723     if(itsPixel & BIT(1))cout << "2nd layer hit" << endl;
 724     if(itsPixel & BIT(2))cout << "3rd layer hit" << endl;
 725     if(itsPixel & BIT(3))cout << "4th layer hit" << endl;
 726     if(itsPixel & BIT(4))cout << "5th layer hit" << endl;
 727     */
 728
 729     // RecKine: ITSTPC cuts
 730     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 731
 732     //RecKink
 733     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 734       if(track->GetKinkIndex(0) != 0) continue;
 735     }
 736
 737     // RecPrim
 738     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 739
 740     // HFEcuts: ITS layers cuts
 741     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 742     if(mcPho && iHijing==0)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex
 743     if(mcPho && iHijing==1)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex
 744
 745     // HFE cuts: TPC PID cleanup
 746     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 747
 748     int nTPCcl = track->GetTPCNcls();
 749     //int nTPCclF = track->GetTPCNclsF(); // warnings
 750     //int nITS = track->GetNcls(0);
 751
 752     if(mcPho && iHijing==1)fPhoVertexReco_HFE->Fill(track->Pt(),conv_proR); // check MC vertex
 753
 754     fTrackPtAftTrkCuts->Fill(track->Pt());
 755
 756     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 757     pt = track->Pt();
 758     if(pt<2.5)continue;
 759     if(pt<0.1)continue;
 760
 761     //Int_t charge = track->Charge();
 762     fTrkpt->Fill(pt);
 763     mom = track->P();
 764     phi = track->Phi();
 765     eta = track->Eta();
 766     float dca_xy;
 767     float dca_z;
 768     track->GetImpactParameters(dca_xy,dca_z);
 769
 770     dEdx = track->GetTPCsignal();
 771     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 772
 773     //cout << "nSigma correctoon-----" << endl;
 774     //cout << "org = " << fTPCnSigma << endl;
 775     /*
 776     if(!fmcData) // nsigma eta correction
 777        {
 778         double nSigexpCorr = NsigmaCorrection(eta,cent);
 779         fTPCnSigma -= nSigexpCorr;
 780        }
 781     */
 782     //cout << "correction = " << fTPCnSigma << endl;
 783
 784     //--- track cluster match
 785
 786         double ncells = -1.0;
 787         double m20 = -1.0;
 788         double m02 = -1.0;
 789         double disp = -1.0;
 790         double rmatch = -1.0;
 791         double nmatch = -1.0;
 792         //double oppstatus = 0.0;
 793         double emctof = 0.0;
 794         Bool_t MaxEmatch = kFALSE;
 795
 796     Int_t clsId = track->GetEMCALcluster();
 797     if (clsId>0){
 798       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 799       if(clust && clust->IsEMCAL()){
 800
 801                 double clustE = clust->E();
 802                 if(clustE==maxE)MaxEmatch = kTRUE;
 803                 eop = clustE/fabs(mom);
 804
 805                 //double clustT = clust->GetTOF();
 806                 ncells = clust->GetNCells();
 807                 m02 = clust->GetM02();
 808                 m20 = clust->GetM20();
 809                 disp = clust->GetDispersion();
 810                 double delphi = clust->GetTrackDx();
 811                 double deleta = clust->GetTrackDz();
 812                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 813                 nmatch = clust->GetNTracksMatched();
 814                 emctof = clust->GetTOF();
 815                 //cout << "emctof = " << emctof << endl;
 816                 cout << "eop org = "<< eop << endl;
 817                 double eoporg =  eop;
 818                 if(fmcData)
 819                   {
 820                    double mceopcorr = MCEopMeanCorrection(pt,cent);
 821                    eop += mceopcorr;
 822                   }
 823                 cout << "eop corr = " << eop << endl;
 824
 825                   double valdedx[16];
 826                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 827                   //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
 828                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy,  valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT;
 829                   valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl;
 830                   valdedx[15] = mcele;
 831                   fEleInfo->Fill(valdedx);
 832
 833       }
 834     }
 835
 836     //Get Cal info PID response
 837     /*
 838     double eop2;
 839     double ss[4];
 840     Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
 841     if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
 842       {
 843        double valEop[3];
 844        valEop[0] = cent;
 845        valEop[1] = pt;
 846        valEop[2] = nSigmaEop;
 847        fElenSigma->Fill(valEop);
 848       }
 849     */
 850    // --- tarck cut & e ID
 851
 852     if(nITS<2.5)continue;
 853     if(nTPCcl<100)continue;
 854
 855     CheckNclust->Fill(nTPCcl);
 856     CheckNits->Fill(nITS);
 857     CheckDCA->Fill(dca_xy,dca_z);
 858     // check production vertex of photons
 859
 860
 861     fdEdxBef->Fill(mom,fTPCnSigma);
 862     fTPCnsigma->Fill(mom,fTPCnSigma);
 863     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 864
 865     //--- track accepted by HFE
 866
 867     Int_t pidpassed = 1;
 868     AliHFEpidObject hfetrack;
 869     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 870     hfetrack.SetRecTrack(track);
 871     double binct = 10.5;
 872     if((0.0< cent) && (cent<5.0)) binct = 0.5;
 873     if((5.0< cent) && (cent<10.0)) binct = 1.5;
 874     if((10.0< cent) && (cent<20.0)) binct = 2.5;
 875     if((20.0< cent) && (cent<30.0)) binct = 3.5;
 876     if((30.0< cent) && (cent<40.0)) binct = 4.5;
 877     if((40.0< cent) && (cent<50.0)) binct = 5.5;
 878     if((50.0< cent) && (cent<60.0)) binct = 6.5;
 879     if((60.0< cent) && (cent<70.0)) binct = 7.5;
 880     if((70.0< cent) && (cent<80.0)) binct = 8.5;
 881     if((80.0< cent) && (cent<90.0)) binct = 9.5;
 882     if((90.0< cent) && (cent<100.0)) binct = 10.5;
 883
 884     hfetrack.SetCentrality((int)binct); //added
 885     hfetrack.SetPbPb();
 886     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 887
 888     if(pidpassed==0) continue; // nSigma rejection
 889
 890     //--- photonic ID
 891
 892     Bool_t fFlagPhotonicElec = kFALSE;
 893     Bool_t fFlagConvinatElec = kFALSE;
 894
 895     if(fSetKFpart)
 896       {
 897        SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
 898       }
 899     else
 900      {
 901       SelectPhotonicElectron2(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
 902       }
 903
 904     //--- check reco eff. with TPC
 905     double phoval[5];
 906     phoval[0] = cent;
 907     phoval[1] = pt;
 908     phoval[2] = fTPCnSigma;
 909     phoval[3] = iHijing;
 910     phoval[4] = mcMompT;
 911
 912     //if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && mcele>-1 && mcPho && mcOrgPi0)
 913     if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta))
 914       {
 915         if(iHijing==1)mcWeight = 1.0;
 916         if(mcOrgPi0)
 917           {
 918            fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
 919            if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 920            if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 921           }
 922         if(mcOrgEta)
 923           {
 924            fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight);
 925            if(fFlagPhotonicElec) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 926            if(fFlagConvinatElec) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 927           }
 928       }
 929
 930     //--- matchin check
 931     double emcphimim = 1.396;
 932     double emcphimax = 3.14;
 933     if(phi>emcphimim && phi<emcphimax)
 934       {
 935        if(fFlagPhotonicElec)fMatchV0_0->Fill(pt);  // data
 936        if(mcele>2.1)fMatchMC_0->Fill(pt,mcWeight); // MC
 937
 938        if(eop>=0.0) // have a match
 939          {
 940           if(fFlagPhotonicElec)fMatchV0_1->Fill(pt);
 941           if(mcele>2.1)fMatchMC_1->Fill(pt,mcWeight);
 942          }
 943       }
 944
 945     //+++++++  E/p cut ++++++++++++++++
 946
 947     if(eop<0.9 || eop>1.3)continue;
 948
 949     fTrkEovPAft->Fill(pt,eop);
 950     fdEdxAft->Fill(mom,fTPCnSigma);
 951
 952     // Fill real data
 953     fIncpT->Fill(cent,pt);
 954     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 955     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 956
 957     if(m20>0.0 && m20<0.3)
 958       {
 959        fIncpTM20->Fill(cent,pt);
 960        ftimingEle->Fill(pt,emctof);
 961        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
 962        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
 963        if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt);
 964      }
 965
 966
 967     // MC
 968     // check label for electron candidiates
 969
 970     int idlabel = 1;
 971     if(mcLabel==0)idlabel = 0;
 972     fLabelCheck->Fill(pt,idlabel);
 973     if(mcLabel==0)fgeoFake->Fill(phi,eta);
 974
 975     if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
 976        {
 977         fFakeTrk0->Fill(cent,pt);
 978        }
 979
 980     if(mcele>-1) // select MC electrons
 981       {
 982
 983           fIncpTMChfeAll->Fill(cent,pt);
 984           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
 985           if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
 986
 987        if(mcBtoE || mcDtoE) // select B->e & D->e
 988          {
 989           fIncpTMChfe->Fill(cent,pt);
 990           if(m20>0.0 && m20<0.3)
 991             {
 992                 //cout << "MC label = " << mcLabel << endl;
 993                 fIncpTMCM20hfe->Fill(cent,pt);
 994                 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
 995                 fIncpTMCM20hfeCheck_weight->Fill(phoval);
 996             }
 997          }
 998
 999        if(mcPho) // select photonic electrons
1000         {
1001
1002          fIncpTMCpho->Fill(phoval);
1003          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
1004          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
1005
1006          if(m20>0.0 && m20<0.3)
1007            {
1008             fIncpTMCM20pho->Fill(phoval);
1009             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
1010             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
1011             // pi0->g->e
1012             if(mcWeight>-1)
1013               {
1014                if(iHijing==1)mcWeight = 1.0;
1015                if(mcOrgPi0)
1016                  {
1017                   fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
1018                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1019                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1020
1021                   // check production vertex
1022                   fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR);
1023                   if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR);
1024
1025                   }
1026                // --- eta
1027                if(mcOrgEta)
1028                  {
1029                   fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
1030                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
1031                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
1032                  }
1033                 // check production vertex
1034               }
1035            }
1036         }
1037       }
1038  }
1039  PostData(1, fOutputList);
1040 }
1041 //_________________________________________
1042 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
1043 {
1044   //--- Check MC
1045
1046   //Bool_t mcData = kFALSE;
1047   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
1048     {
1049      fmcData = kTRUE;
1050      printf("+++++ MC Data available");
1051     }
1052   if(fmcData)
1053     {
1054      printf("++++++++= MC analysis \n");
1055     }
1056   else
1057    {
1058      printf("++++++++ real data analysis \n");
1059    }
1060
1061    printf("+++++++ QA hist %d",fqahist);
1062
1063   //---- Geometry
1064   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1065
1066   //--------Initialize PID
1067   //fPID->SetHasMCData(kFALSE);
1068   fPID->SetHasMCData(fmcData);
1069   if(!fPID->GetNumberOfPIDdetectors())
1070     {
1071       fPID->AddDetector("TPC", 0);
1072       //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task
1073     }
1074
1075   Double_t params[4];
1076   const char *cutmodel;
1077   cutmodel = "pol0";
1078   params[0] = -1.0; //sigma min
1079   double maxnSig = 3.0;
1080   if(fmcData)
1081     {
1082      params[0] = -5.0; //sigma min
1083      maxnSig = 5.0;
1084     }
1085
1086   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
1087
1088   fPID->SortDetectors();
1089   fPIDqa = new AliHFEpidQAmanager();
1090   fPIDqa->Initialize(fPID);
1091
1092   //------- fcut --------------
1093   fCuts = new AliHFEcuts();
1094   fCuts->CreateStandardCuts();
1095   //fCuts->SetMinNClustersTPC(100);
1096   fCuts->SetMinNClustersTPC(90);
1097   fCuts->SetMinRatioTPCclusters(0.6);
1098   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1099   //fCuts->SetMinNClustersITS(3);
1100   fCuts->SetMinNClustersITS(2);
1101   fCuts->SetProductionVertex(0,50,0,50);
1102   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1103   fCuts->SetCheckITSLayerStatus(kFALSE);
1104   fCuts->SetVertexRange(10.);
1105   fCuts->SetTOFPIDStep(kFALSE);
1106   //fCuts->SetPtRange(2, 50);
1107   fCuts->SetPtRange(0.1, 50);
1108   fCuts->SetMaxImpactParam(3.,3.);
1109
1110   //--------Initialize correction Framework and Cuts
1111   fCFM = new AliCFManager;
1112   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1113   fCFM->SetNStepParticle(kNcutSteps);
1114   for(Int_t istep = 0; istep < kNcutSteps; istep++)
1115     fCFM->SetParticleCutsList(istep, NULL);
1116
1117   if(!fCuts){
1118     AliWarning("Cuts not available. Default cuts will be used");
1119     fCuts = new AliHFEcuts;
1120     fCuts->CreateStandardCuts();
1121   }
1122   fCuts->Initialize(fCFM);
1123
1124   //---------Output Tlist
1125   fOutputList = new TList();
1126   fOutputList->SetOwner();
1127   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1128
1129   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1130   fOutputList->Add(fNoEvents);
1131
1132   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
1133   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
1134   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
1135   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
1136   //if(fqahist==1)fOutputList->Add(fEMCAccE);
1137
1138   hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1139   fOutputList->Add(hEMCAccE);
1140
1141   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1142   fOutputList->Add(fTrkpt);
1143
1144   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1145   fOutputList->Add(fTrackPtBefTrkCuts);
1146
1147   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1148   fOutputList->Add(fTrackPtAftTrkCuts);
1149
1150   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1151   fOutputList->Add(fTPCnsigma);
1152
1153   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1154   fOutputList->Add(fTrkEovPBef);
1155
1156   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1157   fOutputList->Add(fTrkEovPAft);
1158
1159   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
1160   fOutputList->Add(fdEdxBef);
1161
1162   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
1163   fOutputList->Add(fdEdxAft);
1164
1165   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
1166   fOutputList->Add(fIncpT);
1167
1168   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1169   fOutputList->Add(fIncpTM20);
1170
1171   Int_t nBinspho[9] =  { 10,  30,   600, 60,   50,    4,  40,   8,  30};
1172   Double_t minpho[9] = {  0.,  0., -0.1, 40,  0, -0.5,   0,-1.5,   0};
1173   Double_t maxpho[9] = {100., 30.,  0.5, 100,   1,  3.5,   2, 6.5,  30};
1174
1175   fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1176   if(fqahist==1)fOutputList->Add(fInvmassLS);
1177
1178   fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1179   if(fqahist==1)fOutputList->Add(fInvmassULS);
1180
1181   fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1182   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
1183
1184   fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1185   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
1186
1187   fInvmassLSreco = new TH2D("fInvmassLSreco", "Inv mass of LS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1188   fInvmassLSreco->Sumw2();
1189   fOutputList->Add(fInvmassLSreco);
1190
1191   fInvmassULSreco = new TH2D("fInvmassULSreco", "Inv mass of ULS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1192   fInvmassULSreco->Sumw2();
1193   fOutputList->Add(fInvmassULSreco);
1194
1195   fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1196   fInvmassLSmc0->Sumw2();
1197   fOutputList->Add(fInvmassLSmc0);
1198
1199   fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1200   fInvmassLSmc1->Sumw2();
1201   fOutputList->Add(fInvmassLSmc1);
1202
1203   fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1204   fInvmassLSmc2->Sumw2();
1205   fOutputList->Add(fInvmassLSmc2);
1206
1207   fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1208   fInvmassLSmc3->Sumw2();
1209   fOutputList->Add(fInvmassLSmc3);
1210
1211   fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1212   fInvmassULSmc0->Sumw2();
1213   fOutputList->Add(fInvmassULSmc0);
1214
1215   fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1216   fInvmassULSmc1->Sumw2();
1217   fOutputList->Add(fInvmassULSmc1);
1218
1219   fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1220   fInvmassULSmc2->Sumw2();
1221   fOutputList->Add(fInvmassULSmc2);
1222
1223   fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1224   fInvmassULSmc3->Sumw2();
1225   fOutputList->Add(fInvmassULSmc3);
1226
1227   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1228   fOutputList->Add(fOpeningAngleLS);
1229
1230   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1231   fOutputList->Add(fOpeningAngleULS);
1232
1233   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1234   fOutputList->Add(fPhotoElecPt);
1235
1236   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
1237   fOutputList->Add(fPhoElecPt);
1238
1239   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
1240   fOutputList->Add(fPhoElecPtM20);
1241
1242   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
1243   fOutputList->Add(fSameElecPt);
1244
1245   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
1246   fOutputList->Add(fSameElecPtM20);
1247
1248   fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1249   fOutputList->Add(fSemiElecPtM20);
1250
1251   fCent = new TH1F("fCent","Centrality",200,0,100) ;
1252   fOutputList->Add(fCent);
1253
1254   // Make common binning
1255   const Double_t kMinP = 0.;
1256   const Double_t kMaxP = 20.;
1257
1258   //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
1259   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, dca_xy, dca_z, M20, Disp, Centrality, select)
1260   Int_t nBins[16] =  {  100,     7,  60,    20,    90,  100,   25,     60,    60,  100,    40,  10,  250,  100, 100,    8};
1261   Double_t min[16] = {kMinP,  -0.5, 1.0,  -1.0,  -5.0,    0,    0,   -3.0,  -3.0,  0.0,   0.0,   0,    0,    0,  80, -1.5};
1262   Double_t max[16] = {kMaxP,   6.5, 4.0,   1.0,   4.0,  2.0, 0.05,    3.0,   3.0,  1.0,  20.0, 100,  100,  2.0, 180,  6.5};
1263   //fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
1264   fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;DCA_xy;DCA_z;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
1265   //if(fqahist==1)fOutputList->Add(fEleInfo);
1266
1267   // Make common binning
1268   Int_t nBinsEop[3] =  { 10, 50, 100};
1269   Double_t minEop[3] = {  0,  0,  -5};
1270   Double_t maxEop[3] = {100, 50,   5};
1271   fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1272   fOutputList->Add(fElenSigma);
1273
1274
1275   //<---  Trigger info
1276   /*
1277   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1278   fOutputList->Add(fClsEBftTrigCut);
1279
1280   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1281   fOutputList->Add(fClsEAftTrigCut);
1282
1283   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1284   fOutputList->Add(fClsEAftTrigCut1);
1285
1286   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1287   fOutputList->Add(fClsEAftTrigCut2);
1288
1289   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1290   fOutputList->Add(fClsEAftTrigCut3);
1291
1292   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1293   fOutputList->Add(fClsEAftTrigCut4);
1294
1295   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1296   fOutputList->Add(fClsETime);
1297
1298   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1299   fOutputList->Add(fClsETime1);
1300
1301   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1302   fOutputList->Add(fTrigTimes);
1303
1304   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1305   fOutputList->Add(fCellCheck);
1306   */
1307   //<---------- MC
1308
1309   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1310   fOutputList->Add(fInputHFEMC);
1311
1312   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1313   fOutputList->Add(fInputAlle);
1314
1315   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1316   fOutputList->Add(fIncpTMChfe);
1317
1318   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1319   fOutputList->Add(fIncpTMChfeAll);
1320
1321   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1322   fOutputList->Add(fIncpTMCM20hfe);
1323
1324   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1325   fOutputList->Add(fIncpTMCM20hfeAll);
1326
1327   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1328   fOutputList->Add(fIncpTMCM20hfeCheck);
1329
1330   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1331   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};
1332   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};
1333
1334   fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1335   fOutputList->Add(fInputHFEMC_weight);
1336
1337   fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1338   fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1339
1340   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1341   fOutputList->Add(fIncpTMCpho);
1342
1343   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1344   fOutputList->Add(fIncpTMCM20pho);
1345
1346   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1347   fOutputList->Add(fPhoElecPtMC);
1348
1349   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1350   fOutputList->Add(fPhoElecPtMCM20);
1351
1352   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1353   fOutputList->Add(fSameElecPtMC);
1354
1355   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1356   fOutputList->Add(fSameElecPtMCM20);
1357
1358   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1359   fIncpTMCM20pho_pi0e->Sumw2();
1360   fOutputList->Add(fIncpTMCM20pho_pi0e);
1361
1362   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1363   fPhoElecPtMCM20_pi0e->Sumw2();
1364   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1365
1366   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1367   fSameElecPtMCM20_pi0e->Sumw2();
1368   fOutputList->Add(fSameElecPtMCM20_pi0e);
1369  //
1370   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1371   fIncpTMCM20pho_eta->Sumw2();
1372   fOutputList->Add(fIncpTMCM20pho_eta);
1373
1374   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1375   fPhoElecPtMCM20_eta->Sumw2();
1376   fOutputList->Add(fPhoElecPtMCM20_eta);
1377
1378   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1379   fSameElecPtMCM20_eta->Sumw2();
1380   fOutputList->Add(fSameElecPtMCM20_eta);
1381   // ------------
1382   fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1383   fIncpTMCpho_pi0e_TPC->Sumw2();
1384   fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1385
1386   fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1387   fPhoElecPtMC_pi0e_TPC->Sumw2();
1388   fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1389
1390   fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1391   fSameElecPtMC_pi0e_TPC->Sumw2();
1392   fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1393
1394   fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1395   fIncpTMCpho_eta_TPC->Sumw2();
1396   fOutputList->Add(fIncpTMCpho_eta_TPC);
1397
1398   fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1399   fPhoElecPtMC_eta_TPC->Sumw2();
1400   fOutputList->Add(fPhoElecPtMC_eta_TPC);
1401
1402   fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1403   fSameElecPtMC_eta_TPC->Sumw2();
1404   fOutputList->Add(fSameElecPtMC_eta_TPC);
1405
1406
1407   //-------------
1408
1409   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1410   fOutputList->Add(CheckNclust);
1411
1412   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1413   fOutputList->Add(CheckNits);
1414
1415   CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1416   fOutputList->Add(CheckDCA);
1417   /*
1418   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1419   fOutputList->Add(Hpi0pTcheck);
1420
1421   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1422   fOutputList->Add(HETApTcheck);
1423   */
1424
1425   Int_t nBinspho3[3] =  { 200, 100, 3};
1426   Double_t minpho3[3] = {  0.,  0., -0.5};
1427   Double_t maxpho3[3] = {100., 50., 2.5};
1428
1429   Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1430   fOutputList->Add(Hpi0pTcheck);
1431
1432   HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1433   fOutputList->Add(HETApTcheck);
1434   //--
1435   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1436   fOutputList->Add(HphopTcheck);
1437   //
1438   HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1439   fOutputList->Add(HDpTcheck);
1440
1441   HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1442   fOutputList->Add(HBpTcheck);
1443   //
1444
1445   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1446   fOutputList->Add(fpTCheck);
1447
1448   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1449   fOutputList->Add(fMomDtoE);
1450
1451   fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1452   fOutputList->Add(fLabelCheck);
1453
1454   fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1455   fOutputList->Add(fgeoFake);
1456
1457   fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1458   fOutputList->Add(fFakeTrk0);
1459
1460   fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1461   fOutputList->Add(fFakeTrk1);
1462
1463   ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1464   fOutputList->Add(ftimingEle);
1465
1466   // eta correction
1467   // note: parameters 01/31new.TPCnSigmaEtaDep
1468   // 70-90 delta_eta = 0.2
1469
1470   double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55};
1471   double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done)
1472   double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done)
1473   double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done)
1474   double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done)
1475   double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done)
1476   double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done)
1477   double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done)
1478
1479   fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1480   fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1481   fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1482   fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1483   fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1484   fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1485   fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1486
1487   fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1488   fOutputList->Add(fIncMaxE);
1489
1490   fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
1491   fOutputList->Add(fIncReco);
1492
1493   fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1494   fOutputList->Add(fPhoReco);
1495
1496   fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1497   fOutputList->Add(fSamReco);
1498
1499   fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1500   fOutputList->Add(fIncRecoMaxE);
1501
1502   fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
1503   fOutputList->Add(fPhoRecoMaxE);
1504
1505   fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
1506   fOutputList->Add(fSamRecoMaxE);
1507
1508   fPhoVertexReco_HFE = new TH2D("fPhoVertexReco_HFE","photon production Vertex mass selection",40,0,20,250,0,50);
1509   fOutputList->Add(fPhoVertexReco_HFE);
1510
1511   fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
1512   fOutputList->Add(fPhoVertexReco_EMCal);
1513
1514   fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
1515   fOutputList->Add(fPhoVertexReco_Invmass);
1516
1517   fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
1518   fPhoVertexReco_step0->Sumw2();
1519   fOutputList->Add(fPhoVertexReco_step0);
1520
1521   fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
1522   fPhoVertexReco_step1->Sumw2();
1523   fOutputList->Add(fPhoVertexReco_step1);
1524
1525   fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20);
1526   fOutputList->Add(fMatchV0_0);
1527
1528   fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20);
1529   fOutputList->Add(fMatchV0_1);
1530
1531   fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20);
1532   fOutputList->Add(fMatchMC_0);
1533
1534   fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20);
1535   fOutputList->Add(fMatchMC_1);
1536
1537   PostData(1,fOutputList);
1538 }
1539
1540 //________________________________________________________________________
1541 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1542 {
1543   // Info("Terminate");
1544         AliAnalysisTaskSE::Terminate();
1545 }
1546
1547 //________________________________________________________________________
1548 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1549 {
1550   // Check single track cuts for a given cut step
1551   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1552   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1553   return kTRUE;
1554 }
1555 //_________________________________________
1556 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta)
1557 {
1558   //Identify non-heavy flavour electrons using Invariant mass method
1559
1560   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1561   fTrackCuts->SetRequireTPCRefit(kTRUE);
1562   fTrackCuts->SetRequireITSRefit(kTRUE);
1563   fTrackCuts->SetEtaRange(-0.9,0.9);
1564   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1565   fTrackCuts->SetMinNClustersTPC(90);
1566
1567   //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1568   Double_t bfield = fESD->GetMagneticField();
1569   Double_t emass = 0.51*0.001; // (0.51 MeV)
1570
1571   double ptEle = track->Pt();  //add
1572   if(ibgevent==0 && w > 0.0)
1573      {
1574       fpTCheck->Fill(ptEle,w);
1575      }
1576
1577   Bool_t flagPhotonicElec = kFALSE;
1578   Bool_t flagConvinatElec = kFALSE;
1579
1580   int p1 = 0;
1581   if(mce==3)
1582      {
1583        Int_t label = TMath::Abs(track->GetLabel());
1584        TParticle* particle = stack->Particle(label);
1585        p1 = particle->GetFirstMother();
1586      }
1587
1588   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1589   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1590     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1591     if (!trackAsso) {
1592       printf("ERROR: Could not receive track %d\n", jTracks);
1593       continue;
1594     }
1595     if(itrack==jTracks)continue;
1596     int jbgevent = 0;
1597
1598     int p2 = 0;
1599     if(mce==3)
1600     {
1601       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1602       TParticle* particle2 = stack->Particle(label2);
1603       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1604       if(MChijing_ass)jbgevent =1;
1605       if(particle2->GetFirstMother()>-1)
1606          p2 = particle2->GetFirstMother();
1607     }
1608
1609     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1610     Double_t mass=999., width = -999;
1611     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1612
1613     //ptPrim = track->Pt();
1614     ptPrim = ptEle;
1615
1616     dEdxAsso = trackAsso->GetTPCsignal();
1617     ptAsso = trackAsso->Pt();
1618     Int_t chargeAsso = trackAsso->Charge();
1619     Int_t charge = track->Charge();
1620
1621
1622     if(ptAsso <fMimpTassCut) continue;
1623     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1624     //if(dEdxAsso <65 || dEdxAsso>100) continue;
1625     double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1626     //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1627     //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1628     if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1629     //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
1630
1631     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1632     if(charge>0) fPDGe1 = -11;
1633     if(chargeAsso>0) fPDGe2 = -11;
1634
1635     //printf("chargeAsso = %d\n",chargeAsso);
1636     //printf("charge = %d\n",charge);
1637     if(charge == chargeAsso) fFlagLS = kTRUE;
1638     if(charge != chargeAsso) fFlagULS = kTRUE;
1639
1640     //printf("fFlagLS = %d\n",fFlagLS);
1641     //printf("fFlagULS = %d\n",fFlagULS);
1642     //printf("\n");
1643
1644     AliKFParticle::SetField(bfield);
1645     AliKFParticle ge1(*track, fPDGe1);
1646     AliKFParticle ge2(*trackAsso, fPDGe2);
1647
1648     if(fSetMassNonlinear)
1649       {
1650        ge1.SetNonlinearMassConstraint(emass);
1651        ge2.SetNonlinearMassConstraint(emass);
1652       }
1653
1654     AliKFParticle recg(ge1, ge2);
1655
1656     /*
1657     // vertex
1658     AliKFVertex primV(*pVtx);
1659     primV += recg;
1660     primV -= ge1;
1661     primV -= ge2;
1662     recg.SetProductionVertex(primV);
1663     */
1664
1665     // check chi2
1666     if(recg.GetNDF()<1) continue;
1667     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1668     Double_t chi2cut = 3.0;
1669
1670     // mass.
1671     if(fSetMassConstraint)
1672       {
1673        recg.SetMassConstraint(0,0.0001);
1674        chi2cut = 30.0;
1675       }
1676     recg.GetMass(mass,width);
1677
1678     if(fSetMassWidthCut && width>1e10)continue;
1679
1680         // angle
1681     openingAngle = ge1.GetAngle(ge2);
1682     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1683     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1684
1685     double ishower = 0;
1686     if(shower>0.0 && shower<0.3)ishower = 1;
1687
1688     double phoinfo[9];
1689     phoinfo[0] = cent;
1690     phoinfo[1] = ptPrim;
1691     phoinfo[2] = mass;
1692     phoinfo[3] = nSig;
1693     //phoinfo[3] = dEdxAsso;
1694     phoinfo[4] = openingAngle;
1695     phoinfo[5] = ishower;
1696     phoinfo[6] = ep;
1697     phoinfo[7] = mce;
1698     phoinfo[8] = ptAsso;
1699
1700     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1701     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1702     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1703     if(fFlagULS && ibgevent==0 && jbgevent==0)
1704        {
1705          fInvmassULSmc->Fill(phoinfo,w);
1706        }
1707     //printf("fInvmassCut %f\n",fInvmassCut);
1708     //printf("openingAngle %f\n",fOpeningAngleCut);
1709
1710     // angle cut
1711     if(openingAngle > fOpeningAngleCut) continue;
1712     // chi2 cut
1713     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1714     if(chi2recg>chi2cut) continue;
1715
1716     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1717     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1718
1719     // for real data
1720     //printf("mce =%f\n",mce);
1721     if(mce<-0.5) // mce==-1. is real
1722        {
1723          //printf("Real data\n");
1724          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1725                flagPhotonicElec = kTRUE;
1726               }
1727          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1728                flagConvinatElec = kTRUE;
1729               }
1730         }
1731     // for MC data
1732     else
1733        {
1734          //printf("MC data\n");
1735
1736          if(w>0.0)
1737            {
1738            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1739            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1740            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1741            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1742            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1743            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1744            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1745            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1746            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1747           }
1748
1749          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1750          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
1751                flagPhotonicElec = kTRUE;
1752               }
1753          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1754                flagConvinatElec = kTRUE;
1755               }
1756         }
1757
1758   }
1759   fFlagPhotonicElec = flagPhotonicElec;
1760   fFlagConvinatElec = flagConvinatElec;
1761
1762 }
1763
1764 //_________________________________________
1765 void AliAnalysisTaskHFECal::SelectPhotonicElectron2(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta)
1766 {
1767   //Identify non-heavy flavour electrons using Invariant mass method
1768
1769   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1770   fTrackCuts->SetRequireTPCRefit(kTRUE);
1771   fTrackCuts->SetRequireITSRefit(kTRUE);
1772   fTrackCuts->SetEtaRange(-0.9,0.9);
1773   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1774   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1775   fTrackCuts->SetMinNClustersTPC(90);
1776
1777   Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
1778   Double_t bfield = fESD->GetMagneticField();
1779
1780   double ptEle = track->Pt();  //add
1781   if(ibgevent==0 && w > 0.0)
1782      {
1783       fpTCheck->Fill(ptEle,w);
1784      }
1785
1786   Bool_t flagPhotonicElec = kFALSE;
1787   Bool_t flagConvinatElec = kFALSE;
1788
1789   int p1 = 0;
1790   if(mce==3)
1791      {
1792        Int_t label = TMath::Abs(track->GetLabel());
1793        TParticle* particle = stack->Particle(label);
1794        p1 = particle->GetFirstMother();
1795      }
1796
1797   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1798   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1799     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1800     if (!trackAsso) {
1801       printf("ERROR: Could not receive track %d\n", jTracks);
1802       continue;
1803     }
1804     if(itrack==jTracks)continue;
1805     int jbgevent = 0;
1806
1807     int p2 = 0;
1808     if(mce==3)
1809     {
1810       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1811       TParticle* particle2 = stack->Particle(label2);
1812       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1813       if(MChijing_ass)jbgevent =1;
1814       if(particle2->GetFirstMother()>-1)
1815          p2 = particle2->GetFirstMother();
1816     }
1817
1818     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1819     //Double_t mass=999., width = -999;
1820     Double_t mass=999.;
1821     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1822
1823     //ptPrim = track->Pt();
1824     ptPrim = ptEle;
1825
1826     dEdxAsso = trackAsso->GetTPCsignal();
1827     ptAsso = trackAsso->Pt();
1828     Int_t chargeAsso = trackAsso->Charge();
1829     Int_t charge = track->Charge();
1830
1831
1832     if(ptAsso <0.5) continue;
1833     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1834     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1835
1836     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1837     if(charge>0) fPDGe1 = -11;
1838     if(chargeAsso>0) fPDGe2 = -11;
1839
1840     //printf("chargeAsso = %d\n",chargeAsso);
1841     //printf("charge = %d\n",charge);
1842     if(charge == chargeAsso) fFlagLS = kTRUE;
1843     if(charge != chargeAsso) fFlagULS = kTRUE;
1844
1845     // mass cal 0
1846     double xt1 = 0.0, xt2 = 0.0;
1847     Double_t p1at[3] = {0,0,0};
1848     Double_t p2at[3] = {0,0,0};
1849     //double dca12 = trackAsso->GetDCA(track,bfield,xt2,xt1);           //DCA track1-track2
1850     double kHasdcaT1 = track->GetPxPyPzAt(xt1,bfield,p1at);             //Track1
1851     double kHasdcaT2 = trackAsso->GetPxPyPzAt(xt2,bfield,p2at);         //Track2
1852    if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
1853    //cout << "dca = " << dca12 << endl;
1854    // 3D
1855    TLorentzVector electron1;
1856    TLorentzVector electron2;
1857    TLorentzVector mother;
1858
1859    electron1.SetXYZM(p1at[0], p1at[1], p1at[2], eMass);
1860    electron2.SetXYZM(p2at[0], p2at[1], p2at[2], eMass);
1861    openingAngle = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
1862
1863    mother      = electron1 + electron2;
1864    //double invmassAtDCA  = mother.M();
1865
1866    // 2D
1867    TLorentzVector electron1_2D;
1868    TLorentzVector electron2_2D;
1869    TLorentzVector mother_2D;
1870
1871    double pT1at = sqrt(pow(p1at[0],2)+pow(p1at[1],2));
1872    double pT2at = sqrt(pow(p2at[0],2)+pow(p2at[1],2));
1873
1874    electron1_2D.SetXYZM(pT1at, 0.0, p1at[2], eMass);
1875    electron2_2D.SetXYZM(pT2at, 0.0, p2at[2], eMass);
1876
1877    mother_2D      = electron1_2D + electron2_2D;
1878    double invmassAtDCA_2D  = mother_2D.M();
1879    mass = invmassAtDCA_2D;
1880
1881         // angle
1882     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1883     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1884
1885     double ishower = 0;
1886     if(shower>0.0 && shower<0.3)ishower = 1;
1887
1888     double phoinfo[9];
1889     phoinfo[0] = cent;
1890     phoinfo[1] = ptPrim;
1891     phoinfo[2] = mass;
1892     phoinfo[3] = nSig;
1893     //phoinfo[3] = dEdxAsso;
1894     phoinfo[4] = openingAngle;
1895     phoinfo[5] = ishower;
1896     phoinfo[6] = ep;
1897     phoinfo[7] = mce;
1898     phoinfo[8] = ptAsso;
1899
1900     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1901     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1902     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1903     if(fFlagULS && ibgevent==0 && jbgevent==0)
1904        {
1905          fInvmassULSmc->Fill(phoinfo,w);
1906        }
1907     //printf("fInvmassCut %f\n",fInvmassCut);
1908     //printf("openingAngle %f\n",fOpeningAngleCut);
1909
1910     // angle cut
1911     if(openingAngle > fOpeningAngleCut) continue;
1912     // chi2 cut
1913     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1914     //if(chi2recg>chi2cut) continue;
1915
1916     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1917     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1918
1919     // for real data
1920     //printf("mce =%f\n",mce);
1921     if(mce<-0.5) // mce==-1. is real
1922        {
1923          //printf("Real data\n");
1924          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1925                flagPhotonicElec = kTRUE;
1926               }
1927          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1928                flagConvinatElec = kTRUE;
1929               }
1930         }
1931     // for MC data
1932     else
1933        {
1934          //printf("MC data\n");
1935
1936          if(w>0.0)
1937            {
1938            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1939            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1940            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1941            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1942            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1943            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1944            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1945            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1946            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1947           }
1948
1949          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1950          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1951                flagPhotonicElec = kTRUE;
1952               }
1953          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1954                flagConvinatElec = kTRUE;
1955               }
1956         }
1957
1958   }
1959   fFlagPhotonicElec = flagPhotonicElec;
1960   fFlagConvinatElec = flagConvinatElec;
1961
1962 }
1963
1964 //-------------------------------------------
1965
1966 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1967 {
1968  //int label = 99999;
1969  //int pid = 99999;
1970
1971  if(part->GetFirstMother()>-1)
1972    {
1973     label = part->GetFirstMother();
1974     pid = stack->Particle(label)->GetPdgCode();
1975    }
1976    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1977 }
1978
1979 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1980 {
1981         double weight = 1.0;
1982
1983         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1984         {
1985                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1986         }
1987         else
1988         {
1989                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1990         }
1991   return weight;
1992 }
1993
1994 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1995 {
1996   double weight = 1.0;
1997
1998   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1999   return weight;
2000 }
2001
2002
2003 //_________________________________________
2004 void AliAnalysisTaskHFECal::FindTriggerClusters()
2005 {
2006   //cout << "finding trigger patch" << endl;
2007   // constants
2008   const int nModuleCols = 2;
2009   const int nModuleRows = 5;
2010   const int nColsFeeModule = 48;
2011   const int nRowsFeeModule = 24;
2012   const int nColsFaltroModule = 24;
2013   const int nRowsFaltroModule = 12;
2014   //const int faltroWidthMax = 20;
2015
2016   // part 1, trigger extraction -------------------------------------
2017   Int_t globCol, globRow;
2018   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
2019   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
2020
2021   //Int_t trigtimes[faltroWidthMax];
2022   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2023   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2024   //Double_t fTrigCutLow = 6;
2025   //Double_t fTrigCutHigh = 10;
2026   Double_t fTimeCutLow =  469e-09;
2027   Double_t fTimeCutHigh = 715e-09;
2028
2029   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2030
2031   // erase trigger maps
2032   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
2033   {
2034     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
2035     {
2036       ftriggersCut[i][j] = 0;
2037       ftriggers[i][j] = 0;
2038       ftriggersTime[i][j] = 0;
2039     }
2040   }
2041
2042   Int_t iglobCol=0, iglobRow=0;
2043   // go through triggers
2044   if( fCaloTrigger->GetEntries() > 0 )
2045   {
2046     // needs reset
2047     fCaloTrigger->Reset();
2048     while( fCaloTrigger->Next() )
2049     {
2050       fCaloTrigger->GetPosition( globCol, globRow );
2051       fCaloTrigger->GetNL0Times( ntimes );
2052       /*
2053       // no L0s
2054       if( ntimes < 1 )   continue;
2055       // get precise timings
2056       fCaloTrigger->GetL0Times( trigtimes );
2057       trigInCut = 0;
2058       for(Int_t i = 0; i < ntimes; i++ )
2059       {
2060         // save the first trigger time in channel
2061         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
2062           triggersTime[globCol][globRow] = trigtimes[i];
2063         //printf("trigger times: %d\n",trigtimes[i]);
2064         // check if in cut
2065         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
2066           trigInCut = 1;
2067
2068         fTrigTimes->Fill(trigtimes[i]);
2069       }
2070       */
2071
2072       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
2073       Int_t bit = 0;
2074       fCaloTrigger->GetTriggerBits(bit);
2075       //cout << "bit = " << bit << endl;
2076
2077       Int_t ts = 0;
2078       fCaloTrigger->GetL1TimeSum(ts);
2079       //cout << "ts = " << ts << endl;
2080       if (ts > 0)ftriggers[globCol][globRow] = 1;
2081       // number of triggered channels in event
2082       nTrigChannel++;
2083       // ... inside cut
2084       if(ts>0 && (bit >> 6 & 0x1))
2085       {
2086         iglobCol = globCol;
2087         iglobRow = globRow;
2088         nTrigChannelCut++;
2089         //cout << "ts cut = " << ts << endl;
2090         //cout << "globCol = " << globCol << endl;
2091         //cout << "globRow = " << globRow << endl;
2092         ftriggersCut[globCol][globRow] = 1;
2093       }
2094
2095     } // calo trigger entries
2096   } // has calo trigger entries
2097
2098   // part 2 go through the clusters here -----------------------------------
2099   //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
2100   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
2101   Short_t cellAddr, nSACell;
2102   Int_t mclabel;
2103   //Int_t nSACell, iSACell, mclabel;
2104   Int_t iSACell;
2105   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
2106   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2107
2108   TRefArray *fCaloClusters = new TRefArray();
2109   fESD->GetEMCALClusters( fCaloClusters );
2110   nCluster = fCaloClusters->GetEntries();
2111
2112
2113   // save all cells times if there are clusters
2114   if( nCluster > 0 ){
2115     // erase time map
2116     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2117       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2118         cellTime[i][j] = 0.;
2119         cellEnergy[i][j] = 0.;
2120       }
2121     }
2122
2123     // get cells
2124     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2125     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2126     nSACell = fCaloCells->GetNumberOfCells();
2127     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2128       // get the cell info *fCal
2129       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2130
2131       // get cell position
2132       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2133       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2134
2135       // convert co global phi eta
2136       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2137       geta = ieta + nColsFeeModule*(nSupMod%2);
2138
2139       // save cell time and energy
2140       cellTime[geta][gphi] = cellTimeT;
2141       cellEnergy[geta][gphi] = cellAmp;
2142
2143     }
2144   }
2145
2146   Int_t nClusterTrig, nClusterTrigCut;
2147   UShort_t *cellAddrs;
2148   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2149   Float_t clsPos[3] = {0.,0.,0.};
2150
2151   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2152   {
2153     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2154     if(!cluster || !cluster->IsEMCAL()) continue;
2155
2156     // get cluster cells
2157     nCell = cluster->GetNCells();
2158
2159     // get cluster energy
2160     clsE = cluster->E();
2161
2162     // get cluster position
2163     cluster->GetPosition(clsPos);
2164     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2165     clsEta = clsPosVec.Eta();
2166     clsPhi = clsPosVec.Phi();
2167
2168     // get the cell addresses
2169     cellAddrs = cluster->GetCellsAbsId();
2170
2171     // check if the cluster contains cell, that was marked as triggered
2172     nClusterTrig = 0;
2173     nClusterTrigCut = 0;
2174
2175     // loop the cells to check, if cluser in acceptance
2176     // any cluster with a cell outside acceptance is not considered
2177     for( iCell = 0; iCell < nCell; iCell++ )
2178     {
2179      // check hot cell
2180      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
2181
2182       // get cell position
2183       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2184       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2185
2186       // convert co global phi eta
2187       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2188       geta = ieta + nColsFeeModule*(nSupMod%2);
2189
2190       if( cellTime[geta][gphi] > 0. ){
2191         clusterTime += cellTime[geta][gphi];
2192         gCell++;
2193       }
2194
2195       // get corresponding FALTRO
2196       fphi = gphi / 2;
2197       feta = geta / 2;
2198
2199         //cout << "fphi = " << fphi << endl;
2200         //cout << "feta = " << feta << endl;
2201
2202       // try to match with a triggered
2203       if( ftriggers[feta][fphi]==1)
2204       {  nClusterTrig++;
2205       }
2206       if( ftriggersCut[feta][fphi]==1)
2207       { nClusterTrigCut++;
2208       }
2209
2210       //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2211
2212     } // cells
2213
2214
2215     if( gCell > 0 )
2216       clusterTime = clusterTime / (Double_t)gCell;
2217     // fix the reconstriction code time 100ns jumps
2218     if( fESD->GetBunchCrossNumber() % 4 < 2 )
2219       clusterTime -= 0.0000001;
2220
2221     //fClsETime->Fill(clsE,clusterTime);
2222     //fClsEBftTrigCut->Fill(clsE);
2223
2224     if(nClusterTrig>0){
2225       //fClsETime1->Fill(clsE,clusterTime);
2226     }
2227
2228     if(nClusterTrig>0){
2229       cluster->SetChi2(1);
2230       //fClsEAftTrigCut1->Fill(clsE);
2231     }
2232
2233     if(nClusterTrigCut>0){
2234       cluster->SetChi2(2);
2235       //fClsEAftTrigCut2->Fill(clsE);
2236     }
2237
2238     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2239     {
2240       cluster->SetChi2(3);
2241       //fClsEAftTrigCut3->Fill(clsE);
2242     }
2243
2244     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2245     {
2246       // cluster->SetChi2(4);
2247       //fClsEAftTrigCut4->Fill(clsE);
2248     }
2249     if(nClusterTrigCut<1)
2250     {
2251       cluster->SetChi2(0);
2252
2253       //fClsEAftTrigCut->Fill(clsE);
2254     }
2255
2256   } // clusters
2257 }
2258
2259 // <-------- only MC correction
2260 double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
2261 {
2262   TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2263   TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2264
2265  double shift = 0.0;
2266
2267  if(central>0 && central<=10)
2268    {
2269     fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2270     fcorr1->SetParameters(9.91e-01,3.466,2.344);
2271    }
2272  else if(central>10 && central<=20)
2273    {
2274     fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2275     fcorr1->SetParameters(0.975,2.276,1.501e-01);
2276    }
2277  else if(central>20 && central<=30)
2278    {
2279     fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2280     fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2281    }
2282  else if(central>30 && central<=40)
2283    {
2284     fcorr0->SetParameters(1.00,1.466,2.305e-1);
2285     fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2286    }
2287  else if(central>40 && central<=50)
2288    {
2289     fcorr0->SetParameters(1.00,1.422,1.518e-01);
2290     fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2291    }
2292
2293  else if(central>50 && central<=70)
2294    {
2295     fcorr0->SetParameters(0.989,2.495,2.167);
2296     fcorr1->SetParameters(0.961,1.734,1.438e-01);
2297    }
2298  else if(central>70 && central<=100)
2299    {
2300     fcorr0->SetParameters(0.981,-3.373,3.93327);
2301     fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
2302    }
2303
2304
2305  shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2306
2307  fcorr0->Delete();
2308  fcorr1->Delete();
2309
2310  return shift;
2311 }
2312
2313 // <-------- only Data correction
2314 double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2315 {
2316  int icent = 0;
2317
2318  if(central>=0 && central<10)
2319    {
2320     icent = 0;
2321    }
2322  else if(central>=10 && central<20)
2323   {
2324    icent = 1;
2325   }
2326  else if(central>=20 && central<30)
2327   {
2328    icent = 2;
2329   }
2330  else if(central>=30 && central<40)
2331   {
2332    icent = 3;
2333   }
2334  else if(central>=40 && central<50)
2335   {
2336    icent = 4;
2337   }
2338  else if(central>=50 && central<70)
2339   {
2340    icent = 5;
2341   }
2342  else
2343   {
2344    icent = 6;
2345   }
2346
2347  double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2348
2349  //cout << "eta correction"<< endl;
2350  //cout << "cent = "<< central<< endl;
2351  //cout << "icent = "<< icent << endl;
2352  //cout << "shift = "<< shift << endl;
2353
2354  return shift;
2355
2356 }
2357
2358