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