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