PWGLF/totEt/AliAnalysisEtReconstructed.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for ESD analysis
   4 //  - reconstruction output
   5 //  implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
   8 //_________________________________________________________________________
   9
  10 #include "AliAnalysisEtReconstructed.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include "AliEMCALTrack.h"
  14 #include "AliESDCaloCluster.h"
  15 #include "TVector3.h"
  16 #include "TGeoGlobalMagField.h"
  17 #include "AliMagF.h"
  18 #include "AliVEvent.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliVParticle.h"
  22 #include "TDatabasePDG.h"
  23 #include "TList.h"
  24 #include "AliESDpid.h"
  25 #include <iostream>
  26 #include "TH3F.h"
  27 #include "TH2F.h"
  28 #include "TH2I.h"
  29 #include "TH1I.h"
  30 #include "TFile.h"
  31 #include "AliAnalysisHadEtCorrections.h"
  32 #include "AliAnalysisEtSelector.h"
  33 #include "AliLog.h"
  34 #include "AliCentrality.h"
  35 #include "AliPHOSGeoUtils.h"
  36 #include "AliPHOSGeometry.h"
  37 #include "AliAnalysisEtRecEffCorrection.h"
  38
  39
  40 using namespace std;
  41
  42 ClassImp(AliAnalysisEtReconstructed);
  43
  44
  45 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
  46         AliAnalysisEt()
  47         ,fCorrections(0)
  48         ,fPidCut(0)
  49         ,fHistChargedPionEnergyDeposit(0)
  50         ,fHistProtonEnergyDeposit(0)
  51         ,fHistAntiProtonEnergyDeposit(0)
  52         ,fHistChargedKaonEnergyDeposit(0)
  53         ,fHistMuonEnergyDeposit(0)
  54         ,fHistRemovedEnergy(0)
  55         ,fGeomCorrection(1.0)
  56         ,fEMinCorrection(1.0/0.687)
  57         ,fRecEffCorrection(1.0)
  58         ,fClusterPosition(0)
  59         ,fClusterEnergy(0)
  60         ,fClusterEt(0)
  61         ,fHistChargedEnergyRemoved(0)
  62         ,fHistNeutralEnergyRemoved(0)
  63         ,fHistGammaEnergyAdded(0)
  64         ,fHistMatchedTracksEvspTvsCent(0)
  65         ,fHistMatchedTracksEvspTvsCentEffCorr(0)
  66         ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
  67         ,fHistFoundHadronsvsCent(0)
  68         ,fHistNotFoundHadronsvsCent(0)
  69         ,fHistFoundHadronsEtvsCent(0)
  70         ,fHistNotFoundHadronsEtvsCent(0)
  71         ,fHistNominalRawEt(0)
  72         ,fHistNominalNonLinHighEt(0)
  73         ,fHistNominalNonLinLowEt(0)
  74         ,fHistNominalEffHighEt(0)
  75         ,fHistNominalEffLowEt(0)
  76         ,fHistTotRawEt(0)
  77 {
  78
  79 }
  80
  81 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
  82 {//destructor
  83     delete fCorrections;
  84     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
  85     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
  86     delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
  87     delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
  88     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
  89
  90     delete fHistRemovedEnergy; // removed energy
  91     delete fClusterPosition;
  92     delete fClusterEnergy;
  93     delete fClusterEt;
  94     delete fHistChargedEnergyRemoved;
  95     delete fHistNeutralEnergyRemoved;
  96     delete fHistGammaEnergyAdded;
  97     delete fHistMatchedTracksEvspTvsCent;
  98     delete fHistMatchedTracksEvspTvsCentEffCorr;
  99     delete fHistMatchedTracksEvspTvsCentEffTMCorr;
 100     delete fHistFoundHadronsvsCent;
 101     delete fHistNotFoundHadronsvsCent;
 102     delete fHistFoundHadronsEtvsCent;
 103     delete fHistNotFoundHadronsEtvsCent;
 104     delete fHistNominalRawEt;
 105     delete fHistNominalNonLinHighEt;
 106     delete fHistNominalNonLinLowEt;
 107     delete fHistNominalEffHighEt;
 108     delete fHistNominalEffLowEt;
 109     delete fHistTotRawEt;
 110 }
 111
 112 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 113 {
 114
 115     //AliAnalysisEt::AnalyseEvent(ev);
 116     // analyse ESD event
 117     ResetEventValues();
 118     if (!ev) {
 119         AliFatal("ERROR: Event does not exist");
 120         return 0;
 121     }
 122
 123     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
 124     if (!event) {
 125         AliFatal("ERROR: ESD Event does not exist");
 126         return 0;
 127     }
 128     if(!fSelector){
 129         AliFatal("ERROR: fSelector does not exist");
 130         return 0;
 131     }
 132     fSelector->SetEvent(event);
 133
 134     Int_t cent = -1;
 135     fCentrality = event->GetCentrality();
 136     if (fCentrality && cent)
 137     {
 138         cent = fCentrality->GetCentralityClass5("V0M");
 139         fCentClass = fCentrality->GetCentralityClass5("V0M");
 140     }
 141
 142     TRefArray *caloClusters = fSelector->GetClusters();
 143     Float_t fClusterMult = caloClusters->GetEntries();
 144
 145     Float_t nominalRawEt = 0;
 146     Float_t nonlinHighRawEt = 0;
 147     Float_t nonlinLowRawEt = 0;
 148     Float_t effHighRawEt = 0;
 149     Float_t effLowRawEt = 0;
 150     Float_t uncorrEt = 0;
 151
 152     Float_t nChargedHadronsMeasured = 0.0;
 153     Float_t nChargedHadronsTotal = 0.0;
 154     Float_t nChargedHadronsEtMeasured = 0.0;
 155     Float_t nChargedHadronsEtTotal = 0.0;
 156
 157
 158     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
 159     {
 160         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 161         if (!cluster)
 162         {
 163             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 164             continue;
 165         }
 166         int x = 0;
 167         fCutFlow->Fill(x++);
 168         if(!fSelector->IsDetectorCluster(*cluster)) continue;
 169         fCutFlow->Fill(x++);
 170         if(!fSelector->PassMinEnergyCut(*cluster)) continue;
 171         fCutFlow->Fill(x++);
 172         if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
 173         fCutFlow->Fill(x++);
 174
 175         Float_t pos[3];
 176
 177         cluster->GetPosition(pos);
 178         TVector3 cp(pos);
 179
 180         Bool_t matched = kTRUE;//default to no track matched
 181         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
 182         matched = !(fSelector->PassTrackMatchingCut(*cluster));//PassTrackMatchingCut is false if there is a matched track
 183         if(matched){//if the track match is good (, is the track good?
 184           if(trackMatchedIndex < 0) matched=kFALSE;//If the index is bad, don't count it
 185           if(matched){
 186             AliESDtrack *track = event->GetTrack(trackMatchedIndex);
 187             //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
 188             matched = fEsdtrackCutsTPC->AcceptTrack(track);//If the track is bad, don't count it
 189           }
 190         }
 191
 192
 193         if (matched)
 194         {
 195
 196             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
 197             {
 198                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
 199                 if (!track) {
 200                     AliError("Error: track does not exist");
 201                 }
 202                 else {
 203                   float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),fClusterMult);
 204                   if(TMath::Abs(eff)<1e-5) eff = 1.0;
 205                   //cout<<"pt "<<track->Pt()<<" eff "<<eff<<endl;
 206                   nChargedHadronsMeasured++;
 207                   nChargedHadronsTotal += 1/eff;
 208                   Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fClusterMult);
 209                   nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*effCorrEt;
 210                   //One efficiency is the gamma efficiency and the other is the track matching efficiency.
 211                   nChargedHadronsEtTotal+= 1/eff *effCorrEt;
 212                   fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*cluster->E(),cent);
 213                   fHistMatchedTracksEvspTvsCentEffCorr->Fill(track->P(),CorrectForReconstructionEfficiency(*cluster,fClusterMult),cent);
 214                   //Weighed by the number of tracks we didn't find
 215                   fHistMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
 216                     const Double_t *pidWeights = track->PID();
 217
 218                     Double_t maxpidweight = 0;
 219                     Int_t maxpid = 0;
 220
 221                     if (pidWeights)
 222                     {
 223                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
 224                         {
 225                             if (pidWeights[p] > maxpidweight)
 226                             {
 227                                 maxpidweight = pidWeights[p];
 228                                 maxpid = p;
 229                             }
 230                         }
 231                         if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
 232                         {
 233                             fEnergyDeposited = cluster->E();
 234                             fMomentumTPC = track->P();
 235                             fCharge = track->Charge();
 236                             fParticlePid = maxpid;
 237                             fPidProb = maxpidweight;
 238                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 239                             if (!esdTrack) {
 240                                 AliError("Error: track does not exist");
 241                             }
 242                             else {
 243                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
 244                                 fDepositTree->Fill();
 245                             }
 246                         }
 247
 248                         if (maxpidweight > fPidCut)
 249                         {
 250                             //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 251
 252                             //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 253
 254                             //Float_t et = cluster->E() * TMath::Sin(theta);
 255                             if (maxpid == AliPID::kProton)
 256                             {
 257
 258                                 if (track->Charge() == 1)
 259                                 {
 260                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
 261                                 }
 262                                 else if (track->Charge() == -1)
 263                                 {
 264                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
 265                                 }
 266                             }
 267                             else if (maxpid == AliPID::kPion)
 268                             {
 269                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
 270                             }
 271                             else if (maxpid == AliPID::kKaon)
 272                             {
 273                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
 274                             }
 275                             else if (maxpid == AliPID::kMuon)
 276                             {
 277                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
 278                             }
 279                         }
 280                     }
 281                 }
 282             }
 283             //continue;
 284         } // distance
 285         else{//these are clusters which were not track matched
 286           fCutFlow->Fill(x++);
 287           //std::cout << x++ << std::endl;
 288
 289           //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 290           //if (cluster->E() < fClusterEnergyCut) continue;
 291           cluster->GetPosition(pos);
 292
 293             TVector3 p2(pos);
 294
 295             fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
 296             fClusterEnergy->Fill(cluster->E());
 297             fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());
 298             uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
 299
 300             Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fClusterMult);
 301             //cout<<"cluster energy "<<cluster->E()<<" eff corr Et "<<effCorrEt<<endl;
 302             fTotNeutralEt += effCorrEt;
 303             nominalRawEt += effCorrEt;
 304             nonlinHighRawEt += effCorrEt*GetCorrectionModification(*cluster,1,0,fClusterMult);
 305             nonlinLowRawEt += effCorrEt*GetCorrectionModification(*cluster,-1,0,fClusterMult);
 306             effHighRawEt += effCorrEt*GetCorrectionModification(*cluster,0,1,fClusterMult);
 307             effLowRawEt += effCorrEt*GetCorrectionModification(*cluster,0,-1,fClusterMult);
 308             fNeutralMultiplicity++;
 309         }
 310         fMultiplicity++;
 311     }
 312
 313     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
 314     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
 315     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
 316     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
 317
 318     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
 319     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 320
 321     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
 322     fHistRemovedEnergy->Fill(removedEnergy);
 323
 324     fTotNeutralEtAcc = fTotNeutralEt;
 325     fHistTotRawEt->Fill(fTotNeutralEt,cent);
 326     //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
 327     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
 328     //cout<<" tot corr "<<fTotNeutralEt<<endl;
 329     fTotEt = fTotChargedEt + fTotNeutralEt;
 330 // Fill the histograms...0
 331     FillHistograms();
 332     //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
 333     //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
 334     fHistNominalRawEt->Fill(nominalRawEt,cent);
 335     fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
 336     fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
 337     fHistNominalEffHighEt->Fill(effHighRawEt,cent);
 338     fHistNominalEffLowEt->Fill(effLowRawEt,cent);
 339     fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
 340     fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
 341     fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
 342     fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
 343 //     cout<<"Number of hadrons measured:  "<<nChargedHadronsMeasured<<" Estimated total number of hadrons "<<nChargedHadronsTotal<<" ET in track matched hadrons "<<
 344 //       nChargedHadronsEtMeasured;
 345 //     if(nChargedHadronsMeasured>0)cout<<" ("<<nChargedHadronsEtMeasured/nChargedHadronsMeasured<<") ";
 346 //     cout<<" ET in all hadrons ";
 347 //     cout<<nChargedHadronsEtTotal;
 348 //     if(nChargedHadronsTotal>0) cout<<" ("<<nChargedHadronsEtTotal/nChargedHadronsTotal<<") ";
 349 //     cout<<endl;
 350     return 0;
 351 }
 352
 353 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 354 { // check vertex
 355
 356     Float_t bxy = 999.;
 357     Float_t bz = 999.;
 358     if (!track) {
 359         AliError("ERROR: no track");
 360         return kFALSE;
 361     }
 362     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 363     if (!esdTrack) {
 364         AliError("ERROR: no track");
 365         return kFALSE;
 366     }
 367     esdTrack->GetImpactParametersTPC(bxy,bz);
 368
 369
 370     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 371                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 372                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 373                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 374                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 375
 376     return status;
 377 }
 378
 379 void AliAnalysisEtReconstructed::Init()
 380 { // Init
 381     AliAnalysisEt::Init();
 382     fPidCut = fCuts->GetReconstructedPidCut();
 383     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 384     if (!fCorrections) {
 385         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 386     }
 387 }
 388
 389 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 390 { // propagate track to detector radius
 391     if (!track) {
 392         cout<<"Warning: track empty"<<endl;
 393         return kFALSE;
 394     }
 395     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
 396     if (!esdTrack) {
 397         AliError("ERROR: no ESD track");
 398         return kFALSE;
 399     }
 400     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 401
 402     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 403
 404     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 405     return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
 406 }
 407
 408 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 409 { // add some extra histograms to the ones from base class
 410     AliAnalysisEt::FillOutputList(list);
 411
 412     list->Add(fHistChargedPionEnergyDeposit);
 413     list->Add(fHistProtonEnergyDeposit);
 414     list->Add(fHistAntiProtonEnergyDeposit);
 415     list->Add(fHistChargedKaonEnergyDeposit);
 416     list->Add(fHistMuonEnergyDeposit);
 417
 418     list->Add(fHistRemovedEnergy);
 419     list->Add(fClusterPosition);
 420     list->Add(fClusterEnergy);
 421     list->Add(fClusterEt);
 422
 423     list->Add(fHistChargedEnergyRemoved);
 424     list->Add(fHistNeutralEnergyRemoved);
 425     list->Add(fHistGammaEnergyAdded);
 426     list->Add(fHistMatchedTracksEvspTvsCent);
 427     list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
 428     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
 429     list->Add(fHistFoundHadronsvsCent);
 430     list->Add(fHistNotFoundHadronsvsCent);
 431     list->Add(fHistFoundHadronsEtvsCent);
 432     list->Add(fHistNotFoundHadronsEtvsCent);
 433     list->Add(fHistNominalRawEt);
 434     list->Add(fHistNominalNonLinHighEt);
 435     list->Add(fHistNominalNonLinLowEt);
 436     list->Add(fHistNominalEffHighEt);
 437     list->Add(fHistNominalEffLowEt);
 438     list->Add(fHistTotRawEt);
 439 }
 440
 441 void AliAnalysisEtReconstructed::CreateHistograms()
 442 { // add some extra histograms to the ones from base class
 443     AliAnalysisEt::CreateHistograms();
 444
 445     Int_t nbinsEt = 1000;
 446     Double_t minEt = 0;
 447     Double_t maxEt = 10;
 448
 449     // possibly change histogram limits
 450 //     if (fCuts) {
 451 //         nbinsEt = fCuts->GetHistNbinsParticleEt();
 452 //         minEt = fCuts->GetHistMinParticleEt();
 453 //         maxEt = fCuts->GetHistMaxParticleEt();
 454 //     }
 455
 456     TString histname;
 457     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
 458     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 459     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 460     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
 461
 462     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
 463     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 464     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 465     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
 466
 467     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
 468     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 469     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 470     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
 471
 472     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
 473     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 474     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 475     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
 476
 477     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
 478     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 479     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 480     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 481
 482     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
 483     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
 484     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 485     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 486
 487     histname = "fClusterPosition" + fHistogramNameSuffix;
 488     fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, 0,2*TMath::Pi(), 100, -0.7 , 0.7);
 489     fClusterPosition->SetXTitle("#phi");
 490     fClusterPosition->SetYTitle("#eta");
 491
 492     histname = "fClusterEnergy" + fHistogramNameSuffix;
 493     fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 494     fClusterEnergy->SetXTitle("Number of clusters");
 495     fClusterEnergy->SetYTitle("Energy of cluster");
 496
 497     histname = "fClusterEt" + fHistogramNameSuffix;
 498     fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 499     fClusterEt->SetXTitle("Number of clusters");
 500     fClusterEt->SetYTitle("E_{T} of cluster");
 501
 502     histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
 503     fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 504
 505     histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
 506     fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 507
 508     histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
 509     fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 510
 511     fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,0,20);
 512     fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,0,20);
 513     fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,0,20);
 514     fHistFoundHadronsvsCent = new TH2F("fHistFoundHadronsvsCent","fHistFoundHadronsvsCent",100,0,100,20,0,20);
 515     fHistNotFoundHadronsvsCent = new TH2F("fHistNotFoundHadronsvsCent","fHistNotFoundHadronsvsCent",100,0,200,20,0,20);
 516     fHistFoundHadronsEtvsCent = new TH2F("fHistFoundHadronsEtvsCent","fHistFoundHadronsEtvsCent",100,0,200,20,0,20);
 517     fHistNotFoundHadronsEtvsCent = new TH2F("fHistNotFoundHadronsEtvsCent","fHistNotFoundHadronsEtvsCent",100,0,300,20,0,20);
 518
 519     fHistTotRawEt = new TH2F("fHistTotRawEt","fHistTotRawEt",200,0,400,20,0,20);
 520
 521     maxEt = 500;
 522     histname = "fHistNominalRawEt" + fHistogramNameSuffix;
 523     fHistNominalRawEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 524     histname = "fHistNominalNonLinHighEt" + fHistogramNameSuffix;
 525     fHistNominalNonLinHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 526     histname = "fHistNominalNonLinLowEt" + fHistogramNameSuffix;
 527     fHistNominalNonLinLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 528     histname = "fHistNominalEffHighEt" + fHistogramNameSuffix;
 529     fHistNominalEffHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 530     histname = "fHistNominalEffLowEt" + fHistogramNameSuffix;
 531     fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 532
 533 }
 534 Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult)
 535 {
 536   Float_t pos[3];
 537   cluster.GetPosition(pos);
 538   TVector3 cp(pos);
 539   Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),mult);
 540
 541   Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,mult);
 542
 543   //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
 544   return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
 545 }
 546
 547 Double_t AliAnalysisEtReconstructed::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr  0 = nominal 1 = high -1 = low
 548   if(nonLinCorr==0){
 549     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
 550   }
 551   if(effCorr==0){
 552     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
 553   }
 554   return cluster.E()*mult;
 555 }