PWGHF/hfe/AliAnalysisTaskHFECal.cxx

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