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 nTPCclF = track->GetTPCNclsF();
 440     int nITS = track->GetNcls(0);
 441
 442     fTrackPtAftTrkCuts->Fill(track->Pt());
 443
 444     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 445     pt = track->Pt();
 446     if(pt<2.0)continue;
 447
 448     // Track extrapolation
 449
 450     Int_t charge = track->Charge();
 451     fTrkpt->Fill(pt);
 452     mom = track->P();
 453     phi = track->Phi();
 454     eta = track->Eta();
 455     dEdx = track->GetTPCsignal();
 456     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 457
 458         double ncells = -1.0;
 459         double m20 = -1.0;
 460         double m02 = -1.0;
 461         double disp = -1.0;
 462         double rmatch = -1.0;
 463         double nmatch = -1.0;
 464         double oppstatus = 0.0;
 465
 466     Bool_t fFlagPhotonicElec = kFALSE;
 467     Bool_t fFlagConvinatElec = kFALSE;
 468
 469     Int_t clsId = track->GetEMCALcluster();
 470     if (clsId>0){
 471       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 472       if(clust && clust->IsEMCAL()){
 473
 474                 double clustE = clust->E();
 475                 eop = clustE/fabs(mom);
 476                 //double clustT = clust->GetTOF();
 477                 ncells = clust->GetNCells();
 478                 m02 = clust->GetM02();
 479                 m20 = clust->GetM20();
 480                 disp = clust->GetDispersion();
 481                 double delphi = clust->GetTrackDx();
 482                 double deleta = clust->GetTrackDz();
 483                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 484                 nmatch = clust->GetNTracksMatched();
 485
 486                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
 487                 {
 488                   SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele);
 489                 }
 490                 if(fFlagPhotonicElec)oppstatus = 1.0;
 491                 if(fFlagConvinatElec)oppstatus = 2.0;
 492                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
 493
 494                   double valdedx[16];
 495                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 496                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
 497                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
 498                   valdedx[15] = mcele;
 499                   if(fqahist==1)fEleInfo->Fill(valdedx);
 500
 501
 502       }
 503     }
 504
 505     if(nITS<2.5)continue;
 506     if(nTPCcl<100)continue;
 507
 508     CheckNclust->Fill(nTPCcl);
 509     CheckNits->Fill(nITS);
 510
 511     fdEdxBef->Fill(mom,fTPCnSigma);
 512     fTPCnsigma->Fill(mom,fTPCnSigma);
 513     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 514
 515     Int_t pidpassed = 1;
 516
 517
 518     //--- track accepted
 519     AliHFEpidObject hfetrack;
 520     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 521     hfetrack.SetRecTrack(track);
 522     double binct = 10.5;
 523     if((0.0< cent) && (cent<5.0)) binct = 0.5;
 524     if((5.0< cent) && (cent<10.0)) binct = 1.5;
 525     if((10.0< cent) && (cent<20.0)) binct = 2.5;
 526     if((20.0< cent) && (cent<30.0)) binct = 3.5;
 527     if((30.0< cent) && (cent<40.0)) binct = 4.5;
 528     if((40.0< cent) && (cent<50.0)) binct = 5.5;
 529     if((50.0< cent) && (cent<60.0)) binct = 6.5;
 530     if((60.0< cent) && (cent<70.0)) binct = 7.5;
 531     if((70.0< cent) && (cent<80.0)) binct = 8.5;
 532     if((80.0< cent) && (cent<90.0)) binct = 9.5;
 533     if((90.0< cent) && (cent<100.0)) binct = 10.5;
 534
 535     hfetrack.SetCentrality((int)binct); //added
 536     hfetrack.SetPbPb();
 537     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 538
 539     if(pidpassed==0) continue;
 540
 541     fTrkEovPAft->Fill(pt,eop);
 542     fdEdxAft->Fill(mom,fTPCnSigma);
 543
 544     // Fill real data
 545     fIncpT->Fill(cent,pt);
 546     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 547     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 548
 549     if(m20>0.0 && m20<0.3)
 550       {
 551        fIncpTM20->Fill(cent,pt);
 552        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
 553        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
 554      }
 555
 556     // MC
 557     if(mcele>0)
 558       {
 559
 560           fIncpTMChfeAll->Fill(cent,pt);
 561           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
 562
 563        if(mcBtoE || mcDtoE)
 564          {
 565           fIncpTMChfe->Fill(cent,pt);
 566           if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
 567          }
 568        if(mcPho)
 569         {
 570          double phoval[3];
 571          phoval[0] = cent;
 572          phoval[1] = pt;
 573          phoval[2] = fTPCnSigma;
 574
 575          fIncpTMCpho->Fill(phoval);
 576          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
 577          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
 578
 579          if(m20>0.0 && m20<0.3)
 580            {
 581             fIncpTMCM20pho->Fill(phoval);
 582             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
 583             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
 584            }
 585         }
 586       }
 587  }
 588  PostData(1, fOutputList);
 589 }
 590 //_________________________________________
 591 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
 592 {
 593   //--- Check MC
 594
 595   //Bool_t mcData = kFALSE;
 596   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 597     {
 598      fmcData = kTRUE;
 599      printf("+++++ MC Data available");
 600     }
 601   if(fmcData)
 602     {
 603      printf("++++++++= MC analysis \n");
 604     }
 605   else
 606    {
 607      printf("++++++++ real data analysis \n");
 608    }
 609
 610    printf("+++++++ QA hist %d",fqahist);
 611
 612   //---- Geometry
 613   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 614
 615   //--------Initialize PID
 616   //fPID->SetHasMCData(kFALSE);
 617   fPID->SetHasMCData(fmcData);
 618   if(!fPID->GetNumberOfPIDdetectors())
 619     {
 620       fPID->AddDetector("TPC", 0);
 621       fPID->AddDetector("EMCAL", 1);
 622     }
 623
 624   Double_t params[4];
 625   const char *cutmodel;
 626   cutmodel = "pol0";
 627   params[0] = -1.0; //sigma min
 628   double maxnSig = 3.0;
 629   if(fmcData)
 630     {
 631      params[0] = -5.0; //sigma min
 632      maxnSig = 5.0;
 633     }
 634
 635   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
 636
 637   fPID->SortDetectors();
 638   fPIDqa = new AliHFEpidQAmanager();
 639   fPIDqa->Initialize(fPID);
 640
 641   //------- fcut --------------
 642   fCuts = new AliHFEcuts();
 643   fCuts->CreateStandardCuts();
 644   //fCuts->SetMinNClustersTPC(100);
 645   fCuts->SetMinNClustersTPC(90);
 646   fCuts->SetMinRatioTPCclusters(0.6);
 647   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 648   //fCuts->SetMinNClustersITS(3);
 649   fCuts->SetMinNClustersITS(2);
 650   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
 651   fCuts->SetCheckITSLayerStatus(kFALSE);
 652   fCuts->SetVertexRange(10.);
 653   fCuts->SetTOFPIDStep(kFALSE);
 654   fCuts->SetPtRange(2, 50);
 655   fCuts->SetMaxImpactParam(3.,3.);
 656
 657   //--------Initialize correction Framework and Cuts
 658   fCFM = new AliCFManager;
 659   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 660   fCFM->SetNStepParticle(kNcutSteps);
 661   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 662     fCFM->SetParticleCutsList(istep, NULL);
 663
 664   if(!fCuts){
 665     AliWarning("Cuts not available. Default cuts will be used");
 666     fCuts = new AliHFEcuts;
 667     fCuts->CreateStandardCuts();
 668   }
 669   fCuts->Initialize(fCFM);
 670
 671   //---------Output Tlist
 672   fOutputList = new TList();
 673   fOutputList->SetOwner();
 674   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 675
 676   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
 677   fOutputList->Add(fNoEvents);
 678
 679   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
 680   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
 681   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
 682   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
 683   if(fqahist==1)fOutputList->Add(fEMCAccE);
 684
 685   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 686   fOutputList->Add(fTrkpt);
 687
 688   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 689   fOutputList->Add(fTrackPtBefTrkCuts);
 690
 691   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 692   fOutputList->Add(fTrackPtAftTrkCuts);
 693
 694   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 695   fOutputList->Add(fTPCnsigma);
 696
 697   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 698   fOutputList->Add(fTrkEovPBef);
 699
 700   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 701   fOutputList->Add(fTrkEovPAft);
 702
 703   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
 704   fOutputList->Add(fdEdxBef);
 705
 706   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
 707   fOutputList->Add(fdEdxAft);
 708
 709   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
 710   fOutputList->Add(fIncpT);
 711
 712   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
 713   fOutputList->Add(fIncpTM20);
 714
 715   Int_t nBinspho[9] =  { 200, 100, 500, 12,   50,    4, 200,   8, 100};
 716   Double_t minpho[9] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5,   0};
 717   Double_t maxpho[9] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5,  50};
 718
 719   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);
 720   fOutputList->Add(fInvmassLS);
 721
 722   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);
 723   fOutputList->Add(fInvmassULS);
 724
 725   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 726   fOutputList->Add(fOpeningAngleLS);
 727
 728   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 729   fOutputList->Add(fOpeningAngleULS);
 730
 731   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 732   fOutputList->Add(fPhotoElecPt);
 733
 734   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
 735   fOutputList->Add(fPhoElecPt);
 736
 737   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
 738   fOutputList->Add(fPhoElecPtM20);
 739
 740   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
 741   fOutputList->Add(fSameElecPt);
 742
 743   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
 744   fOutputList->Add(fSameElecPtM20);
 745
 746   fCent = new TH1F("fCent","Centrality",200,0,100) ;
 747   fOutputList->Add(fCent);
 748
 749   // Make common binning
 750   const Double_t kMinP = 0.;
 751   const Double_t kMaxP = 50.;
 752   //const Double_t kTPCSigMim = 40.;
 753   //const Double_t kTPCSigMax = 140.;
 754
 755   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
 756   Int_t nBins[16] =  {  250,    10,   60,    20,   100,  300,  50,   40,   200, 200,  250, 200,    3,    5, 100,    8};
 757   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};
 758   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};
 759   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);
 760   if(fqahist==1)fOutputList->Add(fEleInfo);
 761
 762   //<---  Trigger info
 763   /*
 764   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
 765   fOutputList->Add(fClsEBftTrigCut);
 766
 767   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
 768   fOutputList->Add(fClsEAftTrigCut);
 769
 770   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
 771   fOutputList->Add(fClsEAftTrigCut1);
 772
 773   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
 774   fOutputList->Add(fClsEAftTrigCut2);
 775
 776   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
 777   fOutputList->Add(fClsEAftTrigCut3);
 778
 779   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
 780   fOutputList->Add(fClsEAftTrigCut4);
 781
 782   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 783   fOutputList->Add(fClsETime);
 784
 785   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 786   fOutputList->Add(fClsETime1);
 787
 788   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
 789   fOutputList->Add(fTrigTimes);
 790
 791   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
 792   fOutputList->Add(fCellCheck);
 793   */
 794   //<---------- MC
 795
 796   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
 797   fOutputList->Add(fInputHFEMC);
 798
 799   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
 800   fOutputList->Add(fInputAlle);
 801
 802   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
 803   fOutputList->Add(fIncpTMChfe);
 804
 805   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
 806   fOutputList->Add(fIncpTMChfeAll);
 807
 808   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
 809   fOutputList->Add(fIncpTMCM20hfe);
 810
 811   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
 812   fOutputList->Add(fIncpTMCM20hfeAll);
 813
 814
 815   Int_t nBinspho2[3] =  { 200, 100,    7};
 816   Double_t minpho2[3] = {  0.,  0., -2.5};
 817   Double_t maxpho2[3] = {100., 50.,  4.5};
 818
 819   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",3,nBinspho2,minpho2,maxpho2);
 820   fOutputList->Add(fIncpTMCpho);
 821
 822   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",3,nBinspho2,minpho2,maxpho2);
 823   fOutputList->Add(fIncpTMCM20pho);
 824
 825   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",3,nBinspho2,minpho2,maxpho2);
 826   fOutputList->Add(fPhoElecPtMC);
 827
 828   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",3,nBinspho2,minpho2,maxpho2);
 829   fOutputList->Add(fPhoElecPtMCM20);
 830
 831   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",3,nBinspho2,minpho2,maxpho2);
 832   fOutputList->Add(fSameElecPtMC);
 833
 834   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",3,nBinspho2,minpho2,maxpho2);
 835   fOutputList->Add(fSameElecPtMCM20);
 836
 837   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
 838   fOutputList->Add(CheckNclust);
 839
 840   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
 841   fOutputList->Add(CheckNits);
 842
 843   PostData(1,fOutputList);
 844 }
 845
 846 //________________________________________________________________________
 847 void AliAnalysisTaskHFECal::Terminate(Option_t *)
 848 {
 849   // Info("Terminate");
 850         AliAnalysisTaskSE::Terminate();
 851 }
 852
 853 //________________________________________________________________________
 854 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 855 {
 856   // Check single track cuts for a given cut step
 857   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 858   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 859   return kTRUE;
 860 }
 861 //_________________________________________
 862 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
 863 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
 864 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)
 865 {
 866   //Identify non-heavy flavour electrons using Invariant mass method
 867
 868   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 869   fTrackCuts->SetRequireTPCRefit(kTRUE);
 870   fTrackCuts->SetRequireITSRefit(kTRUE);
 871   fTrackCuts->SetEtaRange(-0.9,0.9);
 872   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 873   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 874   fTrackCuts->SetMinNClustersTPC(90);
 875
 876   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 877
 878   Bool_t flagPhotonicElec = kFALSE;
 879   Bool_t flagConvinatElec = kFALSE;
 880
 881   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 882   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 883     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 884     if (!trackAsso) {
 885       printf("ERROR: Could not receive track %d\n", jTracks);
 886       continue;
 887     }
 888     if(itrack==jTracks)continue;
 889
 890     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
 891     Double_t mass=999., width = -999;
 892     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 893
 894     ptPrim = track->Pt();
 895
 896     dEdxAsso = trackAsso->GetTPCsignal();
 897     ptAsso = trackAsso->Pt();
 898     Int_t chargeAsso = trackAsso->Charge();
 899     Int_t charge = track->Charge();
 900
 901
 902     if(ptAsso <0.5) continue;
 903     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 904     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
 905
 906     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 907     if(charge>0) fPDGe1 = -11;
 908     if(chargeAsso>0) fPDGe2 = -11;
 909
 910     //printf("chargeAsso = %d\n",chargeAsso);
 911     //printf("charge = %d\n",charge);
 912     if(charge == chargeAsso) fFlagLS = kTRUE;
 913     if(charge != chargeAsso) fFlagULS = kTRUE;
 914
 915     //printf("fFlagLS = %d\n",fFlagLS);
 916     //printf("fFlagULS = %d\n",fFlagULS);
 917     //printf("\n");
 918
 919     AliKFParticle ge1(*track, fPDGe1);
 920     AliKFParticle ge2(*trackAsso, fPDGe2);
 921     AliKFParticle recg(ge1, ge2);
 922
 923     if(recg.GetNDF()<1) continue;
 924     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 925     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 926
 927     AliKFVertex primV(*pVtx);
 928     primV += recg;
 929     recg.SetProductionVertex(primV);
 930
 931     recg.SetMassConstraint(0,0.0001);
 932
 933     openingAngle = ge1.GetAngle(ge2);
 934     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 935     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 936
 937
 938     recg.GetMass(mass,width);
 939
 940     double ishower = 0;
 941     if(shower>0.0 && shower<0.3)ishower = 1;
 942
 943     double phoinfo[9];
 944     phoinfo[0] = cent;
 945     phoinfo[1] = ptPrim;
 946     phoinfo[2] = mass;
 947     phoinfo[3] = nSig;
 948     phoinfo[4] = openingAngle;
 949     phoinfo[5] = ishower;
 950     phoinfo[6] = ep;
 951     phoinfo[7] = mce;
 952     phoinfo[8] = ptAsso;
 953
 954     if(fFlagLS) fInvmassLS->Fill(phoinfo);
 955     if(fFlagULS) fInvmassULS->Fill(phoinfo);
 956
 957     //printf("fInvmassCut %f\n",fInvmassCut);
 958     //printf("openingAngle %f\n",fOpeningAngleCut);
 959
 960     if(openingAngle > fOpeningAngleCut) continue;
 961
 962     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 963       flagPhotonicElec = kTRUE;
 964     }
 965     if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
 966       flagConvinatElec = kTRUE;
 967     }
 968
 969   }
 970   fFlagPhotonicElec = flagPhotonicElec;
 971   fFlagConvinatElec = flagConvinatElec;
 972
 973 }
 974
 975
 976 //_________________________________________
 977 void AliAnalysisTaskHFECal::FindTriggerClusters()
 978 {
 979   // constants
 980   const int nModuleCols = 2;
 981   const int nModuleRows = 5;
 982   const int nColsFeeModule = 48;
 983   const int nRowsFeeModule = 24;
 984   const int nColsFaltroModule = 24;
 985   const int nRowsFaltroModule = 12;
 986   //const int faltroWidthMax = 20;
 987
 988   // part 1, trigger extraction -------------------------------------
 989   Int_t globCol, globRow;
 990   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
 991   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
 992
 993   //Int_t trigtimes[faltroWidthMax];
 994   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
 995   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
 996   //Double_t fTrigCutLow = 6;
 997   //Double_t fTrigCutHigh = 10;
 998   Double_t fTimeCutLow =  469e-09;
 999   Double_t fTimeCutHigh = 715e-09;
1000
1001   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1002
1003   // erase trigger maps
1004   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1005   {
1006     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1007     {
1008       ftriggersCut[i][j] = 0;
1009       ftriggers[i][j] = 0;
1010       ftriggersTime[i][j] = 0;
1011     }
1012   }
1013
1014   Int_t iglobCol=0, iglobRow=0;
1015   // go through triggers
1016   if( fCaloTrigger->GetEntries() > 0 )
1017   {
1018     // needs reset
1019     fCaloTrigger->Reset();
1020     while( fCaloTrigger->Next() )
1021     {
1022       fCaloTrigger->GetPosition( globCol, globRow );
1023       fCaloTrigger->GetNL0Times( ntimes );
1024       /*
1025       // no L0s
1026       if( ntimes < 1 )   continue;
1027       // get precise timings
1028       fCaloTrigger->GetL0Times( trigtimes );
1029       trigInCut = 0;
1030       for(Int_t i = 0; i < ntimes; i++ )
1031       {
1032         // save the first trigger time in channel
1033         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1034           triggersTime[globCol][globRow] = trigtimes[i];
1035         //printf("trigger times: %d\n",trigtimes[i]);
1036         // check if in cut
1037         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1038           trigInCut = 1;
1039
1040         fTrigTimes->Fill(trigtimes[i]);
1041       }
1042       */
1043
1044       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1045       Int_t bit = 0;
1046       fCaloTrigger->GetTriggerBits(bit);
1047
1048       Int_t ts = 0;
1049       fCaloTrigger->GetL1TimeSum(ts);
1050       if (ts > 0)ftriggers[globCol][globRow] = 1;
1051       // number of triggered channels in event
1052       nTrigChannel++;
1053       // ... inside cut
1054       if(ts>0 && (bit >> 6 & 0x1))
1055       {
1056         iglobCol = globCol;
1057         iglobRow = globRow;
1058         nTrigChannelCut++;
1059         ftriggersCut[globCol][globRow] = 1;
1060       }
1061
1062     } // calo trigger entries
1063   } // has calo trigger entries
1064
1065   // part 2 go through the clusters here -----------------------------------
1066   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1067   Short_t cellAddr, nSACell, mclabel;
1068   //Int_t nSACell, iSACell, mclabel;
1069   Int_t iSACell;
1070   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1071   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1072
1073   TRefArray *fCaloClusters = new TRefArray();
1074   fESD->GetEMCALClusters( fCaloClusters );
1075   nCluster = fCaloClusters->GetEntries();
1076
1077
1078   // save all cells times if there are clusters
1079   if( nCluster > 0 ){
1080     // erase time map
1081     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1082       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1083         cellTime[i][j] = 0.;
1084         cellEnergy[i][j] = 0.;
1085       }
1086     }
1087
1088     // get cells
1089     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1090     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1091     nSACell = fCaloCells->GetNumberOfCells();
1092     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1093       // get the cell info *fCal
1094       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1095
1096       // get cell position
1097       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1098       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1099
1100       // convert co global phi eta
1101       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1102       geta = ieta + nColsFeeModule*(nSupMod%2);
1103
1104       // save cell time and energy
1105       cellTime[geta][gphi] = cellTimeT;
1106       cellEnergy[geta][gphi] = cellAmp;
1107
1108     }
1109   }
1110
1111   Int_t nClusterTrig, nClusterTrigCut;
1112   UShort_t *cellAddrs;
1113   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1114   Float_t clsPos[3] = {0.,0.,0.};
1115
1116   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1117   {
1118     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1119     if(!cluster || !cluster->IsEMCAL()) continue;
1120
1121     // get cluster cells
1122     nCell = cluster->GetNCells();
1123
1124     // get cluster energy
1125     clsE = cluster->E();
1126
1127     // get cluster position
1128     cluster->GetPosition(clsPos);
1129     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1130     clsEta = clsPosVec.Eta();
1131     clsPhi = clsPosVec.Phi();
1132
1133     // get the cell addresses
1134     cellAddrs = cluster->GetCellsAbsId();
1135
1136     // check if the cluster contains cell, that was marked as triggered
1137     nClusterTrig = 0;
1138     nClusterTrigCut = 0;
1139
1140     // loop the cells to check, if cluser in acceptance
1141     // any cluster with a cell outside acceptance is not considered
1142     for( iCell = 0; iCell < nCell; iCell++ )
1143     {
1144      // check hot cell
1145      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
1146
1147       // get cell position
1148       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1149       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1150
1151       // convert co global phi eta
1152       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1153       geta = ieta + nColsFeeModule*(nSupMod%2);
1154
1155       if( cellTime[geta][gphi] > 0. ){
1156         clusterTime += cellTime[geta][gphi];
1157         gCell++;
1158       }
1159
1160       // get corresponding FALTRO
1161       fphi = gphi / 2;
1162       feta = geta / 2;
1163
1164       // try to match with a triggered
1165       if( ftriggers[feta][fphi]==1)
1166       {  nClusterTrig++;
1167       }
1168       if( ftriggersCut[feta][fphi]==1)
1169       { nClusterTrigCut++;
1170       }
1171
1172     } // cells
1173
1174
1175     if( gCell > 0 )
1176       clusterTime = clusterTime / (Double_t)gCell;
1177     // fix the reconstriction code time 100ns jumps
1178     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1179       clusterTime -= 0.0000001;
1180
1181     //fClsETime->Fill(clsE,clusterTime);
1182     //fClsEBftTrigCut->Fill(clsE);
1183
1184     if(nClusterTrig>0){
1185       //fClsETime1->Fill(clsE,clusterTime);
1186     }
1187
1188     if(nClusterTrig>0){
1189       cluster->SetChi2(1);
1190       //fClsEAftTrigCut1->Fill(clsE);
1191     }
1192
1193     if(nClusterTrigCut>0){
1194       cluster->SetChi2(2);
1195       //fClsEAftTrigCut2->Fill(clsE);
1196     }
1197
1198     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1199     {
1200       cluster->SetChi2(3);
1201       //fClsEAftTrigCut3->Fill(clsE);
1202     }
1203
1204     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1205     {
1206       // cluster->SetChi2(4);
1207       //fClsEAftTrigCut4->Fill(clsE);
1208     }
1209     if(nClusterTrigCut<1)
1210     {
1211       cluster->SetChi2(0);
1212
1213       //fClsEAftTrigCut->Fill(clsE);
1214     }
1215
1216   } // clusters
1217 }
1218
1219
1220
1221
1222
1223
1224