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