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