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