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