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