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 "TMath.h"
  24 #include "TCanvas.h"
  25 #include "THnSparse.h"
  26 #include "TLorentzVector.h"
  27 #include "TString.h"
  28 #include "TFile.h"
  29
  30 #include "AliAnalysisTask.h"
  31 #include "AliAnalysisManager.h"
  32
  33 #include "AliESDEvent.h"
  34 #include "AliESDHandler.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODHandler.h"
  37
  38 #include "AliMCEventHandler.h"
  39 #include "AliMCEvent.h"
  40 #include "AliMCParticle.h"
  41
  42 #include "AliAnalysisTaskHFECal.h"
  43 #include "TGeoGlobalMagField.h"
  44 #include "AliLog.h"
  45 #include "AliAnalysisTaskSE.h"
  46 #include "TRefArray.h"
  47 #include "TVector.h"
  48 #include "AliESDInputHandler.h"
  49 #include "AliESDpid.h"
  50 #include "AliESDtrackCuts.h"
  51 #include "AliPhysicsSelection.h"
  52 #include "AliESDCaloCluster.h"
  53 #include "AliAODCaloCluster.h"
  54 #include "AliEMCALRecoUtils.h"
  55 #include "AliEMCALGeometry.h"
  56 #include "AliGeomManager.h"
  57 #include "stdio.h"
  58 #include "TGeoManager.h"
  59 #include "iostream"
  60 #include "fstream"
  61
  62 #include "AliEMCALTrack.h"
  63 #include "AliMagF.h"
  64
  65 #include "AliKFParticle.h"
  66 #include "AliKFVertex.h"
  67
  68 #include "AliPID.h"
  69 #include "AliPIDResponse.h"
  70 #include "AliHFEcontainer.h"
  71 #include "AliHFEcuts.h"
  72 #include "AliHFEpid.h"
  73 #include "AliHFEpidBase.h"
  74 #include "AliHFEpidQAmanager.h"
  75 #include "AliHFEtools.h"
  76 #include "AliCFContainer.h"
  77 #include "AliCFManager.h"
  78
  79 #include "AliStack.h"
  80
  81 #include "AliCentrality.h"
  82
  83 ClassImp(AliAnalysisTaskHFECal)
  84 //________________________________________________________________________
  85 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
  86   : AliAnalysisTaskSE(name)
  87   ,fESD(0)
  88   ,fMC(0)
  89   ,fGeom(0)
  90   ,fOutputList(0)
  91   ,fqahist(1)
  92   ,fTrackCuts(0)
  93   ,fCuts(0)
  94   ,fIdentifiedAsOutInz(kFALSE)
  95   ,fPassTheEventCut(kFALSE)
  96   ,fRejectKinkMother(kFALSE)
  97   ,fmcData(kFALSE)
  98   ,fVz(0.0)
  99   ,fCFM(0)
 100   ,fPID(0)
 101   ,fPIDqa(0)
 102   ,fOpeningAngleCut(0.1)
 103   ,fInvmassCut(0.01)
 104   ,fNoEvents(0)
 105   ,fEMCAccE(0)
 106   ,fTrkpt(0)
 107   ,fTrkEovPBef(0)
 108   ,fTrkEovPAft(0)
 109   ,fdEdxBef(0)
 110   ,fdEdxAft(0)
 111   ,fIncpT(0)
 112   ,fIncpTM20(0)
 113   ,fInvmassLS(0)
 114   ,fInvmassULS(0)
 115   ,fInvmassLSmc(0)
 116   ,fInvmassULSmc(0)
 117   ,fOpeningAngleLS(0)
 118   ,fOpeningAngleULS(0)
 119   ,fPhotoElecPt(0)
 120   ,fPhoElecPt(0)
 121   ,fPhoElecPtM20(0)
 122   ,fSameElecPt(0)
 123   ,fSameElecPtM20(0)
 124   ,fTrackPtBefTrkCuts(0)
 125   ,fTrackPtAftTrkCuts(0)
 126   ,fTPCnsigma(0)
 127   ,fCent(0)
 128   ,fEleInfo(0)
 129   /*
 130   ,fClsEBftTrigCut(0)
 131   ,fClsEAftTrigCut(0)
 132   ,fClsEAftTrigCut1(0)
 133   ,fClsEAftTrigCut2(0)
 134   ,fClsEAftTrigCut3(0)
 135   ,fClsEAftTrigCut4(0)
 136   ,fClsETime(0)
 137   ,fClsETime1(0)
 138   ,fTrigTimes(0)
 139   ,fCellCheck(0)
 140   */
 141   ,fInputHFEMC(0)
 142   ,fInputAlle(0)
 143   ,fIncpTMChfe(0)
 144   ,fIncpTMChfeAll(0)
 145   ,fIncpTMCM20hfe(0)
 146   ,fIncpTMCM20hfeAll(0)
 147   ,fIncpTMCpho(0)
 148   ,fIncpTMCM20pho(0)
 149   ,fPhoElecPtMC(0)
 150   ,fPhoElecPtMCM20(0)
 151   ,fSameElecPtMC(0)
 152   ,fSameElecPtMCM20(0)
 153   ,fIncpTMCM20pho_pi0e(0)
 154   ,fPhoElecPtMCM20_pi0e(0)
 155   ,fSameElecPtMCM20_pi0e(0)
 156   ,CheckNclust(0)
 157   ,CheckNits(0)
 158   ,Hpi0pTcheck(0)
 159   ,HphopTcheck(0)
 160   ,fMomDtoE(0)
 161 {
 162   //Named constructor
 163
 164   fPID = new AliHFEpid("hfePid");
 165   fTrackCuts = new AliESDtrackCuts();
 166
 167   // Define input and output slots here
 168   // Input slot #0 works with a TChain
 169   DefineInput(0, TChain::Class());
 170   // Output slot #0 id reserved by the base class for AOD
 171   // Output slot #1 writes into a TH1 container
 172   // DefineOutput(1, TH1I::Class());
 173   DefineOutput(1, TList::Class());
 174   //  DefineOutput(3, TTree::Class());
 175 }
 176
 177 //________________________________________________________________________
 178 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
 179   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
 180   ,fESD(0)
 181   ,fMC(0)
 182   ,fGeom(0)
 183   ,fOutputList(0)
 184   ,fqahist(1)
 185   ,fTrackCuts(0)
 186   ,fCuts(0)
 187   ,fIdentifiedAsOutInz(kFALSE)
 188   ,fPassTheEventCut(kFALSE)
 189   ,fRejectKinkMother(kFALSE)
 190   ,fmcData(kFALSE)
 191   ,fVz(0.0)
 192   ,fCFM(0)
 193   ,fPID(0)
 194   ,fPIDqa(0)
 195   ,fOpeningAngleCut(0.1)
 196   ,fInvmassCut(0.01)
 197   ,fNoEvents(0)
 198   ,fEMCAccE(0)
 199   ,fTrkpt(0)
 200   ,fTrkEovPBef(0)
 201   ,fTrkEovPAft(0)
 202   ,fdEdxBef(0)
 203   ,fdEdxAft(0)
 204   ,fIncpT(0)
 205   ,fIncpTM20(0)
 206   ,fInvmassLS(0)
 207   ,fInvmassULS(0)
 208   ,fInvmassLSmc(0)
 209   ,fInvmassULSmc(0)
 210   ,fOpeningAngleLS(0)
 211   ,fOpeningAngleULS(0)
 212   ,fPhotoElecPt(0)
 213   ,fPhoElecPt(0)
 214   ,fPhoElecPtM20(0)
 215   ,fSameElecPt(0)
 216   ,fSameElecPtM20(0)
 217   ,fTrackPtBefTrkCuts(0)
 218   ,fTrackPtAftTrkCuts(0)
 219   ,fTPCnsigma(0)
 220   ,fCent(0)
 221   ,fEleInfo(0)
 222   /*
 223   ,fClsEBftTrigCut(0)
 224   ,fClsEAftTrigCut(0)
 225   ,fClsEAftTrigCut1(0)
 226   ,fClsEAftTrigCut2(0)
 227   ,fClsEAftTrigCut3(0)
 228   ,fClsEAftTrigCut4(0)
 229   ,fClsETime(0)
 230   ,fClsETime1(0)
 231   ,fTrigTimes(0)
 232   ,fCellCheck(0)
 233   */
 234   ,fInputHFEMC(0)
 235   ,fInputAlle(0)
 236   ,fIncpTMChfe(0)
 237   ,fIncpTMChfeAll(0)
 238   ,fIncpTMCM20hfe(0)
 239   ,fIncpTMCM20hfeAll(0)
 240   ,fIncpTMCpho(0)
 241   ,fIncpTMCM20pho(0)
 242   ,fPhoElecPtMC(0)
 243   ,fPhoElecPtMCM20(0)
 244   ,fSameElecPtMC(0)
 245   ,fSameElecPtMCM20(0)
 246   ,fIncpTMCM20pho_pi0e(0)
 247   ,fPhoElecPtMCM20_pi0e(0)
 248   ,fSameElecPtMCM20_pi0e(0)
 249   ,CheckNclust(0)
 250   ,CheckNits(0)
 251   ,Hpi0pTcheck(0)
 252   ,HphopTcheck(0)
 253   ,fMomDtoE(0)
 254 {
 255         //Default constructor
 256         fPID = new AliHFEpid("hfePid");
 257
 258         fTrackCuts = new AliESDtrackCuts();
 259
 260         // Constructor
 261         // Define input and output slots here
 262         // Input slot #0 works with a TChain
 263         DefineInput(0, TChain::Class());
 264         // Output slot #0 id reserved by the base class for AOD
 265         // Output slot #1 writes into a TH1 container
 266         // DefineOutput(1, TH1I::Class());
 267         DefineOutput(1, TList::Class());
 268         //DefineOutput(3, TTree::Class());
 269 }
 270 //_________________________________________
 271
 272 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
 273 {
 274   //Destructor
 275
 276   delete fOutputList;
 277   delete fGeom;
 278   delete fPID;
 279   delete fCuts;
 280   delete fCFM;
 281   delete fPIDqa;
 282   delete fTrackCuts;
 283 }
 284 //_________________________________________
 285
 286 void AliAnalysisTaskHFECal::UserExec(Option_t*)
 287 {
 288   //Main loop
 289   //Called for each event
 290
 291   // create pointer to event
 292   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 293   if (!fESD) {
 294     printf("ERROR: fESD not available\n");
 295     return;
 296   }
 297
 298   if(!fCuts){
 299     AliError("HFE cuts not available");
 300     return;
 301   }
 302
 303   if(!fPID->IsInitialized()){
 304     // Initialize PID with the given run number
 305     AliWarning("PID not initialised, get from Run no");
 306     fPID->InitializePID(fESD->GetRunNumber());
 307   }
 308
 309   if(fmcData)fMC = MCEvent();
 310   AliStack* stack = NULL;
 311   if(fmcData && fMC)stack = fMC->Stack();
 312
 313   Float_t cent = -1.;
 314   AliCentrality *centrality = fESD->GetCentrality();
 315   cent = centrality->GetCentralityPercentile("V0M");
 316
 317   //---- fill MC track info
 318   if(fmcData && fMC)
 319     {
 320     Int_t nParticles = stack->GetNtrack();
 321     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 322       TParticle* particle = stack->Particle(iParticle);
 323       int fPDG = particle->GetPdgCode();
 324       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
 325       double pTMC = particle->Pt();
 326       double proR = particle->R();
 327       double etaMC = particle->Eta();
 328       if(fabs(etaMC)>0.6)continue;
 329       Bool_t mcInDtoE= kFALSE;
 330       Bool_t mcInBtoE= kFALSE;
 331
 332       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 333       //if(!MChijing)printf("not MC hijing");
 334       int iHijing = 1;
 335       if(!MChijing)iHijing = 0;
 336       if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
 337
 338       if(particle->GetFirstMother()>-1 && fPDG==22)
 339         {
 340          int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 341          if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
 342         }
 343
 344       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
 345         {
 346             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 347             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
 348             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
 349             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
 350          }
 351
 352          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
 353
 354       }
 355     }
 356
 357   fNoEvents->Fill(0);
 358
 359   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 360   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 361
 362   Double_t pVtxZ = -999;
 363   pVtxZ = pVtx->GetZ();
 364
 365   if(TMath::Abs(pVtxZ)>10) return;
 366   fNoEvents->Fill(1);
 367
 368   if(fNOtrks<2) return;
 369   fNoEvents->Fill(2);
 370
 371   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 372   if(!pidResponse){
 373     AliDebug(1, "Using default PID Response");
 374     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 375   }
 376
 377   fPID->SetPIDResponse(pidResponse);
 378
 379   fCFM->SetRecEventInfo(fESD);
 380
 381   //Float_t cent = -1.;
 382   //AliCentrality *centrality = fESD->GetCentrality();
 383   //cent = centrality->GetCentralityPercentile("V0M");
 384   fCent->Fill(cent);
 385
 386   //if(cent>90.) return;
 387
 388  // Calorimeter info.
 389
 390    FindTriggerClusters();
 391
 392   // make EMCAL array
 393   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
 394      {
 395       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
 396       if(clust->IsEMCAL())
 397         {
 398          double clustE = clust->E();
 399          float  emcx[3]; // cluster pos
 400          clust->GetPosition(emcx);
 401          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 402          double emcphi = clustpos.Phi();
 403          double emceta = clustpos.Eta();
 404          double calInfo[5];
 405          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
 406          //if(clustE>1.5)fEMCAccE->Fill(calInfo);
 407          //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
 408         }
 409    }
 410
 411   // Track loop
 412   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 413     AliESDtrack* track = fESD->GetTrack(iTracks);
 414     if (!track) {
 415       printf("ERROR: Could not receive track %d\n", iTracks);
 416       continue;
 417     }
 418
 419     Bool_t mcPho = kFALSE;
 420     Bool_t mcDtoE= kFALSE;
 421     Bool_t mcBtoE= kFALSE;
 422     double mcele = -1.;
 423     double mcpT = 0.0;
 424     double mcMompT = 0.0;
 425     //double mcGrandMompT = 0.0;
 426     double mcWeight = -10.0;
 427
 428     int iHijing = 1;
 429
 430     if(fmcData && fMC && stack)
 431       {
 432        Int_t label = TMath::Abs(track->GetLabel());
 433
 434        Bool_t MChijing = fMC->IsFromBGEvent(label);
 435        if(!MChijing)iHijing = 0;
 436
 437        TParticle* particle = stack->Particle(label);
 438        int mcpid = particle->GetPdgCode();
 439        mcpT = particle->Pt();
 440        //printf("MCpid = %d",mcpid);
 441        if(particle->GetFirstMother()>-1)
 442          {
 443           int parentlabel = particle->GetFirstMother();
 444           int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 445           mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
 446           if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
 447           if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
 448           if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
 449
 450           // make D->e pT correlation
 451           if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
 452
 453           // pi->e (Dalitz)
 454           if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
 455               {
 456                /*
 457                  if(mcMompT>0.0 && mcMompT<5.0)
 458                    {
 459                     mcWeight = 0.323*mcMompT/(TMath::Exp(-1.6+0.767*mcMompT+0.0285*mcMompT*mcMompT));
 460                    }
 461                  else
 462                    {
 463                     mcWeight = 115.0/(0.718*mcMompT*TMath::Power(mcMompT,3.65));
 464                    }
 465                 */
 466                mcWeight = GetMCweight(mcMompT);
 467               }
 468
 469           // access grand parent pi0->g->e
 470           TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
 471           if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
 472             {
 473              int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
 474              if(grand_parentPID==111)
 475                 {
 476                  //mcGrandMompT = stack->Particle(particle_parent->GetFirstMother())->Pt();
 477                  double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
 478                  //if(mcGrandMompT>0.0 && mcGrandMompT<5.0)
 479                  if(pTtmp>0.0)
 480                    {
 481                     /*
 482                     if(pTtmp<5.0)
 483                     {
 484                      mcWeight = 0.323*pTtmp/(TMath::Exp(-1.6+0.767*pTtmp+0.0285*pTtmp*pTtmp));
 485                     }
 486                    else
 487                    {
 488                     mcWeight = 115.0/(0.718*pTtmp*TMath::Power(pTtmp,3.65));
 489                    }
 490                    */
 491                    mcWeight = GetMCweight(pTtmp);
 492                   }
 493                 }
 494             }
 495          }
 496        if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
 497        if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
 498        if(fabs(mcpid)==11 && mcPho)mcele= 3.;
 499       }
 500
 501     //printf("weight = %f\n",mcWeight);
 502
 503     if(TMath::Abs(track->Eta())>0.6) continue;
 504     if(TMath::Abs(track->Pt()<2.0)) continue;
 505
 506     fTrackPtBefTrkCuts->Fill(track->Pt());
 507     // RecKine: ITSTPC cuts
 508     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 509
 510     //RecKink
 511     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 512       if(track->GetKinkIndex(0) != 0) continue;
 513     }
 514
 515     // RecPrim
 516     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 517
 518     // HFEcuts: ITS layers cuts
 519     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 520
 521     // HFE cuts: TPC PID cleanup
 522     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 523
 524     int nTPCcl = track->GetTPCNcls();
 525     int nTPCclF = track->GetTPCNclsF();
 526     int nITS = track->GetNcls(0);
 527
 528     fTrackPtAftTrkCuts->Fill(track->Pt());
 529
 530     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 531     pt = track->Pt();
 532     if(pt<2.0)continue;
 533
 534     // Track extrapolation
 535
 536     Int_t charge = track->Charge();
 537     fTrkpt->Fill(pt);
 538     mom = track->P();
 539     phi = track->Phi();
 540     eta = track->Eta();
 541     dEdx = track->GetTPCsignal();
 542     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 543
 544         double ncells = -1.0;
 545         double m20 = -1.0;
 546         double m02 = -1.0;
 547         double disp = -1.0;
 548         double rmatch = -1.0;
 549         double nmatch = -1.0;
 550         double oppstatus = 0.0;
 551
 552     Bool_t fFlagPhotonicElec = kFALSE;
 553     Bool_t fFlagConvinatElec = kFALSE;
 554
 555     Int_t clsId = track->GetEMCALcluster();
 556     if (clsId>0){
 557       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 558       if(clust && clust->IsEMCAL()){
 559
 560                 double clustE = clust->E();
 561                 eop = clustE/fabs(mom);
 562                 //double clustT = clust->GetTOF();
 563                 ncells = clust->GetNCells();
 564                 m02 = clust->GetM02();
 565                 m20 = clust->GetM20();
 566                 disp = clust->GetDispersion();
 567                 double delphi = clust->GetTrackDx();
 568                 double deleta = clust->GetTrackDz();
 569                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 570                 nmatch = clust->GetNTracksMatched();
 571
 572                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
 573                 {
 574                   SelectPhotonicElectron(iTracks,cent,track,stack,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing);
 575                 }
 576                 if(fFlagPhotonicElec)oppstatus = 1.0;
 577                 if(fFlagConvinatElec)oppstatus = 2.0;
 578                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
 579
 580                   double valdedx[16];
 581                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 582                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
 583                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
 584                   valdedx[15] = mcele;
 585                   if(fqahist==1)fEleInfo->Fill(valdedx);
 586
 587
 588       }
 589     }
 590
 591     if(nITS<2.5)continue;
 592     if(nTPCcl<100)continue;
 593
 594     CheckNclust->Fill(nTPCcl);
 595     CheckNits->Fill(nITS);
 596
 597     fdEdxBef->Fill(mom,fTPCnSigma);
 598     fTPCnsigma->Fill(mom,fTPCnSigma);
 599     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 600
 601     Int_t pidpassed = 1;
 602
 603
 604     //--- track accepted
 605     AliHFEpidObject hfetrack;
 606     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 607     hfetrack.SetRecTrack(track);
 608     double binct = 10.5;
 609     if((0.0< cent) && (cent<5.0)) binct = 0.5;
 610     if((5.0< cent) && (cent<10.0)) binct = 1.5;
 611     if((10.0< cent) && (cent<20.0)) binct = 2.5;
 612     if((20.0< cent) && (cent<30.0)) binct = 3.5;
 613     if((30.0< cent) && (cent<40.0)) binct = 4.5;
 614     if((40.0< cent) && (cent<50.0)) binct = 5.5;
 615     if((50.0< cent) && (cent<60.0)) binct = 6.5;
 616     if((60.0< cent) && (cent<70.0)) binct = 7.5;
 617     if((70.0< cent) && (cent<80.0)) binct = 8.5;
 618     if((80.0< cent) && (cent<90.0)) binct = 9.5;
 619     if((90.0< cent) && (cent<100.0)) binct = 10.5;
 620
 621     hfetrack.SetCentrality((int)binct); //added
 622     hfetrack.SetPbPb();
 623     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 624
 625     if(pidpassed==0) continue;
 626
 627     fTrkEovPAft->Fill(pt,eop);
 628     fdEdxAft->Fill(mom,fTPCnSigma);
 629
 630     // Fill real data
 631     fIncpT->Fill(cent,pt);
 632     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 633     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 634
 635     if(m20>0.0 && m20<0.3)
 636       {
 637        fIncpTM20->Fill(cent,pt);
 638        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
 639        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
 640      }
 641
 642     // MC
 643     if(mcele>0) // select MC electrons
 644       {
 645
 646           fIncpTMChfeAll->Fill(cent,pt);
 647           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
 648
 649        if(mcBtoE || mcDtoE) // select B->e & D->e
 650          {
 651           fIncpTMChfe->Fill(cent,pt);
 652           if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
 653          }
 654
 655        if(mcPho) // select photonic electrons
 656         {
 657          double phoval[5];
 658          phoval[0] = cent;
 659          phoval[1] = pt;
 660          phoval[2] = fTPCnSigma;
 661          phoval[3] = iHijing;
 662          phoval[4] = mcMompT;
 663
 664          fIncpTMCpho->Fill(phoval);
 665          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
 666          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
 667
 668          if(m20>0.0 && m20<0.3)
 669            {
 670             fIncpTMCM20pho->Fill(phoval);
 671             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
 672             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
 673             // pi0->g->e
 674             if(mcWeight>-1)
 675               {
 676                if(iHijing==1)mcWeight = 1.0;
 677                fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
 678                if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
 679                if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
 680               }
 681            }
 682         }
 683       }
 684  }
 685  PostData(1, fOutputList);
 686 }
 687 //_________________________________________
 688 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
 689 {
 690   //--- Check MC
 691
 692   //Bool_t mcData = kFALSE;
 693   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 694     {
 695      fmcData = kTRUE;
 696      printf("+++++ MC Data available");
 697     }
 698   if(fmcData)
 699     {
 700      printf("++++++++= MC analysis \n");
 701     }
 702   else
 703    {
 704      printf("++++++++ real data analysis \n");
 705    }
 706
 707    printf("+++++++ QA hist %d",fqahist);
 708
 709   //---- Geometry
 710   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 711
 712   //--------Initialize PID
 713   //fPID->SetHasMCData(kFALSE);
 714   fPID->SetHasMCData(fmcData);
 715   if(!fPID->GetNumberOfPIDdetectors())
 716     {
 717       fPID->AddDetector("TPC", 0);
 718       fPID->AddDetector("EMCAL", 1);
 719     }
 720
 721   Double_t params[4];
 722   const char *cutmodel;
 723   cutmodel = "pol0";
 724   params[0] = -1.0; //sigma min
 725   double maxnSig = 3.0;
 726   if(fmcData)
 727     {
 728      params[0] = -5.0; //sigma min
 729      maxnSig = 5.0;
 730     }
 731
 732   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
 733
 734   fPID->SortDetectors();
 735   fPIDqa = new AliHFEpidQAmanager();
 736   fPIDqa->Initialize(fPID);
 737
 738   //------- fcut --------------
 739   fCuts = new AliHFEcuts();
 740   fCuts->CreateStandardCuts();
 741   //fCuts->SetMinNClustersTPC(100);
 742   fCuts->SetMinNClustersTPC(90);
 743   fCuts->SetMinRatioTPCclusters(0.6);
 744   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 745   //fCuts->SetMinNClustersITS(3);
 746   fCuts->SetMinNClustersITS(2);
 747   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
 748   fCuts->SetCheckITSLayerStatus(kFALSE);
 749   fCuts->SetVertexRange(10.);
 750   fCuts->SetTOFPIDStep(kFALSE);
 751   fCuts->SetPtRange(2, 50);
 752   fCuts->SetMaxImpactParam(3.,3.);
 753
 754   //--------Initialize correction Framework and Cuts
 755   fCFM = new AliCFManager;
 756   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 757   fCFM->SetNStepParticle(kNcutSteps);
 758   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 759     fCFM->SetParticleCutsList(istep, NULL);
 760
 761   if(!fCuts){
 762     AliWarning("Cuts not available. Default cuts will be used");
 763     fCuts = new AliHFEcuts;
 764     fCuts->CreateStandardCuts();
 765   }
 766   fCuts->Initialize(fCFM);
 767
 768   //---------Output Tlist
 769   fOutputList = new TList();
 770   fOutputList->SetOwner();
 771   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 772
 773   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
 774   fOutputList->Add(fNoEvents);
 775
 776   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
 777   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
 778   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
 779   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
 780   if(fqahist==1)fOutputList->Add(fEMCAccE);
 781
 782   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 783   fOutputList->Add(fTrkpt);
 784
 785   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 786   fOutputList->Add(fTrackPtBefTrkCuts);
 787
 788   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 789   fOutputList->Add(fTrackPtAftTrkCuts);
 790
 791   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 792   fOutputList->Add(fTPCnsigma);
 793
 794   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 795   fOutputList->Add(fTrkEovPBef);
 796
 797   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 798   fOutputList->Add(fTrkEovPAft);
 799
 800   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
 801   fOutputList->Add(fdEdxBef);
 802
 803   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
 804   fOutputList->Add(fdEdxAft);
 805
 806   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
 807   fOutputList->Add(fIncpT);
 808
 809   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
 810   fOutputList->Add(fIncpTM20);
 811
 812   Int_t nBinspho[9] =  { 200, 100, 500, 12,   50,    4, 200,   8, 100};
 813   Double_t minpho[9] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5,   0};
 814   Double_t maxpho[9] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5,  50};
 815
 816   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);
 817   fOutputList->Add(fInvmassLS);
 818
 819   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);
 820   fOutputList->Add(fInvmassULS);
 821
 822   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);
 823   fOutputList->Add(fInvmassLSmc);
 824
 825   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);
 826   fOutputList->Add(fInvmassULSmc);
 827
 828   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 829   fOutputList->Add(fOpeningAngleLS);
 830
 831   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 832   fOutputList->Add(fOpeningAngleULS);
 833
 834   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 835   fOutputList->Add(fPhotoElecPt);
 836
 837   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
 838   fOutputList->Add(fPhoElecPt);
 839
 840   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
 841   fOutputList->Add(fPhoElecPtM20);
 842
 843   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
 844   fOutputList->Add(fSameElecPt);
 845
 846   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
 847   fOutputList->Add(fSameElecPtM20);
 848
 849   fCent = new TH1F("fCent","Centrality",200,0,100) ;
 850   fOutputList->Add(fCent);
 851
 852   // Make common binning
 853   const Double_t kMinP = 0.;
 854   const Double_t kMaxP = 50.;
 855   //const Double_t kTPCSigMim = 40.;
 856   //const Double_t kTPCSigMax = 140.;
 857
 858   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
 859   Int_t nBins[16] =  {  250,    10,   60,    20,   100,  300,  50,   40,   200, 200,  250, 200,    3,    5, 100,    8};
 860   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};
 861   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};
 862   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);
 863   if(fqahist==1)fOutputList->Add(fEleInfo);
 864
 865   //<---  Trigger info
 866   /*
 867   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
 868   fOutputList->Add(fClsEBftTrigCut);
 869
 870   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
 871   fOutputList->Add(fClsEAftTrigCut);
 872
 873   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
 874   fOutputList->Add(fClsEAftTrigCut1);
 875
 876   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
 877   fOutputList->Add(fClsEAftTrigCut2);
 878
 879   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
 880   fOutputList->Add(fClsEAftTrigCut3);
 881
 882   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
 883   fOutputList->Add(fClsEAftTrigCut4);
 884
 885   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 886   fOutputList->Add(fClsETime);
 887
 888   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 889   fOutputList->Add(fClsETime1);
 890
 891   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
 892   fOutputList->Add(fTrigTimes);
 893
 894   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
 895   fOutputList->Add(fCellCheck);
 896   */
 897   //<---------- MC
 898
 899   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
 900   fOutputList->Add(fInputHFEMC);
 901
 902   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
 903   fOutputList->Add(fInputAlle);
 904
 905   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
 906   fOutputList->Add(fIncpTMChfe);
 907
 908   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
 909   fOutputList->Add(fIncpTMChfeAll);
 910
 911   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
 912   fOutputList->Add(fIncpTMCM20hfe);
 913
 914   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
 915   fOutputList->Add(fIncpTMCM20hfeAll);
 916
 917
 918   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 100};
 919   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};
 920   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 50.};
 921
 922   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
 923   fOutputList->Add(fIncpTMCpho);
 924
 925   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
 926   fOutputList->Add(fIncpTMCM20pho);
 927
 928   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
 929   fOutputList->Add(fPhoElecPtMC);
 930
 931   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
 932   fOutputList->Add(fPhoElecPtMCM20);
 933
 934   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
 935   fOutputList->Add(fSameElecPtMC);
 936
 937   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
 938   fOutputList->Add(fSameElecPtMCM20);
 939
 940   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
 941   fIncpTMCM20pho_pi0e->Sumw2();
 942   fOutputList->Add(fIncpTMCM20pho_pi0e);
 943
 944   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
 945   fPhoElecPtMCM20_pi0e->Sumw2();
 946   fOutputList->Add(fPhoElecPtMCM20_pi0e);
 947
 948   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
 949   fSameElecPtMCM20_pi0e->Sumw2();
 950   fOutputList->Add(fSameElecPtMCM20_pi0e);
 951
 952   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
 953   fOutputList->Add(CheckNclust);
 954
 955   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
 956   fOutputList->Add(CheckNits);
 957
 958   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
 959   fOutputList->Add(Hpi0pTcheck);
 960
 961   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
 962   fOutputList->Add(HphopTcheck);
 963
 964   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
 965   fOutputList->Add(fMomDtoE);
 966
 967   PostData(1,fOutputList);
 968 }
 969
 970 //________________________________________________________________________
 971 void AliAnalysisTaskHFECal::Terminate(Option_t *)
 972 {
 973   // Info("Terminate");
 974         AliAnalysisTaskSE::Terminate();
 975 }
 976
 977 //________________________________________________________________________
 978 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 979 {
 980   // Check single track cuts for a given cut step
 981   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 982   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 983   return kTRUE;
 984 }
 985 //_________________________________________
 986 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
 987 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
 988 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, AliStack *stack ,Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent)
 989 {
 990   //Identify non-heavy flavour electrons using Invariant mass method
 991
 992   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 993   fTrackCuts->SetRequireTPCRefit(kTRUE);
 994   fTrackCuts->SetRequireITSRefit(kTRUE);
 995   fTrackCuts->SetEtaRange(-0.9,0.9);
 996   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 997   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 998   fTrackCuts->SetMinNClustersTPC(90);
 999
1000   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1001
1002   Bool_t flagPhotonicElec = kFALSE;
1003   Bool_t flagConvinatElec = kFALSE;
1004
1005   int p1 = 0;
1006   if(mce==3)
1007      {
1008        Int_t label = TMath::Abs(track->GetLabel());
1009        TParticle* particle = stack->Particle(label);
1010        p1 = particle->GetFirstMother();
1011      }
1012
1013   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1014   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1015     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1016     if (!trackAsso) {
1017       printf("ERROR: Could not receive track %d\n", jTracks);
1018       continue;
1019     }
1020     if(itrack==jTracks)continue;
1021
1022     int p2 = 0;
1023     if(mce==3)
1024     {
1025       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1026       TParticle* particle2 = stack->Particle(label2);
1027       if(particle2->GetFirstMother()>-1)
1028          p2 = particle2->GetFirstMother();
1029     }
1030
1031     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1032     Double_t mass=999., width = -999;
1033     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1034
1035     ptPrim = track->Pt();
1036
1037     dEdxAsso = trackAsso->GetTPCsignal();
1038     ptAsso = trackAsso->Pt();
1039     Int_t chargeAsso = trackAsso->Charge();
1040     Int_t charge = track->Charge();
1041
1042
1043     if(ptAsso <0.5) continue;
1044     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1045     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1046
1047     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1048     if(charge>0) fPDGe1 = -11;
1049     if(chargeAsso>0) fPDGe2 = -11;
1050
1051     //printf("chargeAsso = %d\n",chargeAsso);
1052     //printf("charge = %d\n",charge);
1053     if(charge == chargeAsso) fFlagLS = kTRUE;
1054     if(charge != chargeAsso) fFlagULS = kTRUE;
1055
1056     //printf("fFlagLS = %d\n",fFlagLS);
1057     //printf("fFlagULS = %d\n",fFlagULS);
1058     //printf("\n");
1059
1060     AliKFParticle ge1(*track, fPDGe1);
1061     AliKFParticle ge2(*trackAsso, fPDGe2);
1062     AliKFParticle recg(ge1, ge2);
1063
1064     if(recg.GetNDF()<1) continue;
1065     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1066     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1067
1068     AliKFVertex primV(*pVtx);
1069     primV += recg;
1070     recg.SetProductionVertex(primV);
1071
1072     recg.SetMassConstraint(0,0.0001);
1073
1074     openingAngle = ge1.GetAngle(ge2);
1075     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1076     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1077
1078
1079     recg.GetMass(mass,width);
1080
1081     double ishower = 0;
1082     if(shower>0.0 && shower<0.3)ishower = 1;
1083
1084     double phoinfo[9];
1085     phoinfo[0] = cent;
1086     phoinfo[1] = ptPrim;
1087     phoinfo[2] = mass;
1088     phoinfo[3] = nSig;
1089     phoinfo[4] = openingAngle;
1090     phoinfo[5] = ishower;
1091     phoinfo[6] = ep;
1092     phoinfo[7] = mce;
1093     phoinfo[8] = ptAsso;
1094
1095     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1096     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1097     if(fFlagLS && ibgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1098     if(fFlagULS && ibgevent==0) fInvmassULSmc->Fill(phoinfo,w);
1099
1100     //printf("fInvmassCut %f\n",fInvmassCut);
1101     //printf("openingAngle %f\n",fOpeningAngleCut);
1102
1103     if(openingAngle > fOpeningAngleCut) continue;
1104
1105     //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1106     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){
1107       flagPhotonicElec = kTRUE;
1108     }
1109     if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1110       flagConvinatElec = kTRUE;
1111     }
1112
1113   }
1114   fFlagPhotonicElec = flagPhotonicElec;
1115   fFlagConvinatElec = flagConvinatElec;
1116
1117 }
1118 //-------------------------------------------
1119 //void AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1120 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1121 {
1122         double weight = 1.0;
1123
1124         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1125         {
1126                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1127         }
1128         else
1129         {
1130                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1131         }
1132   return weight;
1133 }
1134
1135 //_________________________________________
1136 void AliAnalysisTaskHFECal::FindTriggerClusters()
1137 {
1138   // constants
1139   const int nModuleCols = 2;
1140   const int nModuleRows = 5;
1141   const int nColsFeeModule = 48;
1142   const int nRowsFeeModule = 24;
1143   const int nColsFaltroModule = 24;
1144   const int nRowsFaltroModule = 12;
1145   //const int faltroWidthMax = 20;
1146
1147   // part 1, trigger extraction -------------------------------------
1148   Int_t globCol, globRow;
1149   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1150   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1151
1152   //Int_t trigtimes[faltroWidthMax];
1153   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1154   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1155   //Double_t fTrigCutLow = 6;
1156   //Double_t fTrigCutHigh = 10;
1157   Double_t fTimeCutLow =  469e-09;
1158   Double_t fTimeCutHigh = 715e-09;
1159
1160   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1161
1162   // erase trigger maps
1163   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1164   {
1165     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1166     {
1167       ftriggersCut[i][j] = 0;
1168       ftriggers[i][j] = 0;
1169       ftriggersTime[i][j] = 0;
1170     }
1171   }
1172
1173   Int_t iglobCol=0, iglobRow=0;
1174   // go through triggers
1175   if( fCaloTrigger->GetEntries() > 0 )
1176   {
1177     // needs reset
1178     fCaloTrigger->Reset();
1179     while( fCaloTrigger->Next() )
1180     {
1181       fCaloTrigger->GetPosition( globCol, globRow );
1182       fCaloTrigger->GetNL0Times( ntimes );
1183       /*
1184       // no L0s
1185       if( ntimes < 1 )   continue;
1186       // get precise timings
1187       fCaloTrigger->GetL0Times( trigtimes );
1188       trigInCut = 0;
1189       for(Int_t i = 0; i < ntimes; i++ )
1190       {
1191         // save the first trigger time in channel
1192         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1193           triggersTime[globCol][globRow] = trigtimes[i];
1194         //printf("trigger times: %d\n",trigtimes[i]);
1195         // check if in cut
1196         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1197           trigInCut = 1;
1198
1199         fTrigTimes->Fill(trigtimes[i]);
1200       }
1201       */
1202
1203       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1204       Int_t bit = 0;
1205       fCaloTrigger->GetTriggerBits(bit);
1206
1207       Int_t ts = 0;
1208       fCaloTrigger->GetL1TimeSum(ts);
1209       if (ts > 0)ftriggers[globCol][globRow] = 1;
1210       // number of triggered channels in event
1211       nTrigChannel++;
1212       // ... inside cut
1213       if(ts>0 && (bit >> 6 & 0x1))
1214       {
1215         iglobCol = globCol;
1216         iglobRow = globRow;
1217         nTrigChannelCut++;
1218         ftriggersCut[globCol][globRow] = 1;
1219       }
1220
1221     } // calo trigger entries
1222   } // has calo trigger entries
1223
1224   // part 2 go through the clusters here -----------------------------------
1225   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1226   Short_t cellAddr, nSACell, mclabel;
1227   //Int_t nSACell, iSACell, mclabel;
1228   Int_t iSACell;
1229   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1230   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1231
1232   TRefArray *fCaloClusters = new TRefArray();
1233   fESD->GetEMCALClusters( fCaloClusters );
1234   nCluster = fCaloClusters->GetEntries();
1235
1236
1237   // save all cells times if there are clusters
1238   if( nCluster > 0 ){
1239     // erase time map
1240     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1241       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1242         cellTime[i][j] = 0.;
1243         cellEnergy[i][j] = 0.;
1244       }
1245     }
1246
1247     // get cells
1248     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1249     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1250     nSACell = fCaloCells->GetNumberOfCells();
1251     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1252       // get the cell info *fCal
1253       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1254
1255       // get cell position
1256       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1257       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1258
1259       // convert co global phi eta
1260       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1261       geta = ieta + nColsFeeModule*(nSupMod%2);
1262
1263       // save cell time and energy
1264       cellTime[geta][gphi] = cellTimeT;
1265       cellEnergy[geta][gphi] = cellAmp;
1266
1267     }
1268   }
1269
1270   Int_t nClusterTrig, nClusterTrigCut;
1271   UShort_t *cellAddrs;
1272   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1273   Float_t clsPos[3] = {0.,0.,0.};
1274
1275   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1276   {
1277     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1278     if(!cluster || !cluster->IsEMCAL()) continue;
1279
1280     // get cluster cells
1281     nCell = cluster->GetNCells();
1282
1283     // get cluster energy
1284     clsE = cluster->E();
1285
1286     // get cluster position
1287     cluster->GetPosition(clsPos);
1288     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1289     clsEta = clsPosVec.Eta();
1290     clsPhi = clsPosVec.Phi();
1291
1292     // get the cell addresses
1293     cellAddrs = cluster->GetCellsAbsId();
1294
1295     // check if the cluster contains cell, that was marked as triggered
1296     nClusterTrig = 0;
1297     nClusterTrigCut = 0;
1298
1299     // loop the cells to check, if cluser in acceptance
1300     // any cluster with a cell outside acceptance is not considered
1301     for( iCell = 0; iCell < nCell; iCell++ )
1302     {
1303      // check hot cell
1304      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
1305
1306       // get cell position
1307       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1308       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1309
1310       // convert co global phi eta
1311       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1312       geta = ieta + nColsFeeModule*(nSupMod%2);
1313
1314       if( cellTime[geta][gphi] > 0. ){
1315         clusterTime += cellTime[geta][gphi];
1316         gCell++;
1317       }
1318
1319       // get corresponding FALTRO
1320       fphi = gphi / 2;
1321       feta = geta / 2;
1322
1323       // try to match with a triggered
1324       if( ftriggers[feta][fphi]==1)
1325       {  nClusterTrig++;
1326       }
1327       if( ftriggersCut[feta][fphi]==1)
1328       { nClusterTrigCut++;
1329       }
1330
1331     } // cells
1332
1333
1334     if( gCell > 0 )
1335       clusterTime = clusterTime / (Double_t)gCell;
1336     // fix the reconstriction code time 100ns jumps
1337     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1338       clusterTime -= 0.0000001;
1339
1340     //fClsETime->Fill(clsE,clusterTime);
1341     //fClsEBftTrigCut->Fill(clsE);
1342
1343     if(nClusterTrig>0){
1344       //fClsETime1->Fill(clsE,clusterTime);
1345     }
1346
1347     if(nClusterTrig>0){
1348       cluster->SetChi2(1);
1349       //fClsEAftTrigCut1->Fill(clsE);
1350     }
1351
1352     if(nClusterTrigCut>0){
1353       cluster->SetChi2(2);
1354       //fClsEAftTrigCut2->Fill(clsE);
1355     }
1356
1357     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1358     {
1359       cluster->SetChi2(3);
1360       //fClsEAftTrigCut3->Fill(clsE);
1361     }
1362
1363     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1364     {
1365       // cluster->SetChi2(4);
1366       //fClsEAftTrigCut4->Fill(clsE);
1367     }
1368     if(nClusterTrigCut<1)
1369     {
1370       cluster->SetChi2(0);
1371
1372       //fClsEAftTrigCut->Fill(clsE);
1373     }
1374
1375   } // clusters
1376 }
1377
1378
1379
1380
1381
1382
1383