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