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 "AliAnalysisTaskHFECal.h"
  39 #include "TGeoGlobalMagField.h"
  40 #include "AliLog.h"
  41 #include "AliAnalysisTaskSE.h"
  42 #include "TRefArray.h"
  43 #include "TVector.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDpid.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliPhysicsSelection.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliAODCaloCluster.h"
  50 #include "AliEMCALRecoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52 #include "AliGeomManager.h"
  53 #include "stdio.h"
  54 #include "TGeoManager.h"
  55 #include "iostream"
  56 #include "fstream"
  57
  58 #include "AliEMCALTrack.h"
  59 #include "AliMagF.h"
  60
  61 #include "AliKFParticle.h"
  62 #include "AliKFVertex.h"
  63
  64 #include "AliPID.h"
  65 #include "AliPIDResponse.h"
  66 #include "AliHFEcontainer.h"
  67 #include "AliHFEcuts.h"
  68 #include "AliHFEpid.h"
  69 #include "AliHFEpidBase.h"
  70 #include "AliHFEpidQAmanager.h"
  71 #include "AliHFEtools.h"
  72 #include "AliCFContainer.h"
  73 #include "AliCFManager.h"
  74
  75 #include "AliCentrality.h"
  76
  77 ClassImp(AliAnalysisTaskHFECal)
  78 //________________________________________________________________________
  79 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
  80   : AliAnalysisTaskSE(name)
  81   ,fESD(0)
  82   ,fGeom(0)
  83   ,fOutputList(0)
  84   ,fTrackCuts(0)
  85   ,fCuts(0)
  86   ,fIdentifiedAsOutInz(kFALSE)
  87   ,fPassTheEventCut(kFALSE)
  88   ,fRejectKinkMother(kFALSE)
  89   ,fVz(0.0)
  90   ,fCFM(0)
  91   ,fPID(0)
  92   ,fPIDqa(0)
  93   ,fOpeningAngleCut(0.1)
  94   ,fInvmassCut(0.01)
  95   ,fNoEvents(0)
  96   ,fEMCAccE(0)
  97   ,fTrkpt(0)
  98   ,fTrkEovPBef(0)
  99   ,fTrkEovPAft(0)
 100   ,fdEdxBef(0)
 101   ,fdEdxAft(0)
 102   ,fIncpT(0)
 103   ,fIncpTM20(0)
 104   ,fInvmassLS(0)
 105   ,fInvmassULS(0)
 106   ,fOpeningAngleLS(0)
 107   ,fOpeningAngleULS(0)
 108   ,fPhotoElecPt(0)
 109   ,fPhoElecPt(0)
 110   ,fPhoElecPtM20(0)
 111   ,fSameElecPt(0)
 112   ,fSameElecPtM20(0)
 113   ,fTrackPtBefTrkCuts(0)
 114   ,fTrackPtAftTrkCuts(0)
 115   ,fTPCnsigma(0)
 116   ,fCent(0)
 117   ,fEleInfo(0)
 118   ,fClsEBftTrigCut(0)
 119   ,fClsEAftTrigCut(0)
 120   ,fClsEAftTrigCut1(0)
 121   ,fClsEAftTrigCut2(0)
 122   ,fClsEAftTrigCut3(0)
 123   ,fClsEAftTrigCut4(0)
 124   ,fClsETime(0)
 125   ,fClsETime1(0)
 126   ,fTrigTimes(0)
 127   ,fCellCheck(0)
 128 {
 129   //Named constructor
 130
 131   fPID = new AliHFEpid("hfePid");
 132   fTrackCuts = new AliESDtrackCuts();
 133
 134   // Define input and output slots here
 135   // Input slot #0 works with a TChain
 136   DefineInput(0, TChain::Class());
 137   // Output slot #0 id reserved by the base class for AOD
 138   // Output slot #1 writes into a TH1 container
 139   // DefineOutput(1, TH1I::Class());
 140   DefineOutput(1, TList::Class());
 141   //  DefineOutput(3, TTree::Class());
 142 }
 143
 144 //________________________________________________________________________
 145 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
 146   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
 147   ,fESD(0)
 148   ,fGeom(0)
 149   ,fOutputList(0)
 150   ,fTrackCuts(0)
 151   ,fCuts(0)
 152   ,fIdentifiedAsOutInz(kFALSE)
 153   ,fPassTheEventCut(kFALSE)
 154   ,fRejectKinkMother(kFALSE)
 155   ,fVz(0.0)
 156   ,fCFM(0)
 157   ,fPID(0)
 158   ,fPIDqa(0)
 159   ,fOpeningAngleCut(0.1)
 160   ,fInvmassCut(0.01)
 161   ,fNoEvents(0)
 162   ,fEMCAccE(0)
 163   ,fTrkpt(0)
 164   ,fTrkEovPBef(0)
 165   ,fTrkEovPAft(0)
 166   ,fdEdxBef(0)
 167   ,fdEdxAft(0)
 168   ,fIncpT(0)
 169   ,fIncpTM20(0)
 170   ,fInvmassLS(0)
 171   ,fInvmassULS(0)
 172   ,fOpeningAngleLS(0)
 173   ,fOpeningAngleULS(0)
 174   ,fPhotoElecPt(0)
 175   ,fPhoElecPt(0)
 176   ,fPhoElecPtM20(0)
 177   ,fSameElecPt(0)
 178   ,fSameElecPtM20(0)
 179   ,fTrackPtBefTrkCuts(0)
 180   ,fTrackPtAftTrkCuts(0)
 181   ,fTPCnsigma(0)
 182   ,fCent(0)
 183   ,fEleInfo(0)
 184   ,fClsEBftTrigCut(0)
 185   ,fClsEAftTrigCut(0)
 186   ,fClsEAftTrigCut1(0)
 187   ,fClsEAftTrigCut2(0)
 188   ,fClsEAftTrigCut3(0)
 189   ,fClsEAftTrigCut4(0)
 190   ,fClsETime(0)
 191   ,fClsETime1(0)
 192   ,fTrigTimes(0)
 193   ,fCellCheck(0)
 194 {
 195         //Default constructor
 196         fPID = new AliHFEpid("hfePid");
 197
 198         fTrackCuts = new AliESDtrackCuts();
 199
 200         // Constructor
 201         // Define input and output slots here
 202         // Input slot #0 works with a TChain
 203         DefineInput(0, TChain::Class());
 204         // Output slot #0 id reserved by the base class for AOD
 205         // Output slot #1 writes into a TH1 container
 206         // DefineOutput(1, TH1I::Class());
 207         DefineOutput(1, TList::Class());
 208         //DefineOutput(3, TTree::Class());
 209 }
 210 //_________________________________________
 211
 212 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
 213 {
 214   //Destructor
 215
 216   delete fOutputList;
 217   delete fGeom;
 218   delete fPID;
 219   delete fCFM;
 220   delete fPIDqa;
 221   delete fTrackCuts;
 222 }
 223 //_________________________________________
 224
 225 void AliAnalysisTaskHFECal::UserExec(Option_t*)
 226 {
 227   //Main loop
 228   //Called for each event
 229
 230   // create pointer to event
 231   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 232   if (!fESD) {
 233     printf("ERROR: fESD not available\n");
 234     return;
 235   }
 236
 237   if(!fCuts){
 238     AliError("HFE cuts not available");
 239     return;
 240   }
 241
 242   if(!fPID->IsInitialized()){
 243     // Initialize PID with the given run number
 244     AliWarning("PID not initialised, get from Run no");
 245     fPID->InitializePID(fESD->GetRunNumber());
 246   }
 247
 248   fNoEvents->Fill(0);
 249
 250   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 251   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 252
 253   Double_t pVtxZ = -999;
 254   pVtxZ = pVtx->GetZ();
 255
 256   if(TMath::Abs(pVtxZ)>10) return;
 257   fNoEvents->Fill(1);
 258
 259   if(fNOtrks<2) return;
 260   fNoEvents->Fill(2);
 261
 262   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 263   if(!pidResponse){
 264     AliDebug(1, "Using default PID Response");
 265     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 266   }
 267
 268   fPID->SetPIDResponse(pidResponse);
 269
 270   fCFM->SetRecEventInfo(fESD);
 271
 272   Float_t cent = -1.;
 273   AliCentrality *centrality = fESD->GetCentrality();
 274   cent = centrality->GetCentralityPercentile("V0M");
 275   fCent->Fill(cent);
 276
 277   if(cent>90.) return;
 278
 279  // Calorimeter info.
 280
 281    FindTriggerClusters();
 282
 283   // make EMCAL array
 284   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
 285      {
 286       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
 287       if(clust->IsEMCAL())
 288         {
 289          double clustE = clust->E();
 290          float  emcx[3]; // cluster pos
 291          clust->GetPosition(emcx);
 292          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 293          double emcphi = clustpos.Phi();
 294          double emceta = clustpos.Eta();
 295          double calInfo[5];
 296          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
 297          if(clustE>1.5)fEMCAccE->Fill(calInfo);
 298         }
 299    }
 300
 301   // Track loop
 302   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 303     AliESDtrack* track = fESD->GetTrack(iTracks);
 304     if (!track) {
 305       printf("ERROR: Could not receive track %d\n", iTracks);
 306       continue;
 307     }
 308
 309     if(TMath::Abs(track->Eta())>0.7) continue;
 310     if(TMath::Abs(track->Pt()<2.0)) continue;
 311
 312     fTrackPtBefTrkCuts->Fill(track->Pt());
 313     // RecKine: ITSTPC cuts
 314     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 315
 316     //RecKink
 317     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 318       if(track->GetKinkIndex(0) != 0) continue;
 319     }
 320
 321     // RecPrim
 322     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 323
 324     // HFEcuts: ITS layers cuts
 325     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 326
 327     // HFE cuts: TPC PID cleanup
 328     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 329
 330
 331     fTrackPtAftTrkCuts->Fill(track->Pt());
 332
 333     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 334     //pt = track->Pt();
 335     //if(pt<2.0)continue;
 336
 337     // Track extrapolation
 338
 339     Int_t charge = track->Charge();
 340     fTrkpt->Fill(pt);
 341     mom = track->P();
 342     phi = track->Phi();
 343     eta = track->Eta();
 344     dEdx = track->GetTPCsignal();
 345     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 346
 347         double ncells = -1.0;
 348         double m20 = -1.0;
 349         double m02 = -1.0;
 350         double disp = -1.0;
 351         double rmatch = -1.0;
 352         double nmatch = -1.0;
 353         double oppstatus = 0.0;
 354         double samestatus = 0.0;
 355
 356     Bool_t fFlagPhotonicElec = kFALSE;
 357     Bool_t fFlagConvinatElec = kFALSE;
 358     SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec);
 359     if(fFlagPhotonicElec)oppstatus = 1.0;
 360     if(fFlagConvinatElec)samestatus = 1.0;
 361
 362     Int_t clsId = track->GetEMCALcluster();
 363     if (clsId>0){
 364       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 365       if(clust && clust->IsEMCAL()){
 366
 367                 double clustE = clust->E();
 368                 eop = clustE/fabs(mom);
 369                 //double clustT = clust->GetTOF();
 370                 ncells = clust->GetNCells();
 371                 m02 = clust->GetM02();
 372                 m20 = clust->GetM20();
 373                 disp = clust->GetDispersion();
 374                 double delphi = clust->GetTrackDx();
 375                 double deleta = clust->GetTrackDz();
 376                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 377                 nmatch = clust->GetNTracksMatched();
 378
 379                   double valdedx[16];
 380                   valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 381                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = disp;
 382                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = samestatus; valdedx[15] = clust->Chi2();
 383                   fEleInfo->Fill(valdedx);
 384
 385
 386       }
 387     }
 388
 389     fdEdxBef->Fill(mom,fTPCnSigma);
 390     fTPCnsigma->Fill(mom,fTPCnSigma);
 391     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 392
 393     Int_t pidpassed = 1;
 394
 395
 396     //--- track accepted
 397     AliHFEpidObject hfetrack;
 398     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 399     hfetrack.SetRecTrack(track);
 400     int centf = (int)cent;
 401     hfetrack.SetCentrality(centf); //added
 402     hfetrack.SetPbPb();
 403     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 404
 405     if(pidpassed==0) continue;
 406
 407     fTrkEovPAft->Fill(pt,eop);
 408     fdEdxAft->Fill(mom,fTPCnSigma);
 409
 410     fIncpT->Fill(cent,pt);
 411     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 412     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 413
 414     if(m20>0.0 && m20<0.3)
 415       {
 416        fIncpTM20->Fill(cent,pt);
 417        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
 418        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
 419      }
 420  }
 421  PostData(1, fOutputList);
 422 }
 423 //_________________________________________
 424 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
 425 {
 426   //--- Check MC
 427
 428   Bool_t mcData = kFALSE;
 429   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 430     {
 431      mcData = kTRUE;
 432      printf("+++++ MC Data available");
 433     }
 434   if(mcData)
 435     {
 436      printf("++++++++= MC analysis \n");
 437     }
 438   else
 439    {
 440      printf("++++++++ real data analysis \n");
 441    }
 442
 443   //---- Geometry
 444   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 445
 446   //--------Initialize PID
 447   //fPID->SetHasMCData(kFALSE);
 448   fPID->SetHasMCData(mcData);
 449   if(!fPID->GetNumberOfPIDdetectors())
 450     {
 451       fPID->AddDetector("TPC", 0);
 452       fPID->AddDetector("EMCAL", 1);
 453     }
 454
 455   fPID->SortDetectors();
 456   fPIDqa = new AliHFEpidQAmanager();
 457   fPIDqa->Initialize(fPID);
 458
 459   //--------Initialize correction Framework and Cuts
 460   fCFM = new AliCFManager;
 461   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 462   fCFM->SetNStepParticle(kNcutSteps);
 463   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 464     fCFM->SetParticleCutsList(istep, NULL);
 465
 466   if(!fCuts){
 467     AliWarning("Cuts not available. Default cuts will be used");
 468     fCuts = new AliHFEcuts;
 469     fCuts->CreateStandardCuts();
 470   }
 471   fCuts->Initialize(fCFM);
 472
 473   //---------Output Tlist
 474   fOutputList = new TList();
 475   fOutputList->SetOwner();
 476   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 477
 478   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
 479   fOutputList->Add(fNoEvents);
 480
 481   Int_t binsE[5] =    {250, 100,  1000, 100,   10};
 482   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
 483   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
 484   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
 485   fOutputList->Add(fEMCAccE);
 486
 487   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 488   fOutputList->Add(fTrkpt);
 489
 490   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 491   fOutputList->Add(fTrackPtBefTrkCuts);
 492
 493   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 494   fOutputList->Add(fTrackPtAftTrkCuts);
 495
 496   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 497   fOutputList->Add(fTPCnsigma);
 498
 499   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 500   fOutputList->Add(fTrkEovPBef);
 501
 502   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 503   fOutputList->Add(fTrkEovPAft);
 504
 505   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
 506   fOutputList->Add(fdEdxBef);
 507
 508   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
 509   fOutputList->Add(fdEdxAft);
 510
 511   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
 512   fOutputList->Add(fIncpT);
 513
 514   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
 515   fOutputList->Add(fIncpTM20);
 516
 517   Int_t nBinspho[3] =  { 100, 100, 500};
 518   Double_t minpho[3] = {  0.,  0., 0.};
 519   Double_t maxpho[3] = {100., 50., 0.5};
 520
 521   fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
 522   fOutputList->Add(fInvmassLS);
 523
 524   fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
 525   fOutputList->Add(fInvmassULS);
 526
 527   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 528   fOutputList->Add(fOpeningAngleLS);
 529
 530   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 531   fOutputList->Add(fOpeningAngleULS);
 532
 533   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 534   fOutputList->Add(fPhotoElecPt);
 535
 536   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
 537   fOutputList->Add(fPhoElecPt);
 538
 539   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",100,0,100,100,0,50);
 540   fOutputList->Add(fPhoElecPtM20);
 541
 542   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
 543   fOutputList->Add(fSameElecPt);
 544
 545   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",100,0,100,100,0,50);
 546   fOutputList->Add(fSameElecPtM20);
 547
 548   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 549   fOutputList->Add(fCent);
 550
 551   // Make common binning
 552   const Double_t kMinP = 2.;
 553   const Double_t kMaxP = 50.;
 554   const Double_t kTPCSigMim = 40.;
 555   const Double_t kTPCSigMax = 140.;
 556
 557   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
 558   Int_t nBins[16] =  {  480,        200,   60,    20,   600,  300, 100,   40,   200, 200, 200, 100,    3,    3,    3,   10};
 559   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, -0.5};
 560   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,  2.5,  2.5,  9.5};
 561   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;", 16, nBins, min, max);
 562   fOutputList->Add(fEleInfo);
 563
 564   //<---  Trigger info
 565
 566   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
 567   fOutputList->Add(fClsEBftTrigCut);
 568
 569   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
 570   fOutputList->Add(fClsEAftTrigCut);
 571
 572   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
 573   fOutputList->Add(fClsEAftTrigCut1);
 574
 575   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
 576   fOutputList->Add(fClsEAftTrigCut2);
 577
 578   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
 579   fOutputList->Add(fClsEAftTrigCut3);
 580
 581   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
 582   fOutputList->Add(fClsEAftTrigCut4);
 583
 584   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 585   fOutputList->Add(fClsETime);
 586
 587   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
 588   fOutputList->Add(fClsETime1);
 589
 590   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
 591   fOutputList->Add(fTrigTimes);
 592
 593   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
 594   fOutputList->Add(fCellCheck);
 595
 596   PostData(1,fOutputList);
 597 }
 598
 599 //________________________________________________________________________
 600 void AliAnalysisTaskHFECal::Terminate(Option_t *)
 601 {
 602   // Info("Terminate");
 603         AliAnalysisTaskSE::Terminate();
 604 }
 605
 606 //________________________________________________________________________
 607 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 608 {
 609   // Check single track cuts for a given cut step
 610   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 611   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 612   return kTRUE;
 613 }
 614 //_________________________________________
 615 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
 616 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec)
 617 {
 618   //Identify non-heavy flavour electrons using Invariant mass method
 619
 620   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 621   //fTrackCuts->SetRequireTPCRefit(kTRUE);
 622   //fTrackCuts->SetEtaRange(-0.7,0.7);
 623   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 624   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 625   fTrackCuts->SetMinNClustersTPC(70);
 626
 627   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 628
 629   Bool_t flagPhotonicElec = kFALSE;
 630   Bool_t flagConvinatElec = kFALSE;
 631
 632   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 633   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 634     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 635     if (!trackAsso) {
 636       printf("ERROR: Could not receive track %d\n", jTracks);
 637       continue;
 638     }
 639     if(itrack==jTracks)continue;
 640
 641     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
 642     Double_t mass=999., width = -999;
 643     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 644
 645     ptPrim = track->Pt();
 646
 647     dEdxAsso = trackAsso->GetTPCsignal();
 648     ptAsso = trackAsso->Pt();
 649     Int_t chargeAsso = trackAsso->Charge();
 650     Int_t charge = track->Charge();
 651
 652     //if(ptAsso <0.3) continue;
 653     if(ptAsso <0.5) continue;
 654     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 655     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
 656
 657     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 658     if(charge>0) fPDGe1 = -11;
 659     if(chargeAsso>0) fPDGe2 = -11;
 660
 661     if(charge == chargeAsso) fFlagLS = kTRUE;
 662     if(charge != chargeAsso) fFlagULS = kTRUE;
 663
 664     AliKFParticle ge1(*track, fPDGe1);
 665     AliKFParticle ge2(*trackAsso, fPDGe2);
 666     AliKFParticle recg(ge1, ge2);
 667
 668     if(recg.GetNDF()<1) continue;
 669     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 670     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 671
 672     AliKFVertex primV(*pVtx);
 673     primV += recg;
 674     recg.SetProductionVertex(primV);
 675
 676     recg.SetMassConstraint(0,0.0001);
 677
 678     openingAngle = ge1.GetAngle(ge2);
 679     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 680     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 681
 682     if(openingAngle > fOpeningAngleCut) continue;
 683
 684     recg.GetMass(mass,width);
 685
 686     double phoinfo[3];
 687     phoinfo[0] = cent;
 688     phoinfo[1] = ptPrim;
 689     phoinfo[2] = mass;
 690
 691     if(fFlagLS) fInvmassLS->Fill(phoinfo);
 692     if(fFlagULS) fInvmassULS->Fill(phoinfo);
 693
 694     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 695       flagPhotonicElec = kTRUE;
 696     }
 697     if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
 698       flagConvinatElec = kTRUE;
 699     }
 700
 701   }
 702   fFlagPhotonicElec = flagPhotonicElec;
 703   fFlagConvinatElec = flagConvinatElec;
 704
 705 }
 706
 707
 708 //_________________________________________
 709 void AliAnalysisTaskHFECal::FindTriggerClusters()
 710 {
 711   // constants
 712   const int nModuleCols = 2;
 713   const int nModuleRows = 5;
 714   const int nColsFeeModule = 48;
 715   const int nRowsFeeModule = 24;
 716   const int nColsFaltroModule = 24;
 717   const int nRowsFaltroModule = 12;
 718   //const int faltroWidthMax = 20;
 719
 720   // part 1, trigger extraction -------------------------------------
 721   Int_t globCol, globRow;
 722   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
 723   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
 724
 725   //Int_t trigtimes[faltroWidthMax];
 726   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
 727   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
 728   //Double_t fTrigCutLow = 6;
 729   //Double_t fTrigCutHigh = 10;
 730   Double_t fTimeCutLow =  469e-09;
 731   Double_t fTimeCutHigh = 715e-09;
 732
 733   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
 734
 735   // erase trigger maps
 736   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
 737   {
 738     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
 739     {
 740       ftriggersCut[i][j] = 0;
 741       ftriggers[i][j] = 0;
 742       ftriggersTime[i][j] = 0;
 743     }
 744   }
 745
 746   Int_t iglobCol=0, iglobRow=0;
 747   // go through triggers
 748   if( fCaloTrigger->GetEntries() > 0 )
 749   {
 750     // needs reset
 751     fCaloTrigger->Reset();
 752     while( fCaloTrigger->Next() )
 753     {
 754       fCaloTrigger->GetPosition( globCol, globRow );
 755       fCaloTrigger->GetNL0Times( ntimes );
 756       /*
 757       // no L0s
 758       if( ntimes < 1 )   continue;
 759       // get precise timings
 760       fCaloTrigger->GetL0Times( trigtimes );
 761       trigInCut = 0;
 762       for(Int_t i = 0; i < ntimes; i++ )
 763       {
 764         // save the first trigger time in channel
 765         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
 766           triggersTime[globCol][globRow] = trigtimes[i];
 767         //printf("trigger times: %d\n",trigtimes[i]);
 768         // check if in cut
 769         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
 770           trigInCut = 1;
 771
 772         fTrigTimes->Fill(trigtimes[i]);
 773       }
 774       */
 775
 776       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
 777       Int_t bit = 0;
 778       fCaloTrigger->GetTriggerBits(bit);
 779
 780       Int_t ts = 0;
 781       fCaloTrigger->GetL1TimeSum(ts);
 782       if (ts > 0)ftriggers[globCol][globRow] = 1;
 783       // number of triggered channels in event
 784       nTrigChannel++;
 785       // ... inside cut
 786       if(ts>0 && (bit >> 6 & 0x1))
 787       {
 788         iglobCol = globCol;
 789         iglobRow = globRow;
 790         nTrigChannelCut++;
 791         ftriggersCut[globCol][globRow] = 1;
 792       }
 793
 794     } // calo trigger entries
 795   } // has calo trigger entries
 796
 797   // part 2 go through the clusters here -----------------------------------
 798   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
 799   Short_t cellAddr, nSACell, mclabel;
 800   //Int_t nSACell, iSACell, mclabel;
 801   Int_t iSACell;
 802   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
 803   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
 804
 805   TRefArray *fCaloClusters = new TRefArray();
 806   fESD->GetEMCALClusters( fCaloClusters );
 807   nCluster = fCaloClusters->GetEntries();
 808
 809
 810   // save all cells times if there are clusters
 811   if( nCluster > 0 ){
 812     // erase time map
 813     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
 814       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
 815         cellTime[i][j] = 0.;
 816         cellEnergy[i][j] = 0.;
 817       }
 818     }
 819
 820     // get cells
 821     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
 822     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
 823     nSACell = fCaloCells->GetNumberOfCells();
 824     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
 825       // get the cell info *fCal
 826       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
 827
 828       // get cell position
 829       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
 830       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
 831
 832       // convert co global phi eta
 833       gphi = iphi + nRowsFeeModule*(nSupMod/2);
 834       geta = ieta + nColsFeeModule*(nSupMod%2);
 835
 836       // save cell time and energy
 837       cellTime[geta][gphi] = cellTimeT;
 838       cellEnergy[geta][gphi] = cellAmp;
 839
 840     }
 841   }
 842
 843   Int_t nClusterTrig, nClusterTrigCut;
 844   UShort_t *cellAddrs;
 845   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
 846   Float_t clsPos[3] = {0.,0.,0.};
 847
 848   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
 849   {
 850     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
 851     if(!cluster || !cluster->IsEMCAL()) continue;
 852
 853     // get cluster cells
 854     nCell = cluster->GetNCells();
 855
 856     // get cluster energy
 857     clsE = cluster->E();
 858
 859     // get cluster position
 860     cluster->GetPosition(clsPos);
 861     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
 862     clsEta = clsPosVec.Eta();
 863     clsPhi = clsPosVec.Phi();
 864
 865     // get the cell addresses
 866     cellAddrs = cluster->GetCellsAbsId();
 867
 868     // check if the cluster contains cell, that was marked as triggered
 869     nClusterTrig = 0;
 870     nClusterTrigCut = 0;
 871
 872     // loop the cells to check, if cluser in acceptance
 873     // any cluster with a cell outside acceptance is not considered
 874     for( iCell = 0; iCell < nCell; iCell++ )
 875     {
 876      // check hot cell
 877      if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
 878
 879       // get cell position
 880       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
 881       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
 882
 883       // convert co global phi eta
 884       gphi = iphi + nRowsFeeModule*(nSupMod/2);
 885       geta = ieta + nColsFeeModule*(nSupMod%2);
 886
 887       if( cellTime[geta][gphi] > 0. ){
 888         clusterTime += cellTime[geta][gphi];
 889         gCell++;
 890       }
 891
 892       // get corresponding FALTRO
 893       fphi = gphi / 2;
 894       feta = geta / 2;
 895
 896       // try to match with a triggered
 897       if( ftriggers[feta][fphi]==1)
 898       {  nClusterTrig++;
 899       }
 900       if( ftriggersCut[feta][fphi]==1)
 901       { nClusterTrigCut++;
 902       }
 903
 904     } // cells
 905
 906
 907     if( gCell > 0 )
 908       clusterTime = clusterTime / (Double_t)gCell;
 909     // fix the reconstriction code time 100ns jumps
 910     if( fESD->GetBunchCrossNumber() % 4 < 2 )
 911       clusterTime -= 0.0000001;
 912
 913     fClsETime->Fill(clsE,clusterTime);
 914     fClsEBftTrigCut->Fill(clsE);
 915
 916     if(nClusterTrig>0){
 917       fClsETime1->Fill(clsE,clusterTime);
 918     }
 919
 920     if(nClusterTrig>0){
 921       cluster->SetChi2(1);
 922       fClsEAftTrigCut1->Fill(clsE);
 923     }
 924
 925     if(nClusterTrigCut>0){
 926       cluster->SetChi2(2);
 927       fClsEAftTrigCut2->Fill(clsE);
 928     }
 929
 930     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
 931     {
 932       cluster->SetChi2(3);
 933       fClsEAftTrigCut3->Fill(clsE);
 934     }
 935
 936     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
 937     {
 938       // cluster->SetChi2(4);
 939       fClsEAftTrigCut4->Fill(clsE);
 940     }
 941     if(nClusterTrigCut<1)
 942     {
 943       cluster->SetChi2(0);
 944
 945       fClsEAftTrigCut->Fill(clsE);
 946     }
 947
 948   } // clusters
 949 }
 950
 951
 952
 953
 954
 955
 956