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