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