PWGGA/EMCALJetTasks/AliHadCorrTask.cxx

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