PWGHF/hfe/AliAnalysisTaskHFECal.cxx

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