PWGGA/EMCALJetTasks/AliHadCorrTask.cxx

   1 // $Id$
   2 //
   3 // Hadronic correction task.
   4 //
   5 // Author: R.Reed, C.Loizides, S.Aiola
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12
  13 #include "AliAnalysisManager.h"
  14 #include "AliAODEvent.h"
  15 #include "AliAODCaloCluster.h"
  16 #include "AliESDCaloCluster.h"
  17 #include "AliVTrack.h"
  18 #include "AliPicoTrack.h"
  19 #include "AliVEventHandler.h"
  20 #include "AliEmcalParticle.h"
  21
  22 #include "AliHadCorrTask.h"
  23
  24 ClassImp(AliHadCorrTask)
  25
  26 //________________________________________________________________________
  27 AliHadCorrTask::AliHadCorrTask() :
  28   AliAnalysisTaskEmcal("AliHadCorrTask", kFALSE),
  29   fOutCaloName(),
  30   fPhiMatch(0.05),
  31   fEtaMatch(0.025),
  32   fDoTrackClus(0),
  33   fHadCorr(0),
  34   fEexclCell(0),
  35   fOutClusters(0),
  36   fHistNclusvsCent(0),
  37   fHistNclusMatchvsCent(0),
  38   fHistEbefore(0),
  39   fHistEafter(0),
  40   fHistEoPCent(0),
  41   fHistNMatchCent(0),
  42   fHistCentrality(0)
  43 {
  44   // Default constructor.
  45
  46   for(Int_t i=0; i<8; i++) {
  47       fHistEsubPch[i]    = 0;
  48       fHistEsubPchRat[i] = 0;
  49     for(Int_t j=0; j<4; j++) {
  50       fHistNCellsEnergy[i][j] = 0;
  51     }
  52     if (i<4) {
  53       fHistMatchEvsP[i]    = 0;
  54       fHistMatchdRvsEP[i]  = 0;
  55       fHistNMatchEnergy[i] = 0;
  56     }
  57     for(Int_t j=0; j<9; j++) {
  58       for(Int_t k=0; k<2; k++) {
  59         fHistMatchEtaPhi[i][j][k] = 0;
  60       }
  61     }
  62   }
  63 }
  64
  65 //________________________________________________________________________
  66 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
  67   AliAnalysisTaskEmcal(name, histo),
  68   fOutCaloName("CaloClustersCorr"),
  69   fPhiMatch(0.05),
  70   fEtaMatch(0.025),
  71   fDoTrackClus(1),
  72   fHadCorr(0),
  73   fEexclCell(0),
  74   fOutClusters(0),
  75   fHistNclusvsCent(0),
  76   fHistNclusMatchvsCent(0),
  77   fHistEbefore(0),
  78   fHistEafter(0),
  79   fHistEoPCent(0),
  80   fHistNMatchCent(0),
  81   fHistCentrality(0)
  82 {
  83   // Standard constructor.
  84
  85    for(Int_t i=0; i<8; i++) {
  86     if (i<4) {
  87       fHistMatchEvsP[i]   = 0;
  88       fHistMatchdRvsEP[i] = 0;
  89       for(Int_t j=0; j<3; j++) {
  90         fHistNCellsEnergy[i][j] = 0;
  91       }
  92     }
  93     fHistEsubPch[i]    = 0;
  94     fHistEsubPchRat[i] = 0;
  95     for(Int_t j=0; j<9; j++) {
  96       for(Int_t k=0; k<2; k++) {
  97         fHistMatchEtaPhi[i][j][k] = 0;
  98       }
  99     }
 100   }
 101
 102   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
 103 }
 104
 105 //________________________________________________________________________
 106 AliHadCorrTask::~AliHadCorrTask()
 107 {
 108   // Destructor
 109 }
 110
 111 //________________________________________________________________________
 112 Int_t AliHadCorrTask::GetMomBin(Double_t p) const
 113 {
 114   // Get momenum bin.
 115
 116   Int_t pbin=-1;
 117   if (p<0.5)
 118     pbin=0;
 119   else if (p>=0.5 && p<1.0)
 120     pbin=1;
 121   else if (p>=1.0 && p<1.5)
 122     pbin=2;
 123   else if (p>=1.5 && p<2.)
 124     pbin=3;
 125   else if (p>=2. && p<3.)
 126     pbin=4;
 127   else if (p>=3. && p<4.)
 128     pbin=5;
 129   else if (p>=4. && p<5.)
 130     pbin=6;
 131   else if (p>=5. && p<8.)
 132     pbin=7;
 133   else if (p>=8.)
 134     pbin=8;
 135
 136   return pbin;
 137 }
 138
 139 //________________________________________________________________________
 140 Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
 141 {
 142   // Get sigma in eta.
 143
 144   Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
 145   return 2.0*EtaSigma[pbin];
 146 }
 147
 148 //________________________________________________________________________
 149 Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
 150 {
 151   // Get phi mean.
 152
 153   if (centbin==0) {
 154     Double_t PhiMean[9]={0.036,
 155                          0.021,
 156                          0.0121,
 157                          0.0084,
 158                          0.0060,
 159                          0.0041,
 160                          0.0031,
 161                          0.0022,
 162                          0.001};
 163     return PhiMean[pbin];
 164   } else if (centbin==1) {
 165     Double_t PhiMean[9]={0.036,
 166                          0.021,
 167                          0.0121,
 168                          0.0084,
 169                          0.0060,
 170                          0.0041,
 171                          0.0031,
 172                          0.0022,
 173                          0.001};
 174     return PhiMean[pbin];
 175   } else if (centbin==2) {
 176     Double_t PhiMean[9]={0.036,
 177                          0.021,
 178                          0.0121,
 179                          0.0084,
 180                          0.0060,
 181                          0.0041,
 182                          0.0031,
 183                          0.0022,
 184                          0.001};
 185     return PhiMean[pbin];
 186   } else if (centbin==3) {
 187     Double_t PhiMean[9]={0.036,
 188                          0.021,
 189                          0.0121,
 190                          0.0084,
 191                          0.0060,
 192                          0.0041,
 193                          0.0031,
 194                          0.0022,
 195                          0.001};
 196     return PhiMean[pbin];
 197   } else if (centbin==4) {
 198     Double_t PhiMean[9]={0.036,
 199                          0.021,
 200                          0.0127,
 201                          0.0089,
 202                          0.0068,
 203                          0.0049,
 204                          0.0038,
 205                          0.0028,
 206                          0.0018};
 207     return PhiMean[pbin]*(-1.);
 208   } else if (centbin==5) {
 209     Double_t PhiMean[9]={0.036,
 210                          0.021,
 211                          0.0127,
 212                          0.0089,
 213                          0.0068,
 214                          0.0048,
 215                          0.0038,
 216                          0.0028,
 217                          0.0018};
 218     return PhiMean[pbin]*(-1.);
 219   } else if (centbin==6) {
 220     Double_t PhiMean[9]={0.036,
 221                          0.021,
 222                          0.0127,
 223                          0.0089,
 224                          0.0068,
 225                          0.0045,
 226                          0.0035,
 227                          0.0028,
 228                          0.0018};
 229     return PhiMean[pbin]*(-1.);
 230   } else if (centbin==7) {
 231     Double_t PhiMean[9]={0.036,
 232                          0.021,
 233                          0.0127,
 234                          0.0089,
 235                          0.0068,
 236                          0.0043,
 237                          0.0035,
 238                          0.0028,
 239                          0.0018};
 240     return PhiMean[pbin]*(-1.);
 241   }
 242
 243   return 0;
 244 }
 245
 246 //________________________________________________________________________
 247 Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
 248 {
 249   // Get phi sigma.
 250
 251   if (centbin==0) {
 252     Double_t PhiSigma[9]={0.0221,
 253                           0.0128,
 254                           0.0074,
 255                           0.0064,
 256                           0.0059,
 257                           0.0055,
 258                           0.0052,
 259                           0.0049,
 260                           0.0045};
 261     return 2.*PhiSigma[pbin];
 262   } else if (centbin==1) {
 263     Double_t PhiSigma[9]={0.0217,
 264                           0.0120,
 265                           0.0076,
 266                           0.0066,
 267                           0.0062,
 268                           0.0058,
 269                           0.0054,
 270                           0.0054,
 271 0.0045};
 272     return 2.*PhiSigma[pbin];
 273   } else if (centbin==2) {
 274     Double_t PhiSigma[9]={0.0211,
 275                           0.0124,
 276                           0.0080,
 277                           0.0070,
 278                           0.0067,
 279                           0.0061,
 280                           0.0059,
 281                           0.0054,
 282                           0.0047};
 283     return 2.*PhiSigma[pbin];
 284   } else if (centbin==3) {
 285     Double_t PhiSigma[9]={0.0215,
 286                           0.0124,
 287                           0.0082,
 288                           0.0073,
 289                           0.0069,
 290                           0.0064,
 291                           0.0060,
 292                           0.0055,
 293                           0.0047};
 294     return 2.*PhiSigma[pbin];
 295   } else if (centbin==4) {
 296     Double_t PhiSigma[9]={0.0199,
 297                           0.0108,
 298                           0.0072,
 299                           0.0071,
 300                           0.0060,
 301                           0.0055,
 302                           0.0052,
 303                           0.0049,
 304                           0.0045};
 305     return 2.*PhiSigma[pbin];
 306   } else if (centbin==5) {
 307     Double_t PhiSigma[9]={0.0200,
 308                           0.0110,
 309                           0.0074,
 310                           0.0071,
 311                           0.0064,
 312                           0.0059,
 313                           0.0055,
 314                           0.0052,
 315                           0.0045};
 316     return 2.*PhiSigma[pbin];
 317   } else if (centbin==6) {
 318     Double_t PhiSigma[9]={0.0202,
 319                           0.0113,
 320                           0.0077,
 321                           0.0071,
 322                           0.0069,
 323                           0.0064,
 324                           0.0060,
 325                           0.0055,
 326                           0.0050};
 327     return 2.*PhiSigma[pbin];
 328   } else if (centbin==7) {
 329     Double_t PhiSigma[9]={0.0205,
 330                           0.0113,
 331                           0.0080,
 332                           0.0074,
 333                           0.0078,
 334                           0.0067,
 335                           0.0062,
 336                           0.0055,
 337                           0.0050};
 338     return 2.*PhiSigma[pbin];
 339   }
 340
 341   return 0;
 342 }
 343
 344 //________________________________________________________________________
 345 void AliHadCorrTask::UserCreateOutputObjects()
 346 {
 347   // Create my user objects.
 348
 349   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 350   if (!handler) {
 351     AliError(Form("%s: Input handler not available!", GetName()));
 352     return;
 353   }
 354
 355   if (handler->InheritsFrom("AliESDInputHandler")) {
 356     fOutClusters = new TClonesArray("AliESDCaloCluster");
 357   } else if (handler->InheritsFrom("AliAODInputHandler")) {
 358     fOutClusters = new TClonesArray("AliAODCaloCluster");
 359   } else {
 360     AliError(Form("%s: Input handler not recognized!", GetName()));
 361     return;
 362   }
 363   fOutClusters->SetName(fOutCaloName);
 364
 365   if (!fCreateHisto)
 366     return;
 367
 368   OpenFile(1);
 369   fOutput = new TList();
 370   fOutput->SetOwner();
 371
 372   TString name;
 373
 374   for(Int_t icent=0; icent<8; ++icent) {
 375     for(Int_t ipt=0; ipt<9; ++ipt) {
 376       for(Int_t ieta=0; ieta<2; ++ieta) {
 377         name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
 378         fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, 200, -0.1, 0.1, 400, -0.2, 0.2);
 379         fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
 380       }
 381     }
 382
 383     for(Int_t itrk=0; itrk<4; ++itrk) {
 384       name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
 385       fHistNCellsEnergy[icent][itrk]  = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
 386       fOutput->Add(fHistNCellsEnergy[icent][itrk]);
 387     }
 388
 389     name = Form("fHistEsubPch_%i",icent);
 390     fHistEsubPch[icent]=new TH1F(name, name, 400, 0., 100.);
 391     fOutput->Add(fHistEsubPch[icent]);
 392
 393     name = Form("fHistEsubPchRat_%i",icent);
 394     fHistEsubPchRat[icent]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 395     fOutput->Add(fHistEsubPchRat[icent]);
 396
 397     if (icent<4) {
 398       name = Form("fHistMatchEvsP_%i",icent);
 399       fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 400       fOutput->Add(fHistMatchEvsP[icent]);
 401
 402       name = Form("fHistMatchdRvsEP_%i",icent);
 403       fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
 404       fOutput->Add(fHistMatchdRvsEP[icent]);
 405
 406       name = Form("fHistNMatchEnergy_%i",icent);
 407       fHistNMatchEnergy[icent]  = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
 408       fOutput->Add(fHistNMatchEnergy[icent]);
 409     }
 410   }
 411
 412   fHistCentrality       = new TH1F("fHistCentrality",  "Centrality",       100, 0, 100);
 413   fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
 414   fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
 415   fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
 416   fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
 417   fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
 418   fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 11, -0.5, 10.5);
 419
 420   fOutput->Add(fHistNclusMatchvsCent);
 421   fOutput->Add(fHistNclusvsCent);
 422   fOutput->Add(fHistEbefore);
 423   fOutput->Add(fHistEafter);
 424   fOutput->Add(fHistEoPCent);
 425   fOutput->Add(fHistNMatchCent);
 426   fOutput->Add(fHistCentrality);
 427
 428   PostData(1, fOutput);
 429 }
 430
 431 //________________________________________________________________________
 432 void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
 433 {
 434   // Do the loop over matched tracks for cluster emccluster.
 435
 436   AliVCluster *cluster = emccluster->GetCluster();
 437   Int_t iClus = emccluster->IdInCollection();
 438   Double_t energyclus = cluster->E();
 439
 440   // loop over matched tracks
 441   const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
 442   for (Int_t i = 0; i < Ntrks; ++i) {
 443     Int_t    iTrack = emccluster->GetMatchedObjId(i);
 444     Double_t dR     = emccluster->GetMatchedObjDistance(i);
 445
 446     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
 447     if (!emctrack)
 448       continue;
 449     AliVTrack *track = emctrack->GetTrack();
 450     if (!track)
 451       continue;
 452     if (!AcceptTrack(track))
 453       continue;
 454
 455     // check if track also points to cluster
 456     if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
 457       continue;
 458
 459     Double_t etadiff = 999;
 460     Double_t phidiff = 999;
 461     AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
 462
 463     Double_t mom       = track->P();
 464     Int_t    mombin    = GetMomBin(mom);
 465     Int_t    centbinch = fCentBin;
 466     if (track->Charge()<0)
 467       centbinch += 4;
 468
 469     if (fCreateHisto) {
 470       Int_t etabin = 0;
 471       if(track->Eta() > 0) etabin=1;
 472       fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
 473     }
 474
 475     Double_t etaCut   = 0.0;
 476     Double_t phiCutlo = 0.0;
 477     Double_t phiCuthi = 0.0;
 478
 479     if (fPhiMatch > 0) {
 480       phiCutlo = -fPhiMatch;
 481       phiCuthi = +fPhiMatch;
 482     } else {
 483       phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 484       phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 485     }
 486
 487     if (fEtaMatch > 0) {
 488       etaCut = fEtaMatch;
 489     } else {
 490       etaCut = GetEtaSigma(mombin);
 491     }
 492
 493     if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
 494       ++Nmatches;
 495       totalTrkP += track->P();
 496
 497       if (fCreateHisto) {
 498         if (fHadCorr > 1 && mombin > -1) {
 499           fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
 500         }
 501       }
 502     }
 503   }
 504 }
 505
 506 //________________________________________________________________________
 507 Bool_t AliHadCorrTask::Run()
 508 {
 509   // Run the hadronic correction
 510
 511   // post output in event if not yet present
 512   if (!(InputEvent()->FindListObject(fOutCaloName)))
 513     InputEvent()->AddObject(fOutClusters);
 514
 515   // delete output
 516   fOutClusters->Delete();
 517
 518   // esd or aod mode
 519   Bool_t esdMode = kTRUE;
 520   if (dynamic_cast<AliAODEvent*>(InputEvent()))
 521       esdMode = kFALSE;
 522
 523   if (esdMode) { // optimization in case autobranch loading is off
 524     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 525     if (fCaloName == "CaloClusters")
 526       am->LoadBranch("CaloClusters");
 527     if (fTracksName == "Tracks")
 528       am->LoadBranch("Tracks");
 529     am->LoadBranch("Centrality");
 530   }
 531
 532   if (fCreateHisto) {
 533     fHistCentrality->Fill(fCent);
 534   }
 535
 536    // loop over all clusters
 537   const Int_t Nclus = fCaloClusters->GetEntries();
 538   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
 539     AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 540     if (!emccluster)
 541       continue;
 542
 543     AliVCluster *cluster = emccluster->GetCluster();
 544     if (!cluster)
 545       continue;
 546     if (!AcceptCluster(cluster))
 547       continue;
 548
 549     Double_t energyclus = 0;
 550     if (fCreateHisto)
 551       fHistEbefore->Fill(fCent, cluster->E());
 552
 553     // apply correction / subtraction
 554     if (fHadCorr > 0) {
 555       // to subtract only the closest track set fHadCor to a %
 556       // to subtract all tracks within the cut set fHadCor to %+1
 557       if (fHadCorr > 1)
 558         energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);
 559       else
 560         energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);
 561     }
 562
 563     if (energyclus < 0)
 564       energyclus = 0;
 565
 566     if (energyclus > 0) { // create corrected cluster
 567       AliVCluster *oc;
 568       if (esdMode) {
 569         AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
 570         if (!ec)
 571           continue;
 572         oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
 573       } else {
 574         AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
 575         if (!ac)
 576           continue;
 577         oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
 578       }
 579       oc->SetE(energyclus);
 580
 581       ++clusCount;
 582
 583       if (fCreateHisto)
 584         fHistEafter->Fill(fCent, energyclus);
 585     }
 586   }
 587
 588   return kTRUE;
 589 }
 590
 591 //________________________________________________________________________
 592 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr)
 593 {
 594   // Apply the hadronic correction with one track only.
 595
 596   AliVCluster *cluster = emccluster->GetCluster();
 597   Double_t energyclus  = cluster->E();
 598   Int_t    iMin        = emccluster->GetMatchedObjId();
 599   if (iMin < 0)
 600     return energyclus;
 601
 602   AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 603   if (!emctrack)
 604     return energyclus;
 605
 606   // check if track also points to cluster
 607   Int_t cid = emctrack->GetMatchedObjId();
 608   if (fDoTrackClus && (cid!=emccluster->IdInCollection()))
 609     return energyclus;
 610
 611   AliVTrack *track = emctrack->GetTrack();
 612   if (!track)
 613     return energyclus;
 614
 615   Double_t mom = track->P();
 616   if (mom == 0)
 617     return energyclus;
 618
 619   Double_t dRmin      = emccluster->GetMatchedObjDistance();
 620   Double_t dEtaMin    = 1e9;
 621   Double_t dPhiMin    = 1e9;
 622   AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
 623
 624   Int_t mombin = GetMomBin(mom);
 625   Int_t centbinch = fCentBin;
 626   if (track->Charge()<0)
 627     centbinch += 4;
 628
 629   // plot some histograms if switched on
 630   if (fCreateHisto) {
 631     Int_t etabin = 0;
 632     if(track->Eta() > 0)
 633       etabin = 1;
 634
 635     fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
 636
 637     if (mom > 0) {
 638       fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
 639       fHistEoPCent->Fill(fCent, energyclus / mom);
 640       fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
 641     }
 642   }
 643
 644   // define eta/phi cuts
 645   Double_t etaCut   = 0.0;
 646   Double_t phiCutlo = 0.0;
 647   Double_t phiCuthi = 0.0;
 648   if (fPhiMatch > 0) {
 649     phiCutlo = -fPhiMatch;
 650     phiCuthi = +fPhiMatch;
 651   }
 652   else {
 653     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 654     phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 655   }
 656   if(fEtaMatch > 0) {
 657     etaCut = fEtaMatch;
 658   }
 659   else {
 660     etaCut = GetEtaSigma(mombin);
 661   }
 662
 663   // apply the correction if the track is in the eta/phi window
 664   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
 665     energyclus -= hadCorr * mom;
 666   }
 667
 668   return energyclus;
 669 }
 670
 671 //________________________________________________________________________
 672 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
 673 {
 674   // Apply the hadronic correction with all tracks.
 675
 676   AliVCluster *cluster = emccluster->GetCluster();
 677
 678   Double_t energyclus = cluster->E();
 679   Double_t cNcells = cluster->GetNCells();
 680
 681   Double_t totalTrkP  = 0.0; // count total track momentum
 682   Int_t    Nmatches   = 0;   // count total number of matches
 683
 684   // do the loop over the matched tracks and get the number of matches and the total momentum
 685   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
 686
 687   if(fCreateHisto && Nmatches == 0) {
 688     fHistNclusvsCent->Fill(fCent);
 689     fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
 690   }
 691
 692   if (totalTrkP <= 0)
 693     return energyclus;
 694
 695   Double_t Esub = hadCorr * totalTrkP;
 696
 697   Double_t clusEexcl = fEexclCell * cNcells;
 698
 699   if (energyclus < clusEexcl)
 700     clusEexcl = energyclus;
 701
 702   if (Esub > energyclus)
 703     Esub = energyclus;
 704
 705   // applying Peter's proposed algorithm
 706   // never subtract the full energy of the cluster
 707   if ((energyclus - Esub) < clusEexcl)
 708     Esub = (energyclus - clusEexcl);
 709
 710   Double_t EoP = energyclus / totalTrkP;
 711
 712   // plot some histograms if switched on
 713   if (fCreateHisto) {
 714     fHistNMatchCent->Fill(fCent, Nmatches);
 715     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
 716
 717     if(Nmatches > 0)
 718       fHistNclusMatchvsCent->Fill(fCent);
 719
 720     if(Nmatches == 1)
 721       fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);
 722     else if(Nmatches == 2)
 723       fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
 724     else if(Nmatches > 2)
 725       fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
 726
 727     if (EoP > 0) {
 728       fHistEoPCent->Fill(fCent, EoP);
 729       fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
 730     }
 731
 732     if (Nmatches == 1) {
 733       Int_t iMin = emccluster->GetMatchedObjId();
 734       AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 735       AliVTrack *track = emctrack->GetTrack();
 736       Int_t centbinchm = fCentBin;
 737       if (track->Charge()<0)
 738         centbinchm += 4;
 739
 740       fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
 741       fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
 742     }
 743   }
 744
 745   // apply the correction
 746   energyclus -= Esub;
 747
 748   return energyclus;
 749 }
 750