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