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