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