PWGHF/hfe/AliAnalysisTaskHFECal.cxx

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