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 #include "AliESDpid.h"
  39
  40
  41 using namespace std;
  42
  43 ClassImp(AliAnalysisEtReconstructed);
  44
  45
  46 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
  47         AliAnalysisEt()
  48         ,fCorrections(0)
  49         ,fPidCut(0)
  50         ,fHistChargedPionEnergyDeposit(0)
  51         ,fHistProtonEnergyDeposit(0)
  52         ,fHistAntiProtonEnergyDeposit(0)
  53         ,fHistChargedKaonEnergyDeposit(0)
  54         ,fHistMuonEnergyDeposit(0)
  55         ,fHistRemovedEnergy(0)
  56         ,fGeomCorrection(1.0)
  57         ,fEMinCorrection(1.0/0.687)
  58         ,fRecEffCorrection(1.0)
  59         ,fClusterPositionAccepted(0)
  60         ,fClusterPositionAll(0)
  61         ,fClusterPositionAcceptedEnergy(0)
  62         ,fClusterPositionAllEnergy(0)
  63         ,fClusterEnergy(0)
  64         ,fClusterEnergyCent(0)
  65         ,fClusterEnergyCentMatched(0)
  66         ,fClusterEnergyCentNotMatched(0)
  67         ,fClusterEt(0)
  68         ,fHistChargedEnergyRemoved(0)
  69         ,fHistNeutralEnergyRemoved(0)
  70         ,fHistGammaEnergyAdded(0)
  71         ,fHistMatchedTracksEvspTvsCent(0)
  72         ,fHistMatchedTracksEvspTvsCentEffCorr(0)
  73         ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
  74         ,fHistMatchedTracksEvspTvsCentEffTMCorr500MeV(0)
  75         ,fHistFoundHadronsvsCent(0)
  76         ,fHistNotFoundHadronsvsCent(0)
  77         ,fHistFoundHadronsEtvsCent(0)
  78         ,fHistNotFoundHadronsEtvsCent(0)
  79         ,fHistFoundHadronsvsCent500MeV(0)
  80         ,fHistNotFoundHadronsvsCent500MeV(0)
  81         ,fHistFoundHadronsEtvsCent500MeV(0)
  82         ,fHistNotFoundHadronsEtvsCent500MeV(0)
  83         ,fHistNominalRawEt(0)
  84         ,fHistNominalNonLinHighEt(0)
  85         ,fHistNominalNonLinLowEt(0)
  86         ,fHistNominalEffHighEt(0)
  87         ,fHistNominalEffLowEt(0)
  88         ,fHistTotRawEtEffCorr(0)
  89         ,fHistTotRawEt(0)
  90         ,fHistTotRawEtEffCorr500MeV(0)
  91         ,fHistTotAllRawEt(0)
  92         ,fHistTotAllRawEtEffCorr(0)
  93         ,fHistNClustersPhosVsEmcal(0)
  94         ,fHistClusterSizeVsCent(0)
  95         ,fHistMatchedClusterSizeVsCent(0)
  96         ,fHistTotAllRawEtVsTotalPt(0)
  97         ,fHistTotAllRawEtVsTotalPtVsCent(0)
  98         ,fHistTotMatchedRawEtVsTotalPtVsCent(0)
  99         ,fHistPIDProtonsTrackMatchedDepositedVsNch(0)
 100         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNch(0)
 101         ,fHistPIDProtonsTrackMatchedDepositedVsNcl(0)
 102         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNcl(0)
 103         ,fHistPiKPTrackMatchedDepositedVsNch(0)
 104         ,fHistCentVsNchVsNclReco(0)
 105 {
 106
 107 }
 108
 109 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
 110 {//destructor
 111     delete fCorrections;
 112     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
 113     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
 114     delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
 115     delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
 116     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
 117
 118     delete fHistRemovedEnergy; // removed energy
 119     delete fClusterPositionAccepted;
 120     delete fClusterPositionAll;
 121     delete fClusterPositionAcceptedEnergy;
 122     delete fClusterPositionAllEnergy;
 123     delete fClusterEnergy;
 124     delete fClusterEnergyCent;
 125     delete fClusterEnergyCentMatched;
 126     delete fClusterEnergyCentNotMatched;
 127     delete fClusterEt;
 128     delete fHistChargedEnergyRemoved;
 129     delete fHistNeutralEnergyRemoved;
 130     delete fHistGammaEnergyAdded;
 131     delete fHistMatchedTracksEvspTvsCent;
 132     delete fHistMatchedTracksEvspTvsCentEffCorr;
 133     delete fHistMatchedTracksEvspTvsCentEffTMCorr;
 134     delete fHistMatchedTracksEvspTvsCentEffTMCorr500MeV;
 135     delete fHistFoundHadronsvsCent;
 136     delete fHistNotFoundHadronsvsCent;
 137     delete fHistFoundHadronsEtvsCent;
 138     delete fHistNotFoundHadronsEtvsCent;
 139     delete fHistFoundHadronsvsCent500MeV;
 140     delete fHistNotFoundHadronsvsCent500MeV;
 141     delete fHistFoundHadronsEtvsCent500MeV;
 142     delete fHistNotFoundHadronsEtvsCent500MeV;
 143     delete fHistNominalRawEt;
 144     delete fHistNominalNonLinHighEt;
 145     delete fHistNominalNonLinLowEt;
 146     delete fHistNominalEffHighEt;
 147     delete fHistNominalEffLowEt;
 148     delete fHistTotRawEtEffCorr;
 149     delete fHistTotRawEt;
 150     delete fHistTotAllRawEt;
 151     delete fHistTotAllRawEtEffCorr;
 152     delete fHistTotRawEtEffCorr500MeV;
 153     delete fHistNClustersPhosVsEmcal;
 154     delete fHistClusterSizeVsCent;
 155     delete fHistMatchedClusterSizeVsCent;
 156     delete fHistTotAllRawEtVsTotalPt;
 157     delete fHistTotAllRawEtVsTotalPtVsCent;
 158     delete fHistTotMatchedRawEtVsTotalPtVsCent;
 159     delete fHistPIDProtonsTrackMatchedDepositedVsNch;
 160     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNch;
 161     delete fHistPIDProtonsTrackMatchedDepositedVsNcl;
 162     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNcl;
 163     delete fHistPiKPTrackMatchedDepositedVsNch;
 164     delete fHistCentVsNchVsNclReco;
 165 }
 166
 167 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 168 {
 169
 170     //AliAnalysisEt::AnalyseEvent(ev);
 171     // analyse ESD event
 172     ResetEventValues();
 173     if (!ev) {
 174         AliFatal("ERROR: Event does not exist");
 175         return 0;
 176     }
 177
 178     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
 179     if (!event) {
 180         AliFatal("ERROR: ESD Event does not exist");
 181         return 0;
 182     }
 183     if(!fSelector){
 184         AliFatal("ERROR: fSelector does not exist");
 185         return 0;
 186     }
 187     fSelector->SetEvent(event);
 188
 189     Int_t cent = -1;
 190     fCentrality = event->GetCentrality();
 191     if (fCentrality && cent)
 192     {
 193         cent = fCentrality->GetCentralityClass5("V0M");
 194         fCentClass = fCentrality->GetCentralityClass5("V0M");
 195     }
 196
 197
 198   //for PID
 199   AliESDpid *pID = new AliESDpid();
 200   pID->MakePID(event);
 201   Float_t etPIDProtons = 0.0;
 202   Float_t etPIDAntiProtons = 0.0;
 203   Float_t etPiKPMatched = 0.0;
 204   Float_t multiplicity = fEsdtrackCutsTPC->GetReferenceMultiplicity(event,kTRUE);
 205
 206
 207     Float_t totalMatchedPt = 0.0;
 208     Float_t totalPt = 0.0;
 209     TObjArray* list  = fEsdtrackCutsTPC->GetAcceptedTracks(event);
 210     Int_t nGoodTracks = list->GetEntries();
 211     for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++){
 212         AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
 213         if (!track)
 214           {
 215             Printf("ERROR: Could not get track %d", iTrack);
 216             continue;
 217           }
 218         else{
 219           totalPt +=track->Pt();
 220           pID->MakeITSPID(track);
 221
 222
 223         }
 224     }
 225
 226     //TRefArray *caloClusters = fSelector->GetClusters();//just gets the correct set of clusters - does not apply any cuts
 227     //Float_t fClusterMult = caloClusters->GetEntries();
 228
 229     Float_t nominalRawEt = 0;
 230     Float_t totEt500MeV = 0;
 231     Float_t nonlinHighRawEt = 0;
 232     Float_t nonlinLowRawEt = 0;
 233     Float_t effHighRawEt = 0;
 234     Float_t effLowRawEt = 0;
 235     Float_t uncorrEt = 0;
 236
 237     Float_t nChargedHadronsMeasured = 0.0;
 238     Float_t nChargedHadronsTotal = 0.0;
 239     Float_t nChargedHadronsEtMeasured = 0.0;
 240     Float_t nChargedHadronsEtTotal = 0.0;
 241     Float_t nChargedHadronsMeasured500MeV = 0.0;
 242     Float_t nChargedHadronsTotal500MeV = 0.0;
 243     Float_t nChargedHadronsEtMeasured500MeV = 0.0;
 244     Float_t nChargedHadronsEtTotal500MeV = 0.0;
 245     Float_t fTotAllRawEt = 0.0;
 246     Float_t fTotRawEt = 0.0;
 247     Float_t fTotRawEtEffCorr = 0.0;
 248     Float_t fTotAllRawEtEffCorr = 0.0;
 249     Int_t nPhosClusters = 0;
 250     Int_t nEmcalClusters = 0;
 251
 252     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
 253     {
 254         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 255         if (!cluster)
 256         {
 257             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 258             continue;
 259         }
 260         int x = 0;
 261         fCutFlow->Fill(x++);
 262         if(cluster->IsEMCAL()) nEmcalClusters++;
 263         else nPhosClusters++;
 264         if(!fSelector->IsDetectorCluster(*cluster)) continue;
 265         fCutFlow->Fill(x++);
 266         if(!fSelector->PassMinEnergyCut(*cluster)) continue;
 267         fCutFlow->Fill(x++);
 268         if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
 269         fCutFlow->Fill(x++);
 270         if (!fSelector->CutGeometricalAcceptance(*cluster)) continue;
 271         //fCutFlow->Fill(x++);
 272         Float_t pos[3];
 273
 274         cluster->GetPosition(pos);
 275         TVector3 cp(pos);
 276         fClusterPositionAll->Fill(cp.Phi(), cp.PseudoRapidity());
 277         fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),cluster->E());
 278
 279         //if(TMath::Abs(cp.Eta())> fCuts->fCuts->GetGeometryEmcalEtaAccCut() || cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180. ||  cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.) continue;//Do not accept if cluster is not in the acceptance
 280         fTotAllRawEt += TMath::Sin(cp.Theta())*cluster->E();
 281         fTotAllRawEtEffCorr += CorrectForReconstructionEfficiency(*cluster,cent);
 282
 283         fClusterEnergyCent->Fill(cluster->E(),cent);
 284
 285         Bool_t matched = kTRUE;//default to no track matched
 286         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
 287         matched = !(fSelector->PassTrackMatchingCut(*cluster));//PassTrackMatchingCut is false if there is a matched track
 288         if(matched){//if the track match is good (, is the track good?
 289           if(trackMatchedIndex < 0) matched=kFALSE;//If the index is bad, don't count it
 290           if(matched){
 291             AliESDtrack *track = event->GetTrack(trackMatchedIndex);
 292             //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
 293             matched = fEsdtrackCutsTPC->AcceptTrack(track);//If the track is bad, don't count it
 294           }
 295         }
 296
 297
 298         if (matched)
 299         {
 300
 301             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
 302             {
 303                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
 304                 if (!track) {
 305                     AliError("Error: track does not exist");
 306                 }
 307                 else {
 308                   totalMatchedPt +=track->Pt();
 309                   fClusterEnergyCentMatched->Fill(cluster->E(),cent);
 310                   fHistMatchedClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
 311
 312                   float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cent);
 313                   if(TMath::Abs(eff)<1e-5) eff = 1.0;
 314                   //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
 315                   nChargedHadronsMeasured++;
 316                   nChargedHadronsTotal += 1/eff;
 317                   Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
 318                   nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*cluster->E();
 319                   //One efficiency is the gamma efficiency and the other is the track matching efficiency.
 320                   nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
 321
 322                   Float_t nSigmaPion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
 323                   Float_t nSigmaProton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
 324                   bool isProton = (nSigmaPion>3.0 && nSigmaProton<3.0 && track->Pt()<0.9);
 325                   //cout<<"NSigmaProton "<<nSigmaProton<<endl;
 326                   etPiKPMatched += effCorrEt;
 327                   if(isProton){
 328                     if(track->Charge()>0){
 329                       etPIDProtons += effCorrEt;
 330                     }
 331                     else{
 332                       etPIDAntiProtons += effCorrEt;
 333                     }
 334                   }
 335                   if(TMath::Sin(cp.Theta())*cluster->E()>0.5){
 336                     nChargedHadronsMeasured500MeV++;
 337                     nChargedHadronsTotal500MeV += 1/eff;
 338                     nChargedHadronsEtMeasured500MeV+= TMath::Sin(cp.Theta())*cluster->E();
 339                     nChargedHadronsEtTotal500MeV+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
 340                   }
 341                   fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*cluster->E(),cent);
 342                   fHistMatchedTracksEvspTvsCentEffCorr->Fill(track->P(),effCorrEt,cent);
 343                   //Weighed by the number of tracks we didn't find
 344                   fHistMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
 345                   cluster->GetPosition(pos);
 346                   TVector3 p2(pos);
 347                   uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
 348                   if(uncorrEt>=0.5) fHistMatchedTracksEvspTvsCentEffTMCorr500MeV->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
 349                     const Double_t *pidWeights = track->PID();
 350
 351                     Double_t maxpidweight = 0;
 352                     Int_t maxpid = 0;
 353
 354                     if (pidWeights)
 355                     {
 356                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
 357                         {
 358                             if (pidWeights[p] > maxpidweight)
 359                             {
 360                                 maxpidweight = pidWeights[p];
 361                                 maxpid = p;
 362                             }
 363                         }
 364                         if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
 365                         {
 366                             fEnergyDeposited = cluster->E();
 367                             fMomentumTPC = track->P();
 368                             fCharge = track->Charge();
 369                             fParticlePid = maxpid;
 370                             fPidProb = maxpidweight;
 371                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 372                             if (!esdTrack) {
 373                                 AliError("Error: track does not exist");
 374                             }
 375                             else {
 376                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
 377                                 fDepositTree->Fill();
 378                             }
 379                         }
 380
 381                         if (maxpidweight > fPidCut)
 382                         {
 383                             //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 384
 385                             //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 386
 387                             //Float_t et = cluster->E() * TMath::Sin(theta);
 388                             if (maxpid == AliPID::kProton)
 389                             {
 390
 391                                 if (track->Charge() == 1)
 392                                 {
 393                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
 394                                 }
 395                                 else if (track->Charge() == -1)
 396                                 {
 397                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
 398                                 }
 399                             }
 400                             else if (maxpid == AliPID::kPion)
 401                             {
 402                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
 403                             }
 404                             else if (maxpid == AliPID::kKaon)
 405                             {
 406                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
 407                             }
 408                             else if (maxpid == AliPID::kMuon)
 409                             {
 410                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
 411                             }
 412                         }
 413                     }
 414                 }
 415             }
 416             //continue;
 417         } // distance
 418         else{//these are clusters which were not track matched
 419           fCutFlow->Fill(x++);
 420           //std::cout << x++ << std::endl;
 421
 422           //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 423           //if (cluster->E() < fClusterEnergyCut) continue;
 424           cluster->GetPosition(pos);
 425
 426             TVector3 p2(pos);
 427
 428             fClusterPositionAccepted->Fill(p2.Phi(), p2.PseudoRapidity());
 429             fClusterPositionAcceptedEnergy->Fill(p2.Phi(), p2.PseudoRapidity(),cluster->E());
 430             fClusterEnergy->Fill(cluster->E());
 431             fClusterEnergyCentNotMatched->Fill(cluster->E(),cent);
 432             fHistClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
 433             fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());
 434             uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
 435             float myuncorrEt = TMath::Sin(p2.Theta())*cluster->E();
 436             fTotRawEt += myuncorrEt;
 437
 438             Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
 439             //cout<<"cluster energy "<<cluster->E()<<" eff corr Et "<<effCorrEt<<endl;
 440             fTotRawEtEffCorr += effCorrEt;
 441             fTotNeutralEt += effCorrEt;
 442             nominalRawEt += effCorrEt;
 443             if(myuncorrEt>=0.5){
 444               totEt500MeV += effCorrEt;
 445               //cout<<"test "<<myuncorrEt<<"> 0.5"<<endl;
 446             }
 447             else{
 448               //cout<<"test "<<myuncorrEt<<"< 0.5"<<endl;
 449             }
 450             nonlinHighRawEt += effCorrEt*GetCorrectionModification(*cluster,1,0,cent);
 451             nonlinLowRawEt += effCorrEt*GetCorrectionModification(*cluster,-1,0,cent);
 452             effHighRawEt += effCorrEt*GetCorrectionModification(*cluster,0,1,cent);
 453             effLowRawEt += effCorrEt*GetCorrectionModification(*cluster,0,-1,cent);
 454             fNeutralMultiplicity++;
 455         }
 456         fMultiplicity++;
 457     }
 458
 459     fHistNClustersPhosVsEmcal->Fill(nPhosClusters,nEmcalClusters,cent);
 460     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
 461     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
 462     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
 463     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
 464
 465     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
 466     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 467
 468     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
 469     fHistRemovedEnergy->Fill(removedEnergy);
 470
 471     fTotNeutralEtAcc = fTotNeutralEt;
 472     //fHistTotRawEtEffCorr->Fill(fTotNeutralEt,cent);
 473     fHistTotRawEtEffCorr->Fill(fTotRawEtEffCorr,cent);
 474     fHistTotRawEt->Fill(fTotRawEt,cent);
 475     fHistTotAllRawEt->Fill(fTotAllRawEt,cent);
 476     fHistTotAllRawEtVsTotalPt->Fill(fTotAllRawEt,totalPt);
 477     fHistTotAllRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalPt,cent);
 478     fHistTotMatchedRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalMatchedPt,cent);
 479     fHistTotAllRawEtEffCorr->Fill(fTotAllRawEtEffCorr,cent);
 480     //cout<<"fTotAllRawEtEffCorr "<<fTotAllRawEtEffCorr<<" fTotAllRawEt "<<fTotAllRawEt<<" fTotRawEtEffCorr "<<fTotRawEtEffCorr<<"("<<fTotNeutralEt<<")"<<" fTotRawEt "<<fTotRawEt<<endl;
 481     //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
 482     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
 483     //cout<<" tot corr "<<fTotNeutralEt<<endl;
 484     fTotEt = fTotChargedEt + fTotNeutralEt;
 485 // Fill the histograms...0
 486     FillHistograms();
 487     //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
 488     //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
 489     fHistNominalRawEt->Fill(nominalRawEt,cent);
 490     fHistTotRawEtEffCorr500MeV->Fill(totEt500MeV,cent);
 491     fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
 492     fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
 493     fHistNominalEffHighEt->Fill(effHighRawEt,cent);
 494     fHistNominalEffLowEt->Fill(effLowRawEt,cent);
 495     fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
 496     fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
 497     fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
 498     fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
 499     //cout<<"found "<<nChargedHadronsMeasured<<" total "<<nChargedHadronsTotal<<" not found "<<nChargedHadronsTotal-nChargedHadronsMeasured<<" found "<< nChargedHadronsMeasured500MeV<<" not found "<< nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV <<" total "<<nChargedHadronsTotal500MeV<<endl;
 500     fHistFoundHadronsvsCent500MeV->Fill(nChargedHadronsMeasured500MeV,cent);
 501     fHistNotFoundHadronsvsCent500MeV->Fill(nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV,cent);
 502     fHistFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtMeasured500MeV,cent);
 503     fHistNotFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtTotal500MeV-nChargedHadronsEtMeasured500MeV,cent);
 504 //     cout<<"Number of hadrons measured:  "<<nChargedHadronsMeasured<<" Estimated total number of hadrons "<<nChargedHadronsTotal<<" ET in track matched hadrons "<<
 505 //       nChargedHadronsEtMeasured;
 506 //     if(nChargedHadronsMeasured>0)cout<<" ("<<nChargedHadronsEtMeasured/nChargedHadronsMeasured<<") ";
 507 //     cout<<" ET in all hadrons ";
 508 //     cout<<nChargedHadronsEtTotal;
 509 //     if(nChargedHadronsTotal>0) cout<<" ("<<nChargedHadronsEtTotal/nChargedHadronsTotal<<") ";
 510 //     cout<<endl;
 511     fHistPIDProtonsTrackMatchedDepositedVsNch->Fill(etPIDProtons,multiplicity);
 512     fHistPIDAntiProtonsTrackMatchedDepositedVsNch->Fill(etPIDAntiProtons,multiplicity);
 513     fHistPIDProtonsTrackMatchedDepositedVsNcl->Fill(etPIDProtons,fMultiplicity);
 514     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl->Fill(etPIDAntiProtons,fMultiplicity);
 515     fHistCentVsNchVsNclReco->Fill(cent,multiplicity,fMultiplicity);
 516     fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
 517     delete pID;
 518     return 0;
 519 }
 520
 521 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 522 { // check vertex
 523
 524     Float_t bxy = 999.;
 525     Float_t bz = 999.;
 526     if (!track) {
 527         AliError("ERROR: no track");
 528         return kFALSE;
 529     }
 530     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 531     if (!esdTrack) {
 532         AliError("ERROR: no track");
 533         return kFALSE;
 534     }
 535     esdTrack->GetImpactParametersTPC(bxy,bz);
 536
 537
 538     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 539                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 540                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 541                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 542                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 543
 544     return status;
 545 }
 546
 547 void AliAnalysisEtReconstructed::Init()
 548 { // Init
 549     AliAnalysisEt::Init();
 550     fPidCut = fCuts->GetReconstructedPidCut();
 551     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 552     if (!fCorrections) {
 553         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 554     }
 555 }
 556
 557 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 558 { // propagate track to detector radius
 559     if (!track) {
 560         cout<<"Warning: track empty"<<endl;
 561         return kFALSE;
 562     }
 563     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
 564     if (!esdTrack) {
 565         AliError("ERROR: no ESD track");
 566         return kFALSE;
 567     }
 568     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 569
 570     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 571
 572     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 573     return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
 574 }
 575
 576 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 577 { // add some extra histograms to the ones from base class
 578     AliAnalysisEt::FillOutputList(list);
 579
 580     list->Add(fHistChargedPionEnergyDeposit);
 581     list->Add(fHistProtonEnergyDeposit);
 582     list->Add(fHistAntiProtonEnergyDeposit);
 583     list->Add(fHistChargedKaonEnergyDeposit);
 584     list->Add(fHistMuonEnergyDeposit);
 585
 586     list->Add(fHistRemovedEnergy);
 587     list->Add(fClusterPositionAccepted);
 588     list->Add(fClusterPositionAll);
 589     list->Add(fClusterPositionAcceptedEnergy);
 590     list->Add(fClusterPositionAllEnergy);
 591     list->Add(fClusterEnergy);
 592     list->Add(fClusterEnergyCent);
 593     list->Add(fClusterEnergyCentMatched);
 594     list->Add(fClusterEnergyCentNotMatched);
 595     list->Add(fClusterEt);
 596
 597     list->Add(fHistChargedEnergyRemoved);
 598     list->Add(fHistNeutralEnergyRemoved);
 599     list->Add(fHistGammaEnergyAdded);
 600     list->Add(fHistMatchedTracksEvspTvsCent);
 601     list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
 602     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
 603     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
 604     list->Add(fHistFoundHadronsvsCent);
 605     list->Add(fHistNotFoundHadronsvsCent);
 606     list->Add(fHistFoundHadronsEtvsCent);
 607     list->Add(fHistNotFoundHadronsEtvsCent);
 608     list->Add(fHistFoundHadronsvsCent500MeV);
 609     list->Add(fHistNotFoundHadronsvsCent500MeV);
 610     list->Add(fHistFoundHadronsEtvsCent500MeV);
 611     list->Add(fHistNotFoundHadronsEtvsCent500MeV);
 612     list->Add(fHistNominalRawEt);
 613     list->Add(fHistNominalNonLinHighEt);
 614     list->Add(fHistNominalNonLinLowEt);
 615     list->Add(fHistNominalEffHighEt);
 616     list->Add(fHistNominalEffLowEt);
 617     list->Add(fHistTotRawEtEffCorr);
 618     list->Add(fHistTotRawEtEffCorr500MeV);
 619     list->Add(fHistTotAllRawEtEffCorr);
 620     list->Add(fHistTotRawEt);
 621     list->Add(fHistTotAllRawEt);
 622     list->Add(fHistNClustersPhosVsEmcal);
 623     list->Add(fHistClusterSizeVsCent);
 624     list->Add(fHistMatchedClusterSizeVsCent);
 625     list->Add(fHistTotAllRawEtVsTotalPt);
 626     list->Add(fHistTotAllRawEtVsTotalPtVsCent);
 627     list->Add(fHistTotMatchedRawEtVsTotalPtVsCent);
 628     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNch);
 629     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNch);
 630     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNcl);
 631     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNcl);
 632     list->Add(fHistPiKPTrackMatchedDepositedVsNch);
 633     list->Add(fHistCentVsNchVsNclReco);
 634 }
 635
 636 void AliAnalysisEtReconstructed::CreateHistograms()
 637 { // add some extra histograms to the ones from base class
 638     AliAnalysisEt::CreateHistograms();
 639
 640     Int_t nbinsEt = 1000;
 641     Double_t minEt = 0;
 642     Double_t maxEt = 10;
 643
 644     // possibly change histogram limits
 645 //     if (fCuts) {
 646 //         nbinsEt = fCuts->GetHistNbinsParticleEt();
 647 //         minEt = fCuts->GetHistMinParticleEt();
 648 //         maxEt = fCuts->GetHistMaxParticleEt();
 649 //     }
 650
 651     TString histname;
 652     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
 653     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 654     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 655     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
 656
 657     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
 658     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 659     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 660     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
 661
 662     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
 663     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 664     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 665     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
 666
 667     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
 668     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 669     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 670     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
 671
 672     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
 673     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 674     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 675     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 676
 677     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
 678     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
 679     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 680     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 681
 682     histname = "fClusterPositionAccepted" + fHistogramNameSuffix;
 683     fClusterPositionAccepted = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
 684     fClusterPositionAccepted->SetXTitle("#phi");
 685     fClusterPositionAccepted->SetYTitle("#eta");
 686
 687     histname = "fClusterPositionAll" + fHistogramNameSuffix;
 688     fClusterPositionAll = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
 689     fClusterPositionAll->SetXTitle("#phi");
 690     fClusterPositionAll->SetYTitle("#eta");
 691
 692     histname = "fClusterPositionAcceptedEnergy" + fHistogramNameSuffix;
 693     fClusterPositionAcceptedEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
 694     fClusterPositionAcceptedEnergy->SetXTitle("#phi");
 695     fClusterPositionAcceptedEnergy->SetYTitle("#eta");
 696
 697     histname = "fClusterPositionAllEnergy" + fHistogramNameSuffix;
 698     fClusterPositionAllEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
 699     fClusterPositionAllEnergy->SetXTitle("#phi");
 700     fClusterPositionAllEnergy->SetYTitle("#eta");
 701
 702     histname = "fClusterEnergy" + fHistogramNameSuffix;
 703     fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 704     fClusterEnergy->SetYTitle("Number of clusters");
 705     fClusterEnergy->SetXTitle("Energy of cluster");
 706
 707     histname = "fClusterEnergyCent" + fHistogramNameSuffix;
 708     fClusterEnergyCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
 709     fClusterEnergyCent->SetXTitle("Energy of cluster");
 710     fClusterEnergyCent->SetYTitle("Centrality Bin");
 711     fClusterEnergyCent->SetZTitle("Number of clusters");
 712
 713     histname = "fClusterEnergyCentMatched" + fHistogramNameSuffix;
 714     fClusterEnergyCentMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
 715     fClusterEnergyCentMatched->SetXTitle("Energy of cluster");
 716     fClusterEnergyCentMatched->SetYTitle("Centrality Bin");
 717     fClusterEnergyCentMatched->SetZTitle("Number of Clusters");
 718
 719     histname = "fClusterEnergyCentNotMatched" + fHistogramNameSuffix;
 720     fClusterEnergyCentNotMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
 721     fClusterEnergyCentNotMatched->SetXTitle("Energy of cluster");
 722     fClusterEnergyCentNotMatched->SetYTitle("Centrality Bin");
 723     fClusterEnergyCentNotMatched->SetZTitle("Number of clusters");
 724
 725     histname = "fClusterEt" + fHistogramNameSuffix;
 726     fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 727     fClusterEt->SetXTitle("Number of clusters");
 728     fClusterEt->SetYTitle("E_{T} of cluster");
 729
 730     histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
 731     fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 732
 733     histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
 734     fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 735
 736     histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
 737     fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 738
 739     fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
 740     fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
 741     fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
 742     fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
 743
 744     float max = 200;
 745     if(fHistogramNameSuffix.Contains("P")){max = 100;}
 746     fHistFoundHadronsvsCent = new TH2F("fHistFoundHadronsvsCent","fHistFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
 747        fHistNotFoundHadronsvsCent = new TH2F("fHistNotFoundHadronsvsCent","fHistNotFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
 748     fHistFoundHadronsEtvsCent = new TH2F("fHistFoundHadronsEtvsCent","fHistFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
 749        fHistNotFoundHadronsEtvsCent = new TH2F("fHistNotFoundHadronsEtvsCent","fHistNotFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
 750     fHistFoundHadronsvsCent500MeV = new TH2F("fHistFoundHadronsvsCent500MeV","fHistFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
 751     fHistNotFoundHadronsvsCent500MeV = new TH2F("fHistNotFoundHadronsvsCent500MeV","fHistNotFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
 752     fHistFoundHadronsEtvsCent500MeV = new TH2F("fHistFoundHadronsEtvsCent500MeV","fHistFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
 753     fHistNotFoundHadronsEtvsCent500MeV = new TH2F("fHistNotFoundHadronsEtvsCent500MeV","fHistNotFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
 754
 755     fHistTotRawEtEffCorr = new TH2F("fHistTotRawEtEffCorr","fHistTotRawEtEffCorr",250,0,250,20,-0.5,19.5);
 756     fHistTotRawEt = new TH2F("fHistTotRawEt","fHistTotRawEt",250,0,250,20,-0.5,19.5);
 757     fHistTotRawEtEffCorr500MeV = new TH2F("fHistTotRawEtEffCorr500MeV","fHistTotRawEtEffCorr500MeV",250,0,250,20,-0.5,19.5);
 758     fHistTotAllRawEt = new TH2F("fHistTotAllRawEt","fHistTotAllRawEt",250,0,250,20,-0.5,19.5);
 759     fHistTotAllRawEtEffCorr = new TH2F("fHistTotAllRawEtEffCorr","fHistTotAllRawEtEffCorr",250,0,250,20,-0.5,19.5);
 760     fHistNClustersPhosVsEmcal = new TH3F("fHistNClustersPhosVsEmcal","fHistNClustersPhosVsEmcal",50,0,50,250,0,250,20,-0.5,19);
 761     fHistClusterSizeVsCent = new TH2F("fHistClusterSizeVsCent","fHistClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
 762     fHistMatchedClusterSizeVsCent = new TH2F("fHistMatchedClusterSizeVsCent","fHistMatchedClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
 763     fHistTotAllRawEtVsTotalPt = new TH2F("fHistTotAllRawEtVsTotalPt","fHistTotAllRawEtVsTotalPt",125,0,250,200,0,2000);
 764     fHistTotAllRawEtVsTotalPtVsCent = new TH3F("fHistTotAllRawEtVsTotalPtVsCent","fHistTotAllRawEtVsTotalPtVsCent",125,0,250,200,0,2000,20,-0.5,19.5);
 765     fHistTotMatchedRawEtVsTotalPtVsCent = new TH3F("fHistTotMatchedRawEtVsTotalPtVsCent","fHistTotMatchedRawEtVsTotalPtVsCent",250,0,250,100,0,200,20,-0.5,19.5);
 766
 767     maxEt = 500;
 768     histname = "fHistNominalRawEt" + fHistogramNameSuffix;
 769     fHistNominalRawEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 770     histname = "fHistNominalNonLinHighEt" + fHistogramNameSuffix;
 771     fHistNominalNonLinHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 772     histname = "fHistNominalNonLinLowEt" + fHistogramNameSuffix;
 773     fHistNominalNonLinLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 774     histname = "fHistNominalEffHighEt" + fHistogramNameSuffix;
 775     fHistNominalEffHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 776     histname = "fHistNominalEffLowEt" + fHistogramNameSuffix;
 777     fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
 778
 779       Float_t maxEtRange = 25;
 780       Float_t maxEtRangeHigh = 125;
 781       Float_t minEtRange = 0;
 782       Int_t nbinsMult = 100;
 783       Float_t maxMult = 3000;
 784       Float_t minMult = 0;
 785       Int_t nbinsCl = 150;
 786       Float_t maxCl = 300;
 787       Float_t minCl = 0;
 788     fHistPIDProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNch","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
 789     fHistPIDAntiProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNch","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
 790     fHistPIDProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNcl","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
 791     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
 792     fHistPiKPTrackMatchedDepositedVsNch = new TH2F("fHistPiKPTrackMatchedDepositedVsNch","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
 793     fHistCentVsNchVsNclReco = new TH3F("fHistCentVsNchVsNclReco","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
 794 }
 795 Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)
 796 {
 797   Float_t pos[3];
 798   cluster.GetPosition(pos);
 799   TVector3 cp(pos);
 800   Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),cent);
 801
 802   Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,cent);
 803
 804     cout<<"Warning:  This function should not get called!"<<endl;
 805   //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
 806   return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
 807 }
 808
 809 Double_t AliAnalysisEtReconstructed::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr  0 = nominal 1 = high -1 = low
 810   if(nonLinCorr==0){
 811     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
 812   }
 813   if(effCorr==0){
 814     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
 815   }
 816   return cluster.E()*cent;
 817 }