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