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