PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for MC analysis
   4 //  - MC output
   5 //  implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
   8 //_________________________________________________________________________
   9
  10 #include "AliAnalysisEtMonteCarlo.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliAnalysisEtSelector.h"
  13 #include "AliAnalysisEtSelector.h"
  14 #include "AliESDtrack.h"
  15 #include "AliStack.h"
  16 #include "AliVEvent.h"
  17 #include "AliMCEvent.h"
  18 #include "AliESDEvent.h"
  19 #include "TH2F.h"
  20 #include "TH3F.h"
  21 #include "TParticle.h"
  22 #include "AliGenHijingEventHeader.h"
  23 #include "AliGenPythiaEventHeader.h"
  24 #include "TList.h"
  25 #include "AliESDCaloCluster.h"
  26 #include "AliLog.h"
  27 #include <iostream>
  28 #include <AliCentrality.h>
  29 #include "AliPHOSGeoUtils.h"
  30 #include "AliPHOSGeometry.h"
  31 #include "TFile.h"
  32 #include "AliESDtrackCuts.h"
  33 #include "AliGenEventHeader.h"
  34 #include "AliGenCocktailEventHeader.h"
  35 #include "AliGenEventHeader.h"
  36 using namespace std;
  37
  38 ClassImp(AliAnalysisEtMonteCarlo);
  39
  40
  41 // ctor
  42 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
  43                                                   ,nChargedHadronsMeasured(0)
  44                                                   ,nChargedHadronsTotal(0)
  45                                                   ,fIsMC(kTRUE)
  46                                                   ,checkLabelForHIJING(kFALSE)
  47     ,fImpactParameter(0)
  48     ,fNcoll(0)
  49     ,fNpart(0)
  50     ,fPrimaryTree(0)
  51     ,fTotEtWithSecondaryRemoved(0)
  52     ,fTotEtSecondaryFromEmEtPrimary(0)
  53     ,fTotEtSecondary(0)
  54     ,fPrimaryCode(0)
  55     ,fPrimaryCharge(0)
  56     ,fPrimaryE(0)
  57     ,fPrimaryEt(0)
  58     ,fPrimaryPx(0)
  59     ,fPrimaryPy(0)
  60     ,fPrimaryPz(0)
  61     ,fPrimaryVx(0)
  62     ,fPrimaryVy(0)
  63     ,fPrimaryVz(0)
  64     ,fPrimaryAccepted(0)
  65     ,fPrimaryMatched(0)
  66     ,fDepositedCode(0)
  67     ,fDepositedE(0)
  68     ,fDepositedEt(0)
  69     ,fDepositedCharge(0)
  70     ,fDepositedVx(0)
  71     ,fDepositedVy(0)
  72     ,fDepositedVz(0)
  73     ,fSecondary(kFALSE)
  74     ,fReconstructedE(0)
  75     ,fReconstructedEt(0)
  76     ,fTotPx(0)
  77     ,fTotPy(0)
  78     ,fTotPz(0)
  79     ,fClusterMult(0)
  80     ,fHistDecayVertexNonRemovedCharged(0)
  81     ,fHistDecayVertexRemovedCharged(0)
  82     ,fHistDecayVertexNonRemovedNeutral(0)
  83     ,fHistDecayVertexRemovedNeutral(0)
  84
  85     ,fHistRemovedOrNot(0)
  86     ,fHistEtNonRemovedProtons(0)
  87     ,fHistEtNonRemovedAntiProtons(0)
  88     ,fHistEtNonRemovedPiPlus(0)
  89     ,fHistEtNonRemovedPiMinus(0)
  90     ,fHistEtNonRemovedKaonPlus(0)
  91     ,fHistEtNonRemovedKaonMinus(0)
  92     ,fHistEtNonRemovedK0s(0)
  93     ,fHistEtNonRemovedK0L(0)
  94     ,fHistEtNonRemovedLambdas(0)
  95     ,fHistEtNonRemovedElectrons(0)
  96     ,fHistEtNonRemovedPositrons(0)
  97     ,fHistEtNonRemovedMuPlus(0)
  98     ,fHistEtNonRemovedMuMinus(0)
  99     ,fHistEtNonRemovedNeutrons(0)
 100     ,fHistEtNonRemovedAntiNeutrons(0)
 101     ,fHistEtNonRemovedGammas(0)
 102     ,fHistEtNonRemovedGammasFromPi0(0)
 103     ,fHistEtRemovedGammas(0)
 104     ,fHistEtRemovedNeutrons(0)
 105     ,fHistEtRemovedAntiNeutrons(0)
 106     ,fHistEtRemovedCharged(0)
 107     ,fHistEtRemovedNeutrals(0)
 108     ,fHistEtNonRemovedCharged(0)
 109     ,fHistEtNonRemovedNeutrals(0)
 110     ,fHistMultNonRemovedProtons(0)
 111     ,fHistMultNonRemovedAntiProtons(0)
 112     ,fHistMultNonRemovedPiPlus(0)
 113     ,fHistMultNonRemovedPiMinus(0)
 114     ,fHistMultNonRemovedKaonPlus(0)
 115     ,fHistMultNonRemovedKaonMinus(0)
 116     ,fHistMultNonRemovedK0s(0)
 117     ,fHistMultNonRemovedK0L(0)
 118     ,fHistMultNonRemovedLambdas(0)
 119     ,fHistMultNonRemovedElectrons(0)
 120     ,fHistMultNonRemovedPositrons(0)
 121     ,fHistMultNonRemovedMuPlus(0)
 122     ,fHistMultNonRemovedMuMinus(0)
 123     ,fHistMultNonRemovedNeutrons(0)
 124     ,fHistMultNonRemovedAntiNeutrons(0)
 125     ,fHistMultNonRemovedGammas(0)
 126     ,fHistMultRemovedGammas(0)
 127     ,fHistMultRemovedNeutrons(0)
 128     ,fHistMultRemovedAntiNeutrons(0)
 129     ,fHistMultRemovedCharged(0)
 130     ,fHistMultRemovedNeutrals(0)
 131     ,fHistMultNonRemovedCharged(0)
 132     ,fHistMultNonRemovedNeutrals(0)
 133     ,fHistTrackMultvsNonRemovedCharged(0)
 134     ,fHistTrackMultvsNonRemovedNeutral(0)
 135     ,fHistTrackMultvsRemovedGamma(0)
 136     ,fHistClusterMultvsNonRemovedCharged(0)
 137     ,fHistClusterMultvsNonRemovedNeutral(0)
 138     ,fHistClusterMultvsRemovedGamma(0)
 139     ,fHistMultvsNonRemovedChargedE(0)
 140     ,fHistMultvsNonRemovedNeutralE(0)
 141     ,fHistMultvsRemovedGammaE(0)
 142                                                   ,fCalcForKaonCorrection(kFALSE)
 143                                                   ,fHistK0EDepositsVsPtInAcceptance(0)
 144                                                   ,fHistK0EGammaVsPtInAcceptance(0)
 145                                                   ,fHistK0EDepositsVsPtOutOfAcceptance(0)
 146                                                   ,fHistK0EGammaVsPtOutOfAcceptance(0)
 147                                                   ,fHistSimKaonsInAcceptance(0)
 148                                                   ,fHistSimK0SInAcceptance(0)
 149                                                   ,fHistSimKPlusInAcceptance(0)
 150                                                   ,fHistSimKMinusInAcceptance(0)
 151                                                   ,fHistSimK0LInAcceptance(0)
 152                                                   ,fHistSimKaonsOutOfAcceptance(0)
 153                                                   ,fHistSimKaonsInAcceptanceWithDepositsPrimaries(0)
 154                                                   ,fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries(0)
 155                                                   ,fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries(0)
 156     ,fEtNonRemovedProtons(0)
 157     ,fEtNonRemovedAntiProtons(0)
 158     ,fEtNonRemovedPiPlus(0)
 159     ,fEtNonRemovedPiMinus(0)
 160     ,fEtNonRemovedKaonPlus(0)
 161     ,fEtNonRemovedKaonMinus(0)
 162     ,fEtNonRemovedK0S(0)
 163     ,fEtNonRemovedK0L(0)
 164     ,fEtNonRemovedLambdas(0)
 165     ,fEtNonRemovedElectrons(0)
 166     ,fEtNonRemovedPositrons(0)
 167     ,fEtNonRemovedMuMinus(0)
 168     ,fEtNonRemovedMuPlus(0)
 169     ,fEtNonRemovedGammas(0)
 170     ,fEtNonRemovedGammasFromPi0(0)
 171     ,fEtNonRemovedNeutrons(0)
 172     ,fEtNonRemovedAntiNeutrons(0)
 173     ,fEtRemovedProtons(0)
 174     ,fEtRemovedAntiProtons(0)
 175     ,fEtRemovedPiPlus(0)
 176     ,fEtRemovedPiMinus(0)
 177     ,fEtRemovedKaonPlus(0)
 178     ,fEtRemovedKaonMinus(0)
 179     ,fEtRemovedK0s(0)
 180     ,fEtRemovedK0L(0)
 181     ,fEtRemovedLambdas(0)
 182     ,fEtRemovedElectrons(0)
 183     ,fEtRemovedPositrons(0)
 184     ,fEtRemovedMuMinus(0)
 185     ,fEtRemovedMuPlus(0)
 186     ,fEtRemovedGammasFromPi0(0)
 187     ,fEtRemovedGammas(0)
 188     ,fEtRemovedNeutrons(0)
 189     ,fEtRemovedAntiNeutrons(0)
 190     ,fMultNonRemovedProtons(0)
 191     ,fMultNonRemovedAntiProtons(0)
 192     ,fMultNonRemovedPiPlus(0)
 193     ,fMultNonRemovedPiMinus(0)
 194     ,fMultNonRemovedKaonPlus(0)
 195     ,fMultNonRemovedKaonMinus(0)
 196     ,fMultNonRemovedK0s(0)
 197     ,fMultNonRemovedK0L(0)
 198     ,fMultNonRemovedLambdas(0)
 199     ,fMultNonRemovedElectrons(0)
 200     ,fMultNonRemovedPositrons(0)
 201     ,fMultNonRemovedMuMinus(0)
 202     ,fMultNonRemovedMuPlus(0)
 203     ,fMultNonRemovedGammas(0)
 204     ,fMultNonRemovedNeutrons(0)
 205     ,fMultNonRemovedAntiNeutrons(0)
 206     ,fMultRemovedProtons(0)
 207     ,fMultRemovedAntiProtons(0)
 208     ,fMultRemovedPiPlus(0)
 209     ,fMultRemovedPiMinus(0)
 210     ,fMultRemovedKaonPlus(0)
 211     ,fMultRemovedKaonMinus(0)
 212     ,fMultRemovedK0s(0)
 213     ,fMultRemovedK0L(0)
 214     ,fMultRemovedLambdas(0)
 215     ,fMultRemovedElectrons(0)
 216     ,fMultRemovedPositrons(0)
 217     ,fMultRemovedMuMinus(0)
 218     ,fMultRemovedMuPlus(0)
 219     ,fMultRemovedGammas(0)
 220     ,fMultRemovedNeutrons(0)
 221     ,fMultRemovedAntiNeutrons(0)
 222     ,fTrackMultInAcc(0)
 223     ,fHistDxDzNonRemovedCharged(0)
 224     ,fHistDxDzRemovedCharged(0)
 225     ,fHistDxDzNonRemovedNeutral(0)
 226     ,fHistDxDzRemovedNeutral(0)
 227     ,fHistPiPlusMult(0)
 228     ,fHistPiMinusMult(0)
 229     ,fHistPiZeroMult(0)
 230     ,fHistPiPlusMultAcc(0)
 231     ,fHistPiMinusMultAcc(0)
 232     ,fHistPiZeroMultAcc(0)
 233 //     ,fPiPlusMult(0)
 234 //     ,fPiMinusMult(0)
 235     ,fPiZeroMult(0)
 236     ,fPiPlusMultAcc(0)
 237     ,fPiMinusMultAcc(0)
 238     ,fPiZeroMultAcc(0)
 239     ,fNeutralRemoved(0)
 240     ,fChargedRemoved(0)
 241     ,fChargedNotRemoved(0)
 242     ,fNeutralNotRemoved(0)
 243     ,fGammaRemoved(0)
 244     ,fSecondaryNotRemoved(0)
 245     ,fEnergyNeutralRemoved(0)
 246     ,fEnergyChargedRemoved(0)
 247     ,fEnergyChargedNotRemoved(0)
 248     ,fEnergyNeutralNotRemoved(0)
 249     ,fEnergyGammaRemoved(0)
 250     ,fNClusters(0)
 251     ,fTotNeutralEtAfterMinEnergyCut(0)
 252                                                   ,fCalcTrackMatchVsMult(kFALSE)
 253                                                   ,fHistGammasFound(0)
 254                                                   ,fHistGammasGenerated(0)
 255                                                   ,fHistGammasFoundCent(0)
 256                                                   ,fHistGammasFoundOutOfAccCent(0)
 257                                                   ,fHistGammasFoundAltCent(0)
 258                                                   ,fHistGammasFoundOutOfAccAltCent(0)
 259                                                   ,fHistGammasGeneratedCent(0)
 260                                                   ,fHistGammasFoundRecoEnergyCent(0)
 261                                                   ,fHistAllGammasFoundRecoEnergyCent(0)
 262                                                   ,fHistGammasFoundOutOfAccRecoEnergyCent(0)
 263                                                   ,fHistAllGammasFoundOutOfAccRecoEnergyCent(0)
 264                                                   ,fHistChargedTracksCut(0)
 265                                                   ,fHistChargedTracksAccepted(0)
 266                                                   ,fHistGammasCut(0)
 267                                                   ,fHistGammasAccepted(0)
 268                                                   ,fHistChargedTrackDepositsAcceptedVsPt(0)
 269                                                   ,fHistChargedTrackDepositsAllVsPt(0)
 270                                                   ,fHistChargedTrackDepositsAcceptedVsPtEffCorr(0)
 271                                                   ,fHistChargedTrackDepositsAllVsPtEffCorr(0)
 272                                                   ,fHistChargedTracksCutMult(0)
 273                                                   ,fHistChargedTracksAcceptedMult(0)
 274                                                   ,fHistChargedTracksAcceptedLowPtCentEffCorr(0)
 275                                                   ,fHistChargedTracksAcceptedLowPtCent(0)
 276                                                   ,fHistChargedTracksAcceptedLowPtCent500MeV(0)
 277                                                   ,fHistChargedTracksAcceptedLowPtCentNoAntiProtons(0)
 278                                                   ,fHistChargedTracksAcceptedLowPtCentAntiProtons(0)
 279                                                   ,fHistGammasCutMult(0)
 280                                                   ,fHistGammasAcceptedMult(0)
 281                                                   ,fHistBadTrackMatches(0)
 282                                                   ,fHistMatchedTracksEvspTBkgd(0)
 283                                                   ,fHistMatchedTracksEvspTSignal(0)
 284                                                   ,fHistMatchedTracksEvspTBkgdPeripheral(0)
 285                                                   ,fHistMatchedTracksEvspTSignalPeripheral(0)
 286                                                   ,fHistMatchedTracksEvspTBkgdvsCent(0)
 287                                                   ,fHistMatchedTracksEvspTSignalvsCent(0)
 288                                                   ,fHistMatchedTracksEvspTBkgdvsCentEffCorr(0)
 289                                                   ,fHistMatchedTracksEvspTSignalvsCentEffCorr(0)
 290
 291                                                   ,fHistChargedTracksCutPeripheral(0)
 292                                                   ,fHistChargedTracksAcceptedPeripheral(0)
 293                                                   ,fHistGammasCutPeripheral(0)
 294                                                   ,fHistGammasAcceptedPeripheral(0)
 295                                                   ,fHistBadTrackMatchesdPhidEta(0)
 296                                                   ,fHistGoodTrackMatchesdPhidEta(0)
 297                                                   ,fHistHadronDepositsAll(0)
 298                                                   ,fHistHadronDepositsReco(0)
 299                                                   ,fHistHadronDepositsAllCent(0)
 300                                                   ,fHistHadronDepositsAllCent500MeV(0)
 301                                                   ,fHistHadronDepositsRecoCent(0)
 302                                                   ,fHistHadronDepositsAllvsECent(0)
 303                                                   ,fHistHadronDepositsRecovsECent(0)
 304                                                   ,fHistHadronsAllCent(0)
 305                                                   ,fHistMultChVsSignalVsMult(0)
 306                                                   ,fHistNeutralRemovedSecondaryEtVsCent(0)
 307                                                   ,fHistChargedRemovedSecondaryEtVsCent(0)
 308                                                   ,fHistNeutralNotRemovedSecondaryEtVsCent(0)
 309                                                   ,fHistChargedNotRemovedSecondaryEtVsCent(0)
 310                                                   ,fHistNeutralRemovedSecondaryNumVsNCluster(0)
 311                                                   ,fHistChargedRemovedSecondaryNumVsNCluster(0)
 312                                                   ,fHistNeutralNotRemovedSecondaryNumVsNCluster(0)
 313                                                   ,fHistChargedNotRemovedSecondaryNumVsNCluster(0)
 314                                                   ,fHistNeutralRemovedSecondaryNumVsCent(0)
 315                                                   ,fHistChargedRemovedSecondaryNumVsCent(0)
 316                                                   ,fHistNeutralNotRemovedSecondaryNumVsCent(0)
 317                                                   ,fHistChargedNotRemovedSecondaryNumVsCent(0)
 318                                                   ,fHistNeutronsEtVsCent(0)
 319                                                   ,fHistNeutronsNumVsCent(0)
 320                                                   ,fHistNotNeutronsNumVsCent(0)
 321                                                   ,fHistPiKPDepositedVsNch(0)
 322                                                   ,fHistPiKPNotTrackMatchedDepositedVsNch(0)
 323                                                   ,fHistNeutronsDepositedVsNch(0)//start
 324                                                   ,fHistAntiNeutronsDepositedVsNch(0)
 325                                                   ,fHistProtonsDepositedVsNch(0)
 326                                                   ,fHistAntiProtonsDepositedVsNch(0)
 327                                                   ,fHistProtonsNotTrackMatchedDepositedVsNch(0)
 328                                                   ,fHistAntiProtonsNotTrackMatchedDepositedVsNch(0)
 329                                                   ,fHistNeutronsDepositedVsNcl(0)
 330                                                   ,fHistAntiNeutronsDepositedVsNcl(0)
 331                                                   ,fHistProtonsDepositedVsNcl(0)
 332                                                   ,fHistAntiProtonsDepositedVsNcl(0)
 333                                                   ,fHistProtonsNotTrackMatchedDepositedVsNcl(0)
 334                                                   ,fHistAntiProtonsNotTrackMatchedDepositedVsNcl(0)
 335                                                   ,fHistSecondariesVsNch(0)
 336                                                   ,fHistSecondariesVsNcl(0)
 337                                                   ,fHistSecondariesEffCorrVsNch(0)
 338                                                   ,fHistSecondariesEffCorrVsNcl(0)
 339                                                   ,fHistSecondariesOutOfAccEffCorrVsNch(0)
 340                                                   ,fHistSecondariesDetectorCoverEffCorrVsNch(0)//end
 341                                                   ,fHistNeutronsDepositedVsNchNoEffCorr(0)//start
 342                                                   ,fHistAntiNeutronsDepositedVsNchNoEffCorr(0)
 343                                                   ,fHistProtonsDepositedVsNchNoEffCorr(0)
 344                                                   ,fHistAntiProtonsDepositedVsNchNoEffCorr(0)
 345                                                   ,fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr(0)
 346                                                   ,fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr(0)
 347                                                   ,fHistNeutronsDepositedVsNclNoEffCorr(0)
 348                                                   ,fHistAntiNeutronsDepositedVsNclNoEffCorr(0)
 349                                                   ,fHistProtonsDepositedVsNclNoEffCorr(0)
 350                                                   ,fHistAntiProtonsDepositedVsNclNoEffCorr(0)
 351                                                   ,fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr(0)
 352                                                   ,fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr(0)//end
 353                                                   ,fHistCentVsNchVsNcl(0)
 354                                                         ,fHistSecondaryPositionInDetector(0)
 355                                                   ,fClusterPositionWeird(0)
 356                                                 //,fHistSecondaryPositionInDetectorMultiple(0)
 357                                                   ,fSecondaryClusterEnergy(0)
 358                                                   ,fHistGammaCrossCheck(0)
 359                                                   ,fHistGammaCrossCheckAlt(0)
 360                                                   ,fHistGammaEnergyCrossCheck(0)
 361                                                   ,fHistGammaEnergyCrossCheckAlt(0)
 362     ,fHistNeutronCrossCheck(0)
 363     ,fHistSecondaryCrossCheck(0)
 364     ,fHistHadronCrossCheck(0)
 365     ,fHistKaonCrossCheck(0)
 366     ,fHistAllEnergy(0)
 367     ,fHistSignalEnergy(0)
 368     ,fHistNeutronEnergy(0)
 369     ,fHistKaonEnergy(0)
 370     ,fHistHadronEnergy(0)
 371     ,fHistSecondaryEnergy(0)
 372     ,fHistSecondaryChargedEnergy(0)
 373     ,fHistSecondaryNeutronEnergy(0)
 374     ,fHistSecondaryGammaEnergy(0)
 375     ,fHistSecondaryElectronEnergy(0)
 376     ,fHistSecondaryOtherEnergy(0)
 377     ,fHistSimulatedGammaEnergy(0)
 378     ,fHistReconstructedGammaEnergy(0)
 379                                                   ,fHistSimulatedGammaEnergyAboveThreshold(0)
 380                                                   ,fHistReconstructedSignalEnergy(0)
 381     ,fHistFracSignalVsNClusters(0)
 382     ,fHistFracHadronsVsNClusters(0)
 383     ,fHistFracNeutronsVsNClusters(0)
 384     ,fHistFracKaonsVsNClusters(0)
 385     ,fHistFracSecondariesVsNClusters(0)
 386     ,fHistFracSignalVsNMultiplicity(0)
 387     ,fHistFracHadronsVsNMultiplicity(0)
 388     ,fHistFracNeutronsVsNMultiplicity(0)
 389     ,fHistFracKaonsVsNMultiplicity(0)
 390     ,fHistFracSecondariesVsNMultiplicity(0)
 391     ,fHistFracSignalVsNMatchedTracks(0)
 392     ,fHistFracHadronsVsNMatchedTracks(0)
 393     ,fHistFracNeutronsVsNMatchedTracks(0)
 394     ,fHistFracKaonsVsNMatchedTracks(0)
 395     ,fHistFracSecondariesVsNMatchedTracks(0)
 396     ,fHistFracSignalVsNTotalTracks(0)
 397     ,fHistFracHadronsVsNTotalTracks(0)
 398     ,fHistFracNeutronsVsNTotalTracks(0)
 399     ,fHistFracKaonsVsNTotalTracks(0)
 400     ,fHistFracSecondariesVsNTotalTracks(0)
 401     ,fHistRCorrVsPtVsCent(0)
 402                                                   ,fNMCProducedMin(0)
 403                                                   ,fNMCProducedMax(0)
 404 {
 405 }
 406
 407 // dtor
 408 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
 409 {   //Destructor
 410
 411   if(fPrimaryTree){
 412     fPrimaryTree->Clear();
 413     delete fPrimaryTree;
 414   }
 415     delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
 416     delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
 417     delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
 418     delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
 419
 420     delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
 421
 422     delete fHistEtNonRemovedProtons; // enter comment here
 423     delete fHistEtNonRemovedAntiProtons; // enter comment here
 424     delete fHistEtNonRemovedPiPlus; // enter comment here
 425     delete fHistEtNonRemovedPiMinus; // enter comment here
 426     delete fHistEtNonRemovedKaonPlus; // enter comment here
 427     delete fHistEtNonRemovedKaonMinus; // enter comment here
 428     delete fHistEtNonRemovedK0s; // enter comment here
 429     delete fHistEtNonRemovedK0L; // enter comment here
 430     delete fHistEtNonRemovedLambdas; // enter comment here
 431     delete fHistEtNonRemovedElectrons; // enter comment here
 432     delete fHistEtNonRemovedPositrons; // enter comment here
 433     delete fHistEtNonRemovedMuPlus; // enter comment here
 434     delete fHistEtNonRemovedMuMinus; // enter comment here
 435     delete fHistEtNonRemovedNeutrons; // enter comment here
 436     delete fHistEtNonRemovedAntiNeutrons; // enter comment here
 437     delete fHistEtNonRemovedGammas; // enter comment here
 438     delete fHistEtNonRemovedGammasFromPi0; // enter comment here
 439
 440     delete fHistEtRemovedGammas; // enter comment here
 441     delete fHistEtRemovedNeutrons; // enter comment here
 442     delete fHistEtRemovedAntiNeutrons; // enter comment here
 443
 444
 445     delete fHistMultNonRemovedProtons; // enter comment here
 446     delete fHistMultNonRemovedAntiProtons; // enter comment here
 447     delete fHistMultNonRemovedPiPlus; // enter comment here
 448     delete fHistMultNonRemovedPiMinus; // enter comment here
 449     delete fHistMultNonRemovedKaonPlus; // enter comment here
 450     delete fHistMultNonRemovedKaonMinus; // enter comment here
 451     delete fHistMultNonRemovedK0s; // enter comment here
 452     delete fHistMultNonRemovedK0L; // enter comment here
 453     delete fHistMultNonRemovedLambdas; // enter comment here
 454     delete fHistMultNonRemovedElectrons; // enter comment here
 455     delete fHistMultNonRemovedPositrons; // enter comment here
 456     delete fHistMultNonRemovedMuPlus; // enter comment here
 457     delete fHistMultNonRemovedMuMinus; // enter comment here
 458     delete fHistMultNonRemovedNeutrons; // enter comment here
 459     delete fHistMultNonRemovedAntiNeutrons; // enter comment here
 460     delete fHistMultNonRemovedGammas; // enter comment here
 461
 462     delete fHistMultRemovedGammas; // enter comment here
 463     delete fHistMultRemovedNeutrons; // enter comment here
 464     delete fHistMultRemovedAntiNeutrons; // enter comment here
 465
 466     delete fHistTrackMultvsNonRemovedCharged; // enter comment here
 467     delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
 468     delete fHistTrackMultvsRemovedGamma; // enter comment here
 469
 470     delete fHistClusterMultvsNonRemovedCharged; // enter comment here
 471     delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
 472     delete fHistClusterMultvsRemovedGamma; // enter comment here
 473
 474     delete fHistMultvsNonRemovedChargedE; // enter comment here
 475     delete fHistMultvsNonRemovedNeutralE; // enter comment here
 476     delete fHistMultvsRemovedGammaE; // enter comment here
 477     delete fHistK0EDepositsVsPtInAcceptance;//enter comment here
 478     delete fHistK0EGammaVsPtInAcceptance;//enter comment here
 479     delete fHistK0EDepositsVsPtOutOfAcceptance;
 480     delete fHistK0EGammaVsPtOutOfAcceptance;
 481
 482     delete fHistSimKaonsInAcceptance;
 483     delete fHistSimK0SInAcceptance;
 484     delete fHistSimKPlusInAcceptance;
 485     delete fHistSimKMinusInAcceptance;
 486     delete fHistSimK0LInAcceptance;
 487     delete fHistSimKaonsOutOfAcceptance;
 488     delete fHistSimKaonsInAcceptanceWithDepositsPrimaries;
 489     delete fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries;
 490     delete fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries;
 491
 492     delete fHistDxDzNonRemovedCharged; // enter comment here
 493     delete fHistDxDzRemovedCharged; // enter comment here
 494     delete fHistDxDzNonRemovedNeutral; // enter comment here
 495     delete fHistDxDzRemovedNeutral; // enter comment here
 496
 497     delete fHistPiPlusMult; // enter comment here
 498     delete fHistPiMinusMult; // enter comment here
 499     delete fHistPiZeroMult; // enter comment here
 500
 501     delete fHistPiPlusMultAcc; // enter comment here
 502     delete fHistPiMinusMultAcc; // enter comment here
 503     delete fHistPiZeroMultAcc; // enter comment here
 504     delete fHistGammasFound; // enter comment here
 505     delete fHistGammasGenerated; // enter comment here
 506     delete fHistGammasFoundOutOfAccCent; // enter comment here
 507     delete fHistGammasFoundCent; // enter comment here
 508     delete fHistGammasFoundOutOfAccAltCent; // enter comment here
 509     delete fHistGammasFoundAltCent; // enter comment here
 510     delete fHistGammasGeneratedCent; // enter comment here
 511     delete fHistGammasFoundRecoEnergyCent;
 512     delete fHistAllGammasFoundRecoEnergyCent;
 513     delete fHistGammasFoundOutOfAccRecoEnergyCent;
 514     delete fHistAllGammasFoundOutOfAccRecoEnergyCent;
 515     delete fHistChargedTracksCut;
 516     delete fHistChargedTracksAccepted;
 517     delete fHistGammasCut;
 518     delete fHistGammasAccepted;
 519     delete fHistChargedTracksCutMult;
 520     delete fHistChargedTrackDepositsAcceptedVsPt;
 521     delete fHistChargedTrackDepositsAllVsPt;
 522     delete fHistChargedTrackDepositsAcceptedVsPtEffCorr;
 523     delete fHistChargedTrackDepositsAllVsPtEffCorr;
 524     delete fHistChargedTracksAcceptedMult;
 525     delete fHistChargedTracksAcceptedLowPtCentEffCorr;
 526     delete fHistChargedTracksAcceptedLowPtCent;
 527     delete fHistChargedTracksAcceptedLowPtCent500MeV;
 528     delete fHistChargedTracksAcceptedLowPtCentNoAntiProtons;
 529     delete fHistChargedTracksAcceptedLowPtCentAntiProtons;
 530     delete fHistGammasCutMult;
 531     delete fHistGammasAcceptedMult;
 532     delete fHistBadTrackMatches;
 533     delete fHistMatchedTracksEvspTBkgd;
 534     delete fHistMatchedTracksEvspTSignal;
 535     delete fHistMatchedTracksEvspTBkgdPeripheral;
 536     delete fHistMatchedTracksEvspTSignalPeripheral;
 537     delete fHistMatchedTracksEvspTBkgdvsCent;
 538     delete fHistMatchedTracksEvspTSignalvsCent;
 539     delete fHistMatchedTracksEvspTBkgdvsCentEffCorr;
 540     delete fHistMatchedTracksEvspTSignalvsCentEffCorr;
 541     delete fHistChargedTracksCutPeripheral;
 542     delete fHistChargedTracksAcceptedPeripheral;
 543     delete fHistGammasCutPeripheral;
 544     delete fHistGammasAcceptedPeripheral;
 545     delete fHistBadTrackMatchesdPhidEta;
 546     delete fHistGoodTrackMatchesdPhidEta;
 547     delete fHistHadronDepositsAll;
 548     delete fHistHadronDepositsReco;
 549     delete fHistHadronDepositsAllCent;
 550     delete fHistHadronDepositsAllCent500MeV;
 551     delete fHistHadronDepositsRecoCent;
 552     delete fHistHadronDepositsAllvsECent;
 553     delete fHistHadronDepositsRecovsECent;
 554     delete fHistHadronsAllCent;
 555     delete fHistMultChVsSignalVsMult;
 556     delete fHistNeutralRemovedSecondaryEtVsCent;
 557     delete fHistChargedRemovedSecondaryEtVsCent;
 558     delete fHistNeutralNotRemovedSecondaryEtVsCent;
 559     delete fHistChargedNotRemovedSecondaryEtVsCent;
 560     delete fHistNeutralRemovedSecondaryNumVsNCluster;
 561     delete fHistChargedRemovedSecondaryNumVsNCluster;
 562     delete fHistNeutralNotRemovedSecondaryNumVsNCluster;
 563     delete fHistChargedNotRemovedSecondaryNumVsNCluster;
 564     delete fHistNeutralRemovedSecondaryNumVsCent;
 565     delete fHistChargedRemovedSecondaryNumVsCent;
 566     delete fHistNeutralNotRemovedSecondaryNumVsCent;
 567     delete fHistChargedNotRemovedSecondaryNumVsCent;
 568     delete fHistNeutronsEtVsCent;
 569     delete fHistNeutronsNumVsCent;
 570     delete fHistNotNeutronsNumVsCent;
 571     delete fHistPiKPDepositedVsNch;
 572     delete fHistPiKPNotTrackMatchedDepositedVsNch;
 573
 574     delete fHistNeutronsDepositedVsNch;
 575     delete fHistAntiNeutronsDepositedVsNch;
 576     delete fHistProtonsDepositedVsNch;
 577     delete fHistAntiProtonsDepositedVsNch;
 578     delete fHistProtonsNotTrackMatchedDepositedVsNch;
 579     delete fHistAntiProtonsNotTrackMatchedDepositedVsNch;
 580     delete fHistNeutronsDepositedVsNcl;
 581     delete fHistAntiNeutronsDepositedVsNcl;
 582     delete fHistProtonsDepositedVsNcl;
 583     delete fHistAntiProtonsDepositedVsNcl;
 584     delete fHistProtonsNotTrackMatchedDepositedVsNcl;
 585     delete fHistAntiProtonsNotTrackMatchedDepositedVsNcl;
 586     delete fHistSecondariesVsNch;
 587     delete fHistSecondariesVsNcl;
 588     delete fHistSecondariesEffCorrVsNch;
 589     delete fHistSecondariesEffCorrVsNcl;
 590     delete fHistSecondariesOutOfAccEffCorrVsNch;
 591     delete fHistSecondariesDetectorCoverEffCorrVsNch;
 592
 593
 594     delete fHistNeutronsDepositedVsNchNoEffCorr;
 595     delete fHistAntiNeutronsDepositedVsNchNoEffCorr;
 596     delete fHistProtonsDepositedVsNchNoEffCorr;
 597     delete fHistAntiProtonsDepositedVsNchNoEffCorr;
 598     delete fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr;
 599     delete fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr;
 600     delete fHistNeutronsDepositedVsNclNoEffCorr;
 601     delete fHistAntiNeutronsDepositedVsNclNoEffCorr;
 602     delete fHistProtonsDepositedVsNclNoEffCorr;
 603     delete fHistAntiProtonsDepositedVsNclNoEffCorr;
 604     delete fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr;
 605     delete fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr;
 606
 607     delete fHistCentVsNchVsNcl;
 608     delete fHistSecondaryPositionInDetector;
 609     delete fClusterPositionWeird;
 610     //delete fHistSecondaryPositionInDetectorMultiple;
 611     delete fSecondaryClusterEnergy;
 612     delete fHistGammaCrossCheck;
 613     delete fHistGammaCrossCheckAlt;
 614     delete fHistGammaEnergyCrossCheck;
 615     delete fHistGammaEnergyCrossCheckAlt;
 616     delete fHistNeutronCrossCheck;
 617     delete fHistSecondaryCrossCheck;
 618     delete fHistHadronCrossCheck;
 619     delete fHistKaonCrossCheck;
 620
 621     delete fHistAllEnergy;
 622     delete fHistSignalEnergy;
 623     delete fHistNeutronEnergy;
 624     delete fHistKaonEnergy;
 625     delete fHistHadronEnergy;
 626     delete fHistSecondaryEnergy;
 627     delete fHistSecondaryChargedEnergy;
 628     delete fHistSecondaryNeutronEnergy;
 629     delete fHistSecondaryGammaEnergy;
 630     delete fHistSecondaryElectronEnergy;
 631     delete fHistSecondaryOtherEnergy;
 632     delete fHistSimulatedGammaEnergy;
 633     delete fHistReconstructedGammaEnergy;
 634     delete fHistSimulatedGammaEnergyAboveThreshold;
 635     delete fHistReconstructedSignalEnergy;
 636     delete fHistFracSignalVsNClusters;
 637     delete fHistFracHadronsVsNClusters;
 638     delete fHistFracNeutronsVsNClusters;
 639     delete fHistFracKaonsVsNClusters;
 640     delete fHistFracSecondariesVsNClusters;
 641     delete fHistFracSignalVsNMultiplicity;
 642     delete fHistFracHadronsVsNMultiplicity;
 643     delete fHistFracNeutronsVsNMultiplicity;
 644     delete fHistFracKaonsVsNMultiplicity;
 645     delete fHistFracSecondariesVsNMultiplicity;
 646     delete fHistFracSignalVsNMatchedTracks;
 647     delete fHistFracHadronsVsNMatchedTracks;
 648     delete fHistFracNeutronsVsNMatchedTracks;
 649     delete fHistFracKaonsVsNMatchedTracks;
 650     delete fHistFracSecondariesVsNMatchedTracks;
 651     delete fHistFracSignalVsNTotalTracks;
 652     delete fHistFracHadronsVsNTotalTracks;
 653     delete fHistFracNeutronsVsNTotalTracks;
 654     delete fHistFracKaonsVsNTotalTracks;
 655     delete fHistFracSecondariesVsNTotalTracks;
 656     delete fHistRCorrVsPtVsCent;
 657 }
 658
 659 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 660 {   // analyse MC event
 661   if(!fIsMC) return 0;
 662     ResetEventValues();
 663
 664
 665     fPiPlusMult = 0;
 666     fPiMinusMult = 0;
 667     fPiZeroMult = 0;
 668     if (fCentrality)
 669     {
 670         fCentClass = fCentrality->GetCentralityClass5(fCentralityMethod);
 671
 672     }
 673
 674     // Get us an mc event
 675     if (!ev) {
 676         AliFatal("ERROR: Event does not exist");
 677         return 0;
 678     }
 679     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
 680     if (!event) {
 681         AliFatal("ERROR: MC Event does not exist");
 682         return 0;
 683     }
 684
 685     Double_t protonMass =fgProtonMass;
 686
 687     // Hijing header
 688     AliGenEventHeader* genHeader = event->GenEventHeader();
 689     if (!genHeader) {
 690         Printf("ERROR: Event generation header does not exist");
 691         return 0;
 692     }
 693     AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 694     if (hijingGenHeader) {
 695         fImpactParameter = hijingGenHeader->ImpactParameter();
 696         fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
 697         fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
 698     }
 699
 700     // Let's play with the stack!
 701     AliStack *stack = event->Stack();
 702
 703     Int_t nPrim = stack->GetNtrack();
 704     //cout<<endl<<endl<<"new event simulated nPrim "<<nPrim<<endl;
 705
 706     Int_t partCount = 0;
 707     for (Int_t iPart = 0; iPart < nPrim; iPart++)
 708     {
 709
 710         //Some productions have signals added in.  This switch allows the explicit exclusion of these added in signals when running over the data.
 711         if(checkLabelForHIJING && !IsHIJINGLabel(iPart,event,stack) ) continue;
 712         TParticle *part = stack->Particle(iPart);
 713
 714
 715
 716         if (!part)
 717         {
 718           Printf("ERROR: Could not get particle %d", iPart);
 719             continue;
 720         }
 721         TParticlePDG *pdg = part->GetPDG(0);
 722
 723         if (!pdg)
 724         {
 725           //Printf("ERROR: Could not get particle PDG %d", iPart);
 726             continue;
 727         }
 728
 729         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
 730         Int_t code = pdg->PdgCode();
 731
 732         if(stack->IsPhysicalPrimary(iPart))
 733         {
 734           fTotPx += part->Px();
 735           fTotPy += part->Py();
 736           fTotPz += part->Pz();
 737         }
 738
 739         // Check for reasonable (for now neutral and singly charged) charge on the particle
 740         if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
 741         {
 742
 743             fMultiplicity++;
 744              //PrintFamilyTreeShort(iPart, stack);
 745 //
 746             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 747             {
 748                 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
 749                 // calculate E_T
 750
 751               //fHistGammasFoundRecoEnergyCent->Fill(fReconstructedE,fCentClass);
 752               //cout<<"Messed up filling would have filled with "<<fReconstructedE<<" cent "<<fCentClass<<endl;
 753                   if (
 754                     TMath::Abs(code) == fgProtonCode ||
 755                     TMath::Abs(code) == fgNeutronCode ||
 756                     TMath::Abs(code) == fgLambdaCode ||
 757                     TMath::Abs(code) == fgXiCode ||
 758                     TMath::Abs(code) == fgXi0Code ||
 759                     TMath::Abs(code) == fgOmegaCode
 760                 )
 761                 {
 762                     if (code > 0) {
 763                         particleMassPart = - protonMass;
 764                     }
 765                     if (code < 0) {
 766                         particleMassPart = protonMass;
 767                     }
 768                 }
 769                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 770
 771
 772                 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode || code==fgOmega0Code || code==fgEPlusCode || code==fgEMinusCode)
 773                 {
 774
 775                     if(fSelector->CutGeometricalAcceptance(*part) )
 776                     {
 777                         fNeutralMultiplicity++;
 778                         fTotNeutralEt += et;
 779                         if(fSelector->PassMinEnergyCut(*part))
 780                         {
 781                             fTotNeutralEtAfterMinEnergyCut += et;
 782                         }
 783                         if (part->Energy() > 0.05) partCount++;
 784                     }
 785                 }
 786                 //Charged particles
 787                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
 788                 {
 789
 790                     // inside EMCal acceptance
 791                     if (fSelector->CutGeometricalAcceptance(*part))
 792                     {
 793
 794                         fChargedMultiplicity++;
 795
 796                         fTotChargedEt += et;
 797
 798                         if (code == fgProtonCode || code == fgAntiProtonCode)
 799                         {
 800                         }
 801                         if (code == fgPiPlusCode || code == fgPiMinusCode)
 802                         {
 803                         }
 804                         if (code == fgKPlusCode || code == fgKMinusCode)
 805                         {
 806                         }
 807                         if (code == fgMuPlusCode || code == fgMuMinusCode)
 808                         {
 809                         }
 810
 811                         if (code == fgEPlusCode || code == fgEMinusCode)
 812                         {
 813                             fTotNeutralEt += et; // calling electrons neutral
 814                             fTotChargedEt -= et;
 815                         }
 816                     } // inside EMCal acceptance
 817
 818                     if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 819                     {
 820                         if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 821                         else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 822                     }
 823                 }
 824             }
 825
 826         }
 827     }
 828    // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
 829     fTotEt = fTotChargedEt + fTotNeutralEt;
 830     //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
 831     //std::cout << "Event done! # of particles: " << partCount << std::endl;
 832     return 0;
 833 }
 834 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 835 {   // analyse MC and real event info
 836     //if(!mcEvent || !realEvent){
 837   if(!fIsMC) return 0;
 838     if (!ev || !ev2) {
 839         AliFatal("ERROR: Event does not exist");
 840         return 0;
 841     }
 842     AliAnalysisEt::AnalyseEvent(ev);
 843     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 844     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
 845     if (!mcEvent || !realEvent) {
 846         AliFatal("ERROR: mcEvent or realEvent does not exist");
 847
 848     }
 849     if(checkLabelForHIJING) SetGeneratorMinMaxParticles(mcEvent);
 850     std::vector<Int_t> foundGammas;
 851
 852     fSelector->SetEvent(realEvent);
 853
 854     AnalyseEvent(ev);
 855
 856     AliStack *stack = mcEvent->Stack();
 857
 858     TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(realEvent);
 859     Int_t nGoodTracks = list->GetEntries();
 860
 861     //cout<<"fcuts max phi "<<fCuts->GetGeometryEmcalPhiAccMaxCut()<<endl;
 862     //Note that this only returns clusters for the selected detector.  fSelector actually calls the right GetClusters... for the detector
 863     //It does not apply any cuts on these clusters
 864     TRefArray *caloClusters = fSelector->GetClusters();
 865
 866     Int_t nCluster = caloClusters->GetEntries();
 867     fClusterMult = nCluster;
 868
 869     //cout<<endl<<"new event reconstructed nclusters "<<nCluster<<endl;
 870     fNClusters = 0;
 871     Int_t fClusterMultChargedTracks = 0;
 872     Int_t fClusterMultGammas = 0;
 873     Int_t nChargedSecondariesRemoved = 0;
 874     Int_t nChargedSecondariesNotRemoved = 0;
 875     Int_t nNeutralSecondariesRemoved = 0;
 876     Int_t nNeutralSecondariesNotRemoved = 0;
 877     Int_t nNeutrons = 0;
 878     Int_t nNotNeutrons = 0;
 879     //All efficiency corrected except etSecondaries -->
 880     Float_t etPiKPDeposited = 0.0;//
 881     Float_t etPiKPDepositedNotTrackMatched = 0.0;//
 882     Float_t etProtonDepositedNotTrackMatched = 0.0;//
 883     Float_t etAntiProtonDepositedNotTrackMatched = 0.0;//
 884     Float_t etProtonDepositedNotTrackMatchedNoEffCorr = 0.0;//
 885     Float_t etAntiProtonDepositedNotTrackMatchedNoEffCorr = 0.0;//
 886     //Float_t etPIDProtonDepositedNotTrackMatched = 0.0;//Still has to be filled!
 887     //1Float_t etPIDAntiProtonDepositedNotTrackMatched = 0.0;//Still has to be filled
 888     Float_t etProtonDeposited = 0.0;//
 889     Float_t etAntiProtonDeposited = 0.0;//
 890     Float_t etNeutronDeposited = 0.0;
 891     Float_t etAntiNeutronDeposited = 0.0;
 892     Float_t etSecondaries = 0.0;
 893     Float_t etSecondariesEffCorr = 0.0;
 894     Float_t etSecondariesOutOfAccEffCorr = 0.0;
 895     Float_t etSecondariesDetectorCoverEffCorr = 0.0;
 896     Float_t etProtonDepositedNoEffCorr = 0.0;//
 897     Float_t etAntiProtonDepositedNoEffCorr = 0.0;//
 898     Float_t etNeutronDepositedNoEffCorr = 0.0;
 899     Float_t etAntiNeutronDepositedNoEffCorr = 0.0;
 900     Float_t etGammaCrossCheck = 0.0;
 901     Float_t etGammaCrossCheckAlt = 0.0;
 902 //     Float_t etNeutronCrossCheck = 0.0;
 903 //     Float_t etSecondaryCrossCheck = 0.0;
 904 //     Float_t etHadronCrossCheck = 0.0;
 905 //     Float_t etKaonCrossCheck = 0.0;
 906     Float_t etGammaCrossCheckTrue = 0.0;//--
 907     Float_t etNeutronCrossCheckTrue = 0.0;//--
 908     Float_t etSecondaryCrossCheckTrue = 0.0;//--
 909     Float_t etHadronCrossCheckTrue = 0.0;//--
 910     Float_t etKaonCrossCheckTrue = 0.0;//
 911     Float_t multiplicity = fEsdtrackCutsTPC->GetReferenceMultiplicity(realEvent,kTRUE);
 912
 913     Float_t subtotalAllEnergy = 0.0;//energy of all clusters passing cuts vs centrality
 914     Float_t subtotalSignalEnergy = 0.0;//signal of signal clusters passing cuts vs centrality
 915     Float_t subtotalNeutronEnergy = 0.0;//signal of neutron clusters passing cuts vs centrality
 916     Float_t subtotalKaonEnergy = 0.0;//signal of kaon clusters passing cuts vs centrality
 917     Float_t subtotalHadronEnergy = 0.0;//signal of hadron clusters passing cuts vs centrality
 918     Float_t subtotalSecondaryEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 919     Float_t subtotalSecondaryChargedEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 920     Float_t subtotalSecondaryNeutronEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 921     Float_t subtotalSecondaryGammaEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 922     Float_t subtotalSecondaryElectronEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 923     Float_t subtotalSecondaryOtherEnergy = 0.0;//signal of secondary clusters passing cuts vs centrality
 924     Float_t subtotalSimulatedGammaEnergy = 0.0;//signal of signal clusters passing cuts vs centrality
 925     Float_t subtotalSimulatedGammaEnergyAboveThreshold = 0.0;//signal of signal clusters passing cuts vs centrality
 926     Float_t subtotalReconstructedGammaEnergy = 0.0;//signal of signal clusters passing cuts vs centrality
 927     Float_t subtotalReconstructedSignalEnergy = 0.0;//signal of signal clusters passing cuts vs centrality
 928     // cout<<"fsub "<<fsub<<endl;
 929
 930     // loop the clusters
 931     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 932     {
 933         Int_t cf = 0;
 934         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 935         //Float_t caloE = caloCluster->E()
 936         if (!fSelector->CutGeometricalAcceptance(*caloCluster)) continue;
 937         fNClusters++;
 938         //const UInt_t iPart = (UInt_t)TMath::Abs(fSelector->GetLabel(caloCluster));//->GetLabel());
 939         const UInt_t iPart = fSelector->GetLabel(caloCluster,stack);
 940         //if(checkLabelForHIJING) cerr<<"I am checking the label"<<endl;
 941         //Some productions have signals added in.  This switch allows the explicit exclusion of these added in signals when running over the data.
 942         if(checkLabelForHIJING && !IsHIJINGLabel(iPart,mcEvent,stack) ) continue;
 943         //const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
 944         TParticle *part  =  stack->Particle(iPart);
 945
 946         if (!part)
 947         {
 948             Printf("No MC particle %d", iCluster);
 949             continue;
 950         }
 951
 952         int primIdx = iPart;
 953         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
 954         {
 955             primIdx = GetPrimMother(iPart, stack);
 956             if(primIdx != stack->GetPrimary(iPart))
 957             {
 958               //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
 959               //PrintFamilyTree(iPart, stack);
 960             }
 961             //if it is from a K0S
 962             if(primIdx < 0)
 963             {
 964               //std::cout << "What!? No primary?" << std::endl;
 965               //PrintFamilyTree(iPart, stack);
 966               //continue;
 967               //This is a work around to fix a bug.  For the EMCal when you use the tender supply, the parent particle ID gets messed up.
 968               primIdx = iPart;
 969             }
 970
 971         } // end of primary particle check
 972         //const int primCode = stack->Particle(primIdx)->GetPdgCode();
 973         TParticlePDG *pdg = part->GetPDG();
 974         if (!pdg)
 975         {
 976             Printf("ERROR: Could not get particle PDG %d", iPart);
 977             continue;
 978         }
 979
 980         //Int_t code = pdg->PdgCode();
 981 //      if(primCode == fgGammaCode)
 982 //      {
 983
 984
 985
 986         Bool_t nottrackmatched = kTRUE;//default to no track matched
 987         Bool_t countasmatched = kFALSE;
 988         Float_t matchedTrackp = 0.0;
 989         Float_t matchedTrackpt = 0.0;
 990         fDepositedCode = part->GetPdgCode();
 991         fReconstructedE = caloCluster->E();
 992         Float_t pos[3];
 993         //PrintFamilyTree(
 994         caloCluster->GetPosition(pos);
 995         TVector3 cp(pos);
 996         fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
 997         nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
 998         //by default ALL matched tracks are accepted, whether or not the match is good.  So we check to see if the track is good.
 999         if(!nottrackmatched){//if the track is trackmatched
1000           Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
1001           if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
1002           AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
1003           // cout<<"track code "<<fTrackDepositedCode<<" cluster code "<<fDepositedCode<<" track label "<<trackLabel<<" cluster label "<<iPart<<endl;
1004
1005           //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
1006           nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));//if the track is bad, this is not track matched
1007           if(!nottrackmatched){//if the track that was matched is a good track
1008             matchedTrackp = track->P();
1009             matchedTrackpt = track->Pt();
1010             UInt_t trackLabel = (UInt_t)TMath::Abs(track->GetLabel());
1011             if(checkLabelForHIJING && !IsHIJINGLabel(trackLabel,mcEvent,stack) ) nottrackmatched = kTRUE;//if this was matched to something we don't want to count
1012             TParticle  *trackSimPart  = stack->Particle(trackLabel);
1013             Int_t fTrackDepositedCode = trackSimPart->GetPdgCode();
1014             //if(!nottrackmatched) cout<<"Matched track p: "<<matchedTrackp<<" sim "<<part->P()<<endl;
1015             //now we want to fill our Rcorr histos
1016             Float_t rcorr = fReconstructedE/track->P();
1017             fHistRCorrVsPtVsCent->Fill(rcorr,matchedTrackpt, fCentClass);
1018             //cout<<"rcorr "<<rcorr<<endl;
1019             Int_t n=caloCluster->GetNLabels() ;
1020             //if(fReconstructedE - fsub* track->P() > 0.0){//if more energy was deposited than the momentum of the track  and more than one particle led to the cluster
1021             //cout<<"was matched"<<endl;
1022             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
1023               //then we say the cluster was not track matched but correct the energy
1024               nottrackmatched = kTRUE;
1025               //but we don't want to double count the matched tracks...
1026               countasmatched = kTRUE;
1027               //cout<<" fsub "<<fsub;
1028               //cout<<" Reassigning energy "<<fReconstructedEt;
1029               fReconstructedE = fReconstructedE - fsub* track->P();
1030               fReconstructedEt = fReconstructedE*TMath::Sin(cp.Theta());
1031               //cout<<" to "<<fReconstructedEt<<endl;
1032               if(fDepositedCode==fTrackDepositedCode && n>1){
1033                 //the energy deposited was more than the momentum of the track but the cluster was assigned to the particle that created the track.  We therefore need to reassign the cluster label
1034                 //this is a sub-optimal way to re-assign the particle but it is rare that a particle needs to be reassigned
1035                 //cout<<"Particle was "<<part->GetName()<<" reassigning label from "<<fDepositedCode<<", testing ";
1036                 Int_t iMax=-1;
1037                 Double_t*  Ekin=  new  Double_t[n] ;
1038                 for(Int_t i=0;  i<n;  i++){
1039                   TParticle*  p=  stack->Particle(caloCluster->GetLabelAt(i)) ;
1040                   if(p->GetPdgCode()==fgPiPlusCode || p->GetPdgCode()==fgKPlusCode || p->GetPdgCode()==fgProtonCode || p->GetPdgCode()==fgPiMinusCode || p->GetPdgCode()==fgKMinusCode || p->GetPdgCode()==fgAntiProtonCode){
1041                     Ekin[i]=0.3;//p->P()/3.0 ;  // estimate of MIP peak, more likely to be what was deposited if the MIP peak overlapped with another peak
1042                   }
1043                   else{
1044                     Ekin[i]=p->Energy() ;  // what's deposited by electromagnetic particles
1045                   }
1046                   if(p->GetPdgCode()==fgAntiProtonCode){
1047                     Ekin[i]+=1.8  ;  //due to annihilation
1048                   }
1049                 }
1050                 Double_t eMax=0.;//eSubMax=0. ;
1051                 //cout<<"n="<<n<<", ";
1052                 for(Int_t i=0;  i<n;  i++){
1053                     Int_t label = caloCluster->GetLabelAt(i);
1054                     if(label!=trackMatchedIndex){
1055                       //TParticle*  p=  stack->Particle(caloCluster->GetLabelAt(i)) ;
1056                       //cout<<i<<" "<<p->GetName()<<" with E="<<Ekin[i]<<",";
1057                       if(Ekin[i]>eMax){
1058                         //      eSubMax=eMax;
1059                         eMax=Ekin[i];
1060                         iMax=i;
1061                       }
1062                     }
1063                 }
1064                 delete [] Ekin;
1065                 UInt_t correctLabel = caloCluster->GetLabelAt(iMax);
1066                 if(iMax>0){
1067                   TParticle *newPart  =  stack->Particle(correctLabel);
1068                   if(newPart){
1069                     part = newPart;
1070                     fDepositedCode = part->GetPdgCode();
1071                     //cout<<", to "<<fDepositedCode<<" and particle is now "<<part->GetName();
1072                   }
1073                 }
1074                 //cout<<endl;
1075               }
1076               //              if(fDepositedCode==fgProtonCode ||  fDepositedCode==fgAntiProtonCode || fDepositedCode==fgPiPlusCode || fDepositedCode==fgPiMinusCode || fDepositedCode==fgKPlusCode || fDepositedCode==fgKMinusCode){
1077               //        for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
1078               //        {
1079               //          Int_t pIdx = caloCluster->GetLabelAt(i);
1080               //        const UInt_t iPart = fSelector->GetLabel(caloCluster,stack);
1081               //        //const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
1082               //         TParticle *part  =  stack->Particle(iPart);
1083
1084
1085             }
1086
1087           }
1088         }
1089
1090         Bool_t containsGamma = kFALSE;
1091         //Int_t gammaCode = -1;
1092         Float_t gammaEnergy = 0.0;
1093         for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
1094         {
1095           Int_t pIdx = caloCluster->GetLabelAt(i);
1096           //Some productions have signals added in.  This switch allows the explicit exclusion of these added in signals when running over the data.
1097           if(checkLabelForHIJING && !IsHIJINGLabel(pIdx,mcEvent,stack) ) continue;
1098           //Int_t initialLabel = pIdx;
1099           //TParticle *p = stack->Particle(pIdx);
1100
1101           if(!stack->IsPhysicalPrimary(pIdx))
1102             {//This should count gammas as reconstructed if we found even part of them
1103               pIdx = GetPrimMother(pIdx, stack);
1104               //TParticle *part2  =  stack->Particle(initialLabel);
1105 //            TParticle *part2  =  stack->Particle(pIdx);
1106 //            if(part2){
1107 //              TParticlePDG *pdg2 = part2->GetPDG();
1108 //              if(pdg2){
1109 //                Int_t code2 = pdg2->PdgCode();
1110 //                if(code2 == fgGammaCode){
1111 //                  cout<<"Relabeling primary index initial label "<<initialLabel<<" final label "<<pIdx<<endl;
1112 //                  PrintFamilyTree(initialLabel, stack);
1113 //                }
1114 //              }
1115 //            }
1116             }
1117           if(fSelector->PassDistanceToBadChannelCut(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
1118           {
1119             //if this contained a gamma...
1120             if(pIdx>0){
1121               TParticle *part2  =  stack->Particle(pIdx);
1122               if(part2){
1123                 TParticlePDG *pdg2 = part2->GetPDG();
1124                 if(pdg2){
1125                   Int_t code2 = pdg2->PdgCode();
1126                   if(code2 == fgGammaCode){//if this is counted as a reconstructed gamma, we want to know what was reconstructed about it
1127                     //gammaCode = pIdx;
1128                     containsGamma = kTRUE;
1129                     gammaEnergy += part2->Energy();//do += because if we have multiple gammas we want to add that to our reconstructed gamma energy, even if we didn't separat them
1130                   }
1131                 }
1132               }
1133               //            std::cout << "Gamma primary: " << primIdx << std::endl;
1134               //            foundGammas.push_back(primIdx);
1135               if(nottrackmatched){
1136                 foundGammas.push_back(pIdx);
1137               }
1138             }
1139           }
1140         }
1141         fCutFlow->Fill(cf++);
1142         if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
1143         Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster,fReconstructedE,fCentClass);
1144 //      if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
1145         //if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
1146                  if(!fSelector->PassMinEnergyCut(fReconstructedE)) continue;
1147
1148
1149         fCutFlow->Fill(cf++);
1150
1151         TParticle *primPart = stack->Particle(primIdx);
1152         fPrimaryCode = primPart->GetPdgCode();
1153         if(primPart->GetPDG()) fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
1154
1155         fPrimaryE = primPart->Energy();
1156         fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
1157         fPrimaryPx = primPart->Px();
1158         fPrimaryPy = primPart->Py();
1159         fPrimaryPz = primPart->Pz();
1160         //cout<<"I have a cluster and it's good energy "<<caloCluster->E()<<" simulated "<<fPrimaryE<<endl;
1161         fPrimaryVx = primPart->Vx();
1162         fPrimaryVy = primPart->Vy();
1163         fPrimaryVz = primPart->Vz();
1164
1165         fPrimaryAccepted = false;
1166         fPrimaryMatched = false;
1167
1168         fDepositedE = part->Energy();
1169         fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
1170         if(part->GetPDG()) fDepositedCharge = (Int_t) part->GetPDG()->Charge();
1171
1172         fDepositedVx = part->Vx();
1173         fDepositedVy = part->Vy();
1174         fDepositedVz = part->Vz();
1175
1176         if(nottrackmatched) subtotalAllEnergy += fReconstructedEt;
1177
1178         //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
1179         fSecondary =fSelector->FromSecondaryInteraction(iPart, *stack);
1180         //if(fSecondary && fReconstructedEt<0.3) fSecondary = kFALSE;//patch to do quick cross check THIS SHOULD NOT GET COMMITTED!!!  IF IT DID YELL AT CHRISTINE ASAP
1181 //      if(fSecondary)
1182 //      {
1183 //        //std::cout << "Have secondary!" << std::endl;
1184 //        //PrintFamilyTree(iPart, stack);
1185 //      }
1186
1187         pdg = primPart->GetPDG(0);
1188         if (!pdg)
1189         {
1190           //Printf("ERROR: Could not get particle PDG %d", iPart);
1191             continue;
1192         }
1193         //Int_t code = primPart->GetPdgCode();
1194
1195         Bool_t written = kFALSE;
1196
1197 //      Bool_t nottrackmatched = kTRUE;//default to no track matched
1198 //      nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
1199 //      //by default ALL matched tracks are accepted, whether or not the match is good.  So we check to see if the track is good.
1200 //      if(!nottrackmatched){
1201 //        Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
1202 //        if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
1203 //        AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
1204 //        //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
1205 //        nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));
1206 //      }
1207
1208         if(fSecondary){//all particles from secondary interactions
1209           //===================================BEGIN SECONDARIES==================================================
1210           written = kTRUE;
1211
1212           if (fDepositedCharge != 0 && fDepositedCode!=fgEMinusCode && fDepositedCode!=fgEPlusCode){//if the particle hitting the calorimeter is pi/k/p/mu
1213           //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BEGIN CHARGED SECONDARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214             fHistChargedTrackDepositsAllVsPt->Fill(part->Pt(), fCentClass,fReconstructedEt);
1215             fHistChargedTrackDepositsAllVsPtEffCorr->Fill(part->Pt(), fCentClass,clEt);
1216             if(!nottrackmatched){//if it is track matched
1217 //            cout<<"Charged Secondary ";
1218 //            PrintFamilyTreeShort(iPart, stack);
1219               fHistChargedTrackDepositsAcceptedVsPt->Fill(part->Pt(), fCentClass,fReconstructedEt);
1220               fHistChargedTrackDepositsAcceptedVsPtEffCorr->Fill(part->Pt(), fCentClass,fReconstructedEt);
1221               fHistHadronDepositsReco->Fill(part->Pt());
1222               fHistHadronDepositsRecoCent->Fill(part->Pt(), fCentClass);
1223               fHistHadronDepositsRecovsECent->Fill(fReconstructedE, fCentClass);
1224               fChargedRemoved++;
1225               fEnergyChargedRemoved += clEt;
1226               fHistRemovedOrNot->Fill(0.0, fCentClass);
1227               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1228               fHistChargedRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1229               nChargedSecondariesRemoved++;
1230               fHistChargedTracksCut->Fill(fReconstructedEt);
1231               if(fCalcTrackMatchVsMult){
1232                 fHistChargedTracksCutMult->Fill(fReconstructedEt,fClusterMult);
1233                 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fReconstructedEt);}
1234               }
1235               fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
1236               if(fCalcTrackMatchVsMult){
1237                   if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1238                   fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1239                   fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//Fill with the efficiency corrected energy
1240               }
1241               //Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
1242               UInt_t trackindex = fSelector->GetLabel(caloCluster,stack);//(caloCluster->GetLabelsArray())->At(1);
1243               if(((UInt_t)caloCluster->GetLabel())!=trackindex){
1244                 //if(fSelector->GetLabel(caloCluster,stack) !=trackindex){
1245                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1246                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1247                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1248               }
1249               else{
1250                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1251               }
1252             }
1253             else{
1254               fChargedNotRemoved++;
1255               fEnergyChargedNotRemoved += clEt;
1256               fHistRemovedOrNot->Fill(2.0, fCentClass);
1257               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1258               fHistChargedNotRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1259               nChargedSecondariesNotRemoved++;
1260             }
1261             fHistHadronDepositsAll->Fill(part->Pt());
1262             fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
1263             fHistHadronDepositsAllvsECent->Fill(fReconstructedE, fCentClass);
1264           //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END CHARGED SECONDARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1265           }
1266           else{
1267           //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BEGIN NEUTRAL SECONDARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1268             if(nottrackmatched  && !countasmatched){//secondaries not removed
1269 //            cout<<"Neutral Secondary ";
1270 //            PrintFamilyTreeShort(iPart, stack);
1271
1272               subtotalSecondaryEnergy += fReconstructedEt;
1273               Bool_t tmpIsWritten = kFALSE;
1274               if(fDepositedCode==fgProtonCode ||  fDepositedCode==fgAntiProtonCode || fDepositedCode==fgPiPlusCode || fDepositedCode==fgPiMinusCode || fDepositedCode==fgKPlusCode || fDepositedCode==fgKMinusCode){
1275                 subtotalSecondaryChargedEnergy += fReconstructedEt;
1276                 tmpIsWritten = kTRUE;
1277               }
1278               if(fDepositedCode==fgNeutronCode ||  fDepositedCode==fgAntiNeutronCode){
1279                 subtotalSecondaryNeutronEnergy += fReconstructedEt;
1280                 tmpIsWritten = kTRUE;
1281               }
1282               if(fDepositedCode==fgGammaCode){
1283                 subtotalSecondaryNeutronEnergy += fReconstructedEt;
1284                 tmpIsWritten = kTRUE;
1285               }
1286               if(fDepositedCode==fgEPlusCode || fDepositedCode==fgEMinusCode){
1287                 subtotalSecondaryElectronEnergy += fReconstructedEt;
1288                 tmpIsWritten = kTRUE;
1289               }
1290               if(!tmpIsWritten){
1291                 subtotalSecondaryOtherEnergy += fReconstructedEt;
1292               }
1293
1294               if(!fSelector->CutGeometricalAcceptance(*part)){
1295                 if(TMath::Sqrt(part->Vx() * part->Vx() + part->Vy() * part->Vy() + part->Vz()* part->Vz())>100){
1296                   etSecondariesOutOfAccEffCorr += clEt;
1297                 }
1298               }
1299               else{
1300                 if(TMath::Sqrt(part->Vx() * part->Vx() + part->Vy() * part->Vy() + part->Vz()* part->Vz())>430){
1301                   //clusters which are in the cover of the PHOS/EMCal
1302                   etSecondariesDetectorCoverEffCorr += clEt;
1303                 }
1304               }
1305               fSecondaryClusterEnergy->Fill(fReconstructedEt);
1306               fHistSecondaryPositionInDetector->Fill(part->Vx(),part->Vy(),part->Vz());
1307               fSecondaryNotRemoved++;
1308               etSecondaries += fReconstructedEt;
1309               etSecondariesEffCorr += clEt;
1310               fHistNeutralNotRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1311               nNeutralSecondariesNotRemoved++;
1312             }
1313             else{//secondaries removed
1314                 fNeutralRemoved++;
1315                 fEnergyNeutralRemoved += clEt;
1316                 fHistRemovedOrNot->Fill(1.0, fCentClass);
1317                 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1318                 fHistNeutralRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1319
1320                 nNeutralSecondariesRemoved++;
1321             }
1322           }//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END NEUTRAL SECONDARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1323         }//===================================END SECONDARIES==================================================
1324         else{//not a secondary
1325           //===================================BEGIN CHARGED TRACKS==================================================
1326           if (fDepositedCharge != 0 && fDepositedCode!=fgEMinusCode && fDepositedCode!=fgEPlusCode){//if the particle hitting the calorimeter is pi/k/p/mu
1327 //            cout<<"Hadron ";
1328 //            PrintFamilyTreeShort(iPart, stack);
1329             fHistChargedTrackDepositsAllVsPt->Fill(part->Pt(), fCentClass,fReconstructedEt);
1330             fHistChargedTrackDepositsAllVsPtEffCorr->Fill(part->Pt(), fCentClass,clEt);
1331
1332             written = kTRUE;
1333             fClusterMultChargedTracks++;
1334             etPiKPDeposited += clEt;
1335             if(fDepositedCode==fgProtonCode){
1336               etProtonDeposited += clEt;
1337               etProtonDepositedNoEffCorr += fReconstructedEt;
1338             }
1339             if(fDepositedCode==fgAntiProtonCode){
1340               etAntiProtonDeposited += clEt;
1341               etAntiProtonDepositedNoEffCorr += fReconstructedEt;
1342             }
1343             if(nottrackmatched && !countasmatched){//not removed but should be
1344               subtotalHadronEnergy += fReconstructedEt;
1345               etHadronCrossCheckTrue += clEt;
1346               etPiKPDepositedNotTrackMatched += clEt;
1347               if(fDepositedCode==fgProtonCode){
1348                 etProtonDepositedNotTrackMatched += clEt;
1349                 etProtonDepositedNotTrackMatchedNoEffCorr += fReconstructedEt;
1350               }
1351               if(fDepositedCode==fgAntiProtonCode){
1352                 etAntiProtonDepositedNotTrackMatched += clEt;
1353                 etAntiProtonDepositedNotTrackMatchedNoEffCorr += fReconstructedEt;
1354               }
1355               fHistHadronDepositsAll->Fill(part->Pt());
1356               fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);//denominator in track matching efficiency
1357               //we want this to be all hadrons where we actually reconstructed and would accept the track
1358
1359               for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++){
1360                 AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
1361                 if (!track){
1362                   Printf("ERROR: Could not get track %d", iTrack);
1363                   continue;
1364                 }
1365                 else{
1366                   if((Int_t)fSelector->GetLabel(caloCluster,stack)==track->GetLabel()){//then we found the track from the particle that created this
1367                     fHistHadronDepositsAllvsECent->Fill(fReconstructedE, fCentClass);
1368                   }
1369                 }
1370               }
1371               if(fReconstructedEt>0.5) fHistHadronDepositsAllCent500MeV->Fill(part->Pt(), fCentClass);
1372               fChargedNotRemoved++;
1373               fEnergyChargedNotRemoved += clEt;
1374               fHistRemovedOrNot->Fill(2.0, fCentClass);
1375               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1376               fHistChargedTracksAccepted->Fill(fReconstructedEt);
1377               if(fCalcTrackMatchVsMult){
1378                 if(matchedTrackpt<0.5){//if we could never have matched this because of its pt, how much energy did it deposit?
1379                   fHistChargedTracksAcceptedLowPtCent->Fill(fReconstructedEt, fCentClass);
1380                   fHistChargedTracksAcceptedLowPtCentEffCorr->Fill(clEt, fCentClass);
1381                   if(fReconstructedEt>=0.5) fHistChargedTracksAcceptedLowPtCent500MeV->Fill(fReconstructedEt, fCentClass);
1382                   if(pdg && pdg->PdgCode()!=fgAntiProtonCode){
1383                     fHistChargedTracksAcceptedLowPtCentNoAntiProtons->Fill(fReconstructedEt, fCentClass);
1384                   }
1385                   else{
1386                     fHistChargedTracksAcceptedLowPtCentAntiProtons->Fill(fReconstructedEt, fCentClass);
1387                   }
1388                 }
1389                 fHistChargedTracksAcceptedMult->Fill(fReconstructedEt,fClusterMult);
1390                 if(fClusterMult<25){fHistChargedTracksAcceptedPeripheral->Fill(fReconstructedEt);}
1391               }
1392             }
1393             else{//removed and should have been
1394               //if(countasmatched) cout<<" I was counted as matched even though some of my energy might have been saved."<<endl;
1395               //cout<<" t.m. primary"<<endl;
1396               Int_t trackindex =  fSelector->GetLabel(caloCluster,stack);// (caloCluster->GetLabelsArray())->At(0);
1397               fHistChargedTrackDepositsAcceptedVsPt->Fill(part->Pt(), fCentClass,fReconstructedEt);
1398               fHistChargedTrackDepositsAcceptedVsPtEffCorr->Fill(part->Pt(), fCentClass,clEt);
1399               fHistHadronDepositsReco->Fill(part->Pt());
1400               fHistHadronDepositsRecoCent->Fill(part->Pt(), fCentClass);
1401               fHistHadronDepositsRecovsECent->Fill(fReconstructedE, fCentClass);
1402               fHistHadronDepositsAll->Fill(part->Pt());
1403               fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
1404               for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++){
1405                 AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
1406                 if (!track){
1407                   Printf("ERROR: Could not get track %d", iTrack);
1408                   continue;
1409                 }
1410                 else{
1411                   if((Int_t) fSelector->GetLabel(caloCluster,stack)==track->GetLabel()){//then we found the track from the particle that created this
1412                     fHistHadronDepositsAllvsECent->Fill(fReconstructedE, fCentClass);
1413                   }
1414                 }
1415               }
1416               if(fReconstructedEt>0.5) fHistHadronDepositsAllCent500MeV->Fill(part->Pt(), fCentClass);
1417               if((Int_t)fSelector->GetLabel(caloCluster,stack)!= trackindex){
1418                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1419                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1420                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1421               }
1422               else{
1423                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1424               }
1425               fChargedRemoved++;
1426               fEnergyChargedRemoved += clEt;
1427               fHistRemovedOrNot->Fill(0.0, fCentClass);
1428               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1429               fHistChargedTracksCut->Fill(fReconstructedEt);
1430               if(fCalcTrackMatchVsMult){
1431                 fHistChargedTracksCutMult->Fill(fReconstructedEt,fClusterMult);
1432                 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fReconstructedEt);}
1433               }
1434               fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
1435               if(fCalcTrackMatchVsMult){
1436                 if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1437                 fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1438                 fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//fill with the efficiency corrected energy
1439               }
1440             }
1441           }
1442           //===================================END CHARGED TRACKS==================================================
1443           //===================================BEGIN KAONS==================================================
1444           //K0L and any neutral particles from the decay of K+/- or K0S
1445           if(!written && (fPrimaryCode==fgKPlusCode || fPrimaryCode==fgKMinusCode || fPrimaryCode==fgK0SCode ||fPrimaryCode==fgK0LCode) && nottrackmatched){
1446 //          cout<<"Kaon ";
1447 //          PrintFamilyTreeShort(iPart, stack);
1448
1449             etKaonCrossCheckTrue += clEt;
1450             subtotalKaonEnergy += fReconstructedEt;
1451             //etKaonCrossCheckTrue += fReconstructedEt;
1452             written = kTRUE;//At this point we are not tracking them but we don't count them as neutrals accidentally removed
1453           }
1454           //===================================END KAONS==================================================
1455           //===================================BEGIN SIGNAL==================================================
1456           if(!written && (fDepositedCode==fgGammaCode || fDepositedCode==fgEMinusCode || fDepositedCode ==fgEPlusCode)){//if the particle hitting the calorimeter is gamma, electron and not from a kaon
1457 //          cout<<"Signal ";
1458 //          PrintFamilyTreeShort(iPart, stack);
1459 //efficiencies without cutting track matched gammas
1460             if(containsGamma){
1461               if(!fSelector->CutGeometricalAcceptance(*primPart)){
1462                 fHistAllGammasFoundOutOfAccRecoEnergyCent->Fill(fReconstructedE,fCentClass);
1463               }
1464               else{
1465                 fHistAllGammasFoundRecoEnergyCent->Fill(fReconstructedE,fCentClass);
1466               }
1467               //cout<<"contains gamma"<<endl;
1468             }
1469             //else{cout<<"does not contain gamma"<<endl;}
1470
1471
1472             fClusterMultGammas++;
1473             written = kTRUE;
1474             if(nottrackmatched){//Not removed and not supposed to be removed - signal
1475               subtotalSignalEnergy += fReconstructedEt;
1476               subtotalReconstructedSignalEnergy += fReconstructedEt;
1477               if(fPrimaryCode==fgGammaCode){//if this comes from a primary electron - this is for reconstruction efficiency
1478                 subtotalReconstructedGammaEnergy += fReconstructedEt;
1479               }
1480 //            else{
1481 //              if(fPrimaryCode!=fgEMinusCode && fPrimaryCode !=fgEPlusCode){
1482 //                PrintFamilyTreeShort(iPart, stack);
1483 //              }
1484 //            }
1485               fEtNonRemovedGammas += clEt;
1486               fMultNonRemovedGammas++;
1487               fNeutralNotRemoved--;
1488               fEnergyNeutralNotRemoved -= clEt;
1489               fHistGammasAccepted->Fill(fDepositedEt);
1490               if(fCalcTrackMatchVsMult){
1491                 fHistGammasAcceptedMult->Fill(fDepositedEt,fClusterMult);
1492                 if(fClusterMult<25){fHistGammasAcceptedPeripheral->Fill(fDepositedEt);}
1493               }
1494               if(containsGamma){
1495                 if(!fSelector->CutGeometricalAcceptance(*primPart)){
1496                   fHistGammasFoundOutOfAccRecoEnergyCent->Fill(fReconstructedE,fCentClass);
1497                 }
1498                 else{
1499                   //cout<<"filling reference histogram"<<endl;
1500                   if(gammaEnergy>0) fHistGammaEnergyCrossCheckAlt->Fill(gammaEnergy,(fReconstructedEt-gammaEnergy)/fReconstructedEt);
1501                   fHistGammasFoundRecoEnergyCent->Fill(fReconstructedE,fCentClass);
1502                   etGammaCrossCheckAlt += clEt;
1503                 }
1504               }
1505               if(fDepositedCode==fgGammaCode){
1506                 if(!fSelector->CutGeometricalAcceptance(*primPart)){
1507                   fHistGammasFoundOutOfAccCent->Fill(fReconstructedE,fCentClass);
1508                 }
1509                 else{
1510                   etGammaCrossCheck += clEt;
1511                   // cout<<"Found gamma et "<<fReconstructedEt<<" sim et "<< fPrimaryEt<<endl;
1512                   if(fPrimaryEt>0) fHistGammaEnergyCrossCheck->Fill(fPrimaryEt,(fReconstructedEt-fPrimaryEt)/fReconstructedEt);
1513                   fHistGammasFoundAltCent->Fill(fReconstructedE,fCentClass);
1514                 }
1515               }
1516             }
1517             else{//removed but shouldn't have been
1518               Int_t trackindex = fSelector->GetLabel(caloCluster,stack);// (caloCluster->GetLabelsArray())->At(1);
1519               if(caloCluster->GetLabel()!=trackindex){
1520                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1521                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1522 //              cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1523 //              PrintFamilyTreeShort(trackindex, stack);
1524 //              cout<<"Cluster"<<endl;
1525               }
1526               fGammaRemoved++;
1527               fGammaRemovedEt+=clEt;
1528               fHistGammasCut->Fill(fDepositedEt);
1529               if(fCalcTrackMatchVsMult){
1530                 fHistGammasCutMult->Fill(fDepositedEt,fClusterMult);
1531                 if(fClusterMult<25){fHistGammasCutPeripheral->Fill(fDepositedEt);}
1532               }
1533               fHistMatchedTracksEvspTSignal->Fill(matchedTrackp,fReconstructedE);
1534               if(fCalcTrackMatchVsMult){
1535                 if(fClusterMult<25){fHistMatchedTracksEvspTSignalPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1536                 fHistMatchedTracksEvspTSignalvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1537                 fHistMatchedTracksEvspTSignalvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);
1538               }
1539             }
1540           }
1541           //===================================END SIGNAL==================================================
1542           //===================================BEGIN NEUTRONS==================================================
1543           //all other cases - neutron, anti-neutron, not aware of other cases
1544           if(!written){
1545 //          cout<<"Neutron ";
1546 //          PrintFamilyTreeShort(iPart, stack);
1547             fNeutralNotRemoved++;
1548             fEnergyNeutralNotRemoved += clEt;//this is the efficiency corrected energy
1549             fHistRemovedOrNot->Fill(3.0, fCentClass);
1550             if ((fDepositedCode == fgNeutronCode || fDepositedCode == fgAntiNeutronCode) && nottrackmatched){
1551               etNeutronCrossCheckTrue += clEt;
1552               subtotalNeutronEnergy += fReconstructedEt;
1553               fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1554               fHistNeutronsEtVsCent->Fill(clEt,fCentClass);
1555               nNeutrons++;
1556               if(fDepositedCode == fgNeutronCode){
1557                 etNeutronDeposited += clEt;
1558                 etNeutronDeposited += fReconstructedEt;
1559               }
1560               if(fDepositedCode == fgAntiNeutronCode){
1561                 etAntiNeutronDeposited += clEt;
1562                 etAntiNeutronDeposited += fReconstructedEt;
1563               }
1564               //cout<<"I am a";
1565               //if(fDepositedCode == fgAntiNeutronCode) cout<<"n anti-";
1566               //else{cout<<" ";}
1567               //cout<<"neutron!! pt "<<part->Pt()<<" eta "<<part->Eta()<<" phi "<<part->Phi()<<" Et "<<clEt<<endl;
1568               //PrintFamilyTreeShort(iPart, stack);
1569             }
1570             else{
1571               nNotNeutrons++;
1572 //            cout<<"I am a";
1573 //            if(fDepositedCode == fgAntiNeutronCode) cout<<"n anti-";
1574 //            else{cout<<" ";}
1575 //            cout<<"neutron!! pt "<<part->Pt()<<" eta "<<part->Eta()<<" phi "<<part->Phi()<<" Et "<<clEt<<endl;
1576 //            PrintFamilyTreeShort(iPart, stack);
1577 //            cout<<"I am not a neutron!!"<<endl;
1578 //            PrintFamilyTreeShort(iPart, stack);
1579             }
1580           }
1581           //===================================END NEUTRONS==================================================
1582         }
1583         if(fCuts->GetHistMakeTree()) fPrimaryTree->Fill();
1584     } // end of loop over clusters
1585
1586
1587     fHistAllEnergy->Fill(fCentClass,subtotalAllEnergy);
1588     fHistSignalEnergy->Fill(fCentClass,subtotalSignalEnergy);
1589     fHistNeutronEnergy->Fill(fCentClass,subtotalNeutronEnergy);
1590     fHistKaonEnergy->Fill(fCentClass,subtotalKaonEnergy);
1591     fHistHadronEnergy->Fill(fCentClass,subtotalHadronEnergy);
1592     fHistSecondaryEnergy->Fill(fCentClass,subtotalSecondaryEnergy);
1593     fHistSecondaryChargedEnergy->Fill(fCentClass,subtotalSecondaryChargedEnergy);
1594     fHistSecondaryNeutronEnergy->Fill(fCentClass,subtotalSecondaryNeutronEnergy);
1595     fHistSecondaryGammaEnergy->Fill(fCentClass,subtotalSecondaryGammaEnergy);
1596     fHistSecondaryElectronEnergy->Fill(fCentClass,subtotalSecondaryElectronEnergy);
1597     fHistSecondaryOtherEnergy->Fill(fCentClass,subtotalSecondaryOtherEnergy);
1598     fHistReconstructedGammaEnergy->Fill(fCentClass,subtotalReconstructedGammaEnergy);
1599     fHistReconstructedSignalEnergy->Fill(fCentClass,subtotalReconstructedSignalEnergy);
1600     Float_t fracSignal =0;
1601     Float_t fracHadron = 0;
1602     Float_t fracKaon = 0;
1603     Float_t fracNeutron = 0;
1604     Float_t fracSecondary = 0;
1605     if(subtotalAllEnergy>0){
1606       fracSignal = subtotalReconstructedSignalEnergy/subtotalAllEnergy;
1607       fracHadron = subtotalHadronEnergy/subtotalAllEnergy;
1608       fracKaon = subtotalKaonEnergy/subtotalAllEnergy;
1609       fracNeutron = subtotalNeutronEnergy/subtotalAllEnergy;
1610       fracSecondary = subtotalSecondaryEnergy/subtotalAllEnergy;
1611     }
1612      fHistFracSignalVsNClusters->Fill(fClusterMult,fracSignal);
1613      fHistFracHadronsVsNClusters->Fill(fClusterMult,fracHadron);
1614      fHistFracKaonsVsNClusters->Fill(fClusterMult,fracKaon);
1615      fHistFracNeutronsVsNClusters->Fill(fClusterMult,fracNeutron);
1616      fHistFracSecondariesVsNClusters->Fill(fClusterMult,fracSecondary);
1617      fHistFracSignalVsNMultiplicity->Fill(multiplicity,fracSignal);
1618     fHistFracHadronsVsNMultiplicity->Fill(multiplicity,fracHadron);
1619     fHistFracKaonsVsNMultiplicity->Fill(multiplicity,fracKaon);
1620     fHistFracNeutronsVsNMultiplicity->Fill(multiplicity,fracNeutron);
1621     fHistFracSecondariesVsNMultiplicity->Fill(multiplicity,fracSecondary);
1622      fHistFracSignalVsNMatchedTracks->Fill(nChargedHadronsMeasured,fracSignal);
1623     fHistFracHadronsVsNMatchedTracks->Fill(nChargedHadronsMeasured,fracHadron);
1624     fHistFracKaonsVsNMatchedTracks->Fill(nChargedHadronsMeasured,fracKaon);
1625     fHistFracNeutronsVsNMatchedTracks->Fill(nChargedHadronsMeasured,fracNeutron);
1626     fHistFracSecondariesVsNMatchedTracks->Fill(nChargedHadronsMeasured,fracSecondary);
1627      fHistFracSignalVsNTotalTracks->Fill(nChargedHadronsTotal,fracSignal);
1628     fHistFracHadronsVsNTotalTracks->Fill(nChargedHadronsTotal,fracHadron);
1629     fHistFracKaonsVsNTotalTracks->Fill(nChargedHadronsTotal,fracKaon);
1630     fHistFracNeutronsVsNTotalTracks->Fill(nChargedHadronsTotal,fracNeutron);
1631     fHistFracSecondariesVsNTotalTracks->Fill(nChargedHadronsTotal,fracSecondary);
1632
1633
1634
1635     fHistPiKPNotTrackMatchedDepositedVsNch->Fill(etPiKPDepositedNotTrackMatched,multiplicity);
1636     fHistPiKPDepositedVsNch->Fill(etPiKPDeposited,multiplicity);
1637
1638     fHistNeutronsDepositedVsNch->Fill(etNeutronDeposited,multiplicity);
1639     fHistAntiNeutronsDepositedVsNch->Fill(etAntiNeutronDeposited,multiplicity);
1640     fHistProtonsDepositedVsNch->Fill(etProtonDeposited,multiplicity);
1641     fHistAntiProtonsDepositedVsNch->Fill(etAntiProtonDeposited,multiplicity);
1642     fHistProtonsNotTrackMatchedDepositedVsNch->Fill(etProtonDepositedNotTrackMatched,multiplicity);
1643     fHistAntiProtonsNotTrackMatchedDepositedVsNch->Fill(etAntiProtonDepositedNotTrackMatched,multiplicity);
1644
1645     fHistNeutronsDepositedVsNcl->Fill(etNeutronDeposited,fClusterMult);
1646     fHistAntiNeutronsDepositedVsNcl->Fill(etAntiNeutronDeposited,fClusterMult);
1647     fHistProtonsDepositedVsNcl->Fill(etProtonDeposited,fClusterMult);
1648     fHistAntiProtonsDepositedVsNcl->Fill(etAntiProtonDeposited,fClusterMult);
1649     fHistProtonsNotTrackMatchedDepositedVsNcl->Fill(etProtonDepositedNotTrackMatched,fClusterMult);
1650     fHistAntiProtonsNotTrackMatchedDepositedVsNcl->Fill(etAntiProtonDepositedNotTrackMatched,fClusterMult);
1651
1652     fHistSecondariesVsNch->Fill(etSecondaries,multiplicity);
1653     fHistSecondariesVsNcl->Fill(etSecondaries,fClusterMult);
1654     fHistSecondariesEffCorrVsNch->Fill(etSecondariesEffCorr,multiplicity);
1655     fHistSecondariesEffCorrVsNcl->Fill(etSecondariesEffCorr,fClusterMult);
1656     fHistSecondariesOutOfAccEffCorrVsNch->Fill(etSecondariesOutOfAccEffCorr,multiplicity);
1657     fHistSecondariesDetectorCoverEffCorrVsNch->Fill(etSecondariesDetectorCoverEffCorr,multiplicity);
1658     fHistCentVsNchVsNcl->Fill(fCentClass,multiplicity, fClusterMult);
1659
1660     fHistNeutronsDepositedVsNchNoEffCorr->Fill(etNeutronDepositedNoEffCorr,multiplicity);
1661     fHistAntiNeutronsDepositedVsNchNoEffCorr->Fill(etAntiNeutronDepositedNoEffCorr,multiplicity);
1662     fHistNeutronsDepositedVsNclNoEffCorr->Fill(etNeutronDepositedNoEffCorr,fClusterMult);
1663     fHistAntiNeutronsDepositedVsNclNoEffCorr->Fill(etAntiNeutronDepositedNoEffCorr,fClusterMult);
1664
1665     fHistProtonsDepositedVsNchNoEffCorr->Fill(etProtonDepositedNoEffCorr,multiplicity);
1666     fHistAntiProtonsDepositedVsNchNoEffCorr->Fill(etAntiProtonDepositedNoEffCorr,multiplicity);
1667     fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr->Fill(etProtonDepositedNotTrackMatchedNoEffCorr,multiplicity);
1668     fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr->Fill(etAntiProtonDepositedNotTrackMatchedNoEffCorr,multiplicity);
1669     fHistProtonsDepositedVsNclNoEffCorr->Fill(etProtonDepositedNoEffCorr,fClusterMult);
1670     fHistAntiProtonsDepositedVsNclNoEffCorr->Fill(etAntiProtonDepositedNoEffCorr,fClusterMult);
1671     fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr->Fill(etProtonDepositedNotTrackMatchedNoEffCorr,fClusterMult);
1672     fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr->Fill(etAntiProtonDepositedNotTrackMatchedNoEffCorr,fClusterMult);
1673
1674
1675     fHistNeutronsNumVsCent->Fill(nNeutrons,fCentClass);
1676     fHistNotNeutronsNumVsCent->Fill(nNotNeutrons,fCentClass);
1677
1678     std::sort(foundGammas.begin(), foundGammas.end());
1679     for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
1680     {
1681
1682         if(!stack->IsPhysicalPrimary(iPart)) continue;
1683         //Some productions have signals added in.  This switch allows the explicit exclusion of these added in signals when running over the data.
1684         if(checkLabelForHIJING && !IsHIJINGLabel(iPart,mcEvent,stack) ) continue;
1685
1686         TParticle *part = stack->Particle(iPart);
1687
1688         if (!part)
1689         {
1690           //Printf("ERROR: Could not get particle %d", iPart);
1691             continue;
1692         }
1693         TParticlePDG *pdg = part->GetPDG(0);
1694
1695         if (!pdg)
1696         {
1697           //Printf("ERROR: Could not get particle PDG %d", iPart);
1698             continue;
1699         }
1700
1701         if(pdg->PdgCode()==fgGammaCode)// TMath::Abs(part->Eta()) < 0.12)
1702         {
1703           if(fSelector->CutGeometricalAcceptance(*part)){
1704             subtotalSimulatedGammaEnergy += part->Energy()*TMath::Sin(part->Theta());
1705             if(fSelector->PassMinEnergyCut(*part)){
1706               //cout<<"Sim gamma et "<< part->Energy()*TMath::Sin(part->Theta())<<endl;
1707               etGammaCrossCheckTrue += part->Energy()*TMath::Sin(part->Theta());
1708               subtotalSimulatedGammaEnergyAboveThreshold += part->Energy()*TMath::Sin(part->Theta());
1709             }
1710             fHistGammasGenerated->Fill(part->Energy());
1711             fHistGammasGeneratedCent->Fill(part->Energy(),fCentClass);
1712           }
1713           if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
1714           {
1715             if(!fSelector->CutGeometricalAcceptance(*part)){
1716               //cout<<"Gamma NOT in acceptance"<<endl;
1717               fHistGammasFoundOutOfAccCent->Fill(part->Energy(),fCentClass);
1718             }
1719             else{
1720               fHistGammasFound->Fill(part->Energy());
1721               fHistGammasFoundCent->Fill(part->Energy(),fCentClass);
1722               //cout<<"Gamma IN acceptance"<<endl;
1723             }
1724           }
1725         }
1726         if(pdg->PdgCode()==fgPiPlusCode || pdg->PdgCode()==fgPiMinusCode || pdg->PdgCode()==fgProtonCode || pdg->PdgCode()==fgAntiProtonCode){//section here for all hadrons generated
1727           //if(fSelector->CutGeometricalAcceptance(*part)){
1728             fHistHadronsAllCent->Fill(part->Pt(), fCentClass);
1729             //}
1730         }
1731
1732
1733     }
1734     if(fCalcForKaonCorrection){
1735       Float_t etCuts[11] = {0.0,0.05,0.1,0.15,0.2,0.25, 0.3,0.35,0.4,0.45,0.5};
1736       Int_t nEtCuts = 11;
1737       //loop over simulated particles in order to find K0S
1738       for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++){
1739         //Some productions have signals added in.  This switch allows the explicit exclusion of these added in signals when running over the data.
1740         if(checkLabelForHIJING && !IsHIJINGLabel(iPart,mcEvent,stack) ) continue;
1741         TParticle *part = stack->Particle(iPart);
1742         if (!part){
1743           //Printf("ERROR: Could not get particle %d", iPart);
1744           continue;
1745         }
1746         TParticlePDG *pdg = part->GetPDG(0);
1747         if (!pdg){
1748           //Printf("ERROR: Could not get particle PDG %d", iPart);
1749           continue;
1750         }
1751         //if(stack->IsPhysicalPrimary(iPart)){//if it is a K0 it might have decayed into four pions
1752         //fgK0SCode, fgGammaCode, fgPi0Code
1753         Int_t code = pdg->PdgCode();
1754         if(code == fgK0SCode || code==fgKPlusCode || code==fgKMinusCode ||code==fgK0LCode || code==fgK0Code){//this is a kaon
1755           //cout<<"I am a kaon too! "<<stack->Particle(iPart)->GetName()<<" "<<code<<endl;
1756           Float_t pTk = stack->Particle(iPart)->Pt();
1757           if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//these are particles which would be included in our spectra measurements
1758             fHistSimKaonsInAcceptance->Fill(pTk);
1759             if(code == fgK0SCode){fHistSimK0SInAcceptance->Fill(pTk);}
1760             if(code == fgK0LCode){fHistSimK0LInAcceptance->Fill(pTk);}
1761             if(code == fgKPlusCode){fHistSimKPlusInAcceptance->Fill(pTk);}
1762             if(code == fgKMinusCode){fHistSimKMinusInAcceptance->Fill(pTk);}
1763             if(code == fgK0Code){//Split K0's between the two
1764               fHistSimK0SInAcceptance->Fill(pTk,0.5);
1765               fHistSimK0LInAcceptance->Fill(pTk,0.5);
1766             }
1767           }
1768           else{
1769             fHistSimKaonsOutOfAcceptance->Fill(pTk);
1770             //      if(!stack->IsPhysicalPrimary(iPart)){
1771             //        PrintFamilyTree(iPart, stack);
1772             //      }
1773           }
1774           Float_t totalGammaEts[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,  0.0};
1775           Float_t totalClusterEts[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,  0.0};
1776           for (int iCluster = 0; iCluster < nCluster; iCluster++ ){//if this cluster is from any of the decay daughters of any kaon...  but there is no easy way to look at this so we loop over clusters...
1777             AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
1778             if (!fSelector->CutGeometricalAcceptance(*caloCluster)) continue;
1779             const Int_t myPart = fSelector->GetLabel(caloCluster,stack);
1780             //const Int_t myPart = TMath::Abs(caloCluster->GetLabel());
1781             //identify the primary particle which created this cluster
1782             int primIdx = myPart;
1783             if (!stack->IsPhysicalPrimary(myPart)){
1784               primIdx = GetPrimMother(iPart, stack);
1785             } // end of primary particle check
1786             TParticle *hitPart = stack->Particle(myPart);
1787             Bool_t hitsAsChargedKaon = kFALSE;
1788             if(hitPart->GetPdgCode()== fgKPlusCode || hitPart->GetPdgCode()== fgKPlusCode){
1789               if(myPart==primIdx){
1790                 //The particle hits as a charged kaon and that kaon is a primary kaon - do not count because this is counted in the hadronic correction!
1791                 hitsAsChargedKaon = kTRUE;
1792                 //cout<<"Found primary charged kaon cluster!"<<endl;
1793               }
1794             }
1795             if(primIdx==iPart && primIdx>0 && !hitsAsChargedKaon){//This cluster is from our primary particle and our primary particle is a kaon
1796               //cout<<"I have a particle match! prim code"<<code<<" id "<<primIdx <<endl;
1797               Float_t pos[3];
1798               caloCluster->GetPosition(pos);
1799               TVector3 cp(pos);
1800               Double_t clEt = caloCluster->E()*TMath::Sin(cp.Theta());
1801               Double_t clEtCorr = CorrectForReconstructionEfficiency(*caloCluster,fCentClass);
1802               for(int l=0;l<nEtCuts;l++){//loop over cut values
1803                 if(clEt>=etCuts[l]){
1804                   //cout<<", "<<clEt<<">="<<etCuts[l];
1805                   totalClusterEts[l] += clEtCorr;//if cluster et is above the cut off energy add it
1806                   totalGammaEts[l] += clEt;//if cluster et is above the cut off energy add it
1807                 }
1808               }
1809             }
1810           }
1811           //      cout<<"Deposits:  pT: "<<pTk;
1812           //      for(int l=0;l<nEtCuts;l++){//loop over cut values
1813           //        cout<<" "<<totalClusterEts[l];
1814           //      }
1815           //      cout<<endl;
1816           if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//within the acceptance of our spectra and is a primary particle
1817             if(totalClusterEts[0]>0.0){fHistSimKaonsInAcceptanceWithDepositsPrimaries->Fill(pTk);}
1818             //cout<<"I have a particle match! prim code"<<code<<" id "<<iPart <<endl;
1819             for(int l=0;l<nEtCuts;l++){
1820               fHistK0EDepositsVsPtInAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1821               fHistK0EGammaVsPtInAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
1822             }
1823           }
1824           else{//outside the acceptance of our spectra
1825             if(totalClusterEts[0]>0.0){
1826               if(stack->IsPhysicalPrimary(iPart)){fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries->Fill(pTk);}
1827               else{fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries->Fill(pTk);}
1828             }
1829             for(int l=0;l<nEtCuts;l++){
1830               fHistK0EDepositsVsPtOutOfAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1831               fHistK0EGammaVsPtOutOfAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
1832             }
1833           }
1834
1835         }
1836       }
1837     }
1838
1839     fHistSimulatedGammaEnergy->Fill(fCentClass,subtotalSimulatedGammaEnergy);
1840     fHistSimulatedGammaEnergyAboveThreshold->Fill(fCentClass,subtotalSimulatedGammaEnergyAboveThreshold);
1841
1842     fHistMultChVsSignalVsMult->Fill(fClusterMultChargedTracks,fClusterMultGammas,fClusterMult);
1843     fHistNeutralRemovedSecondaryNumVsNCluster->Fill(nNeutralSecondariesRemoved,fClusterMult);
1844     fHistChargedRemovedSecondaryNumVsNCluster->Fill(nChargedSecondariesRemoved,fClusterMult);
1845     fHistNeutralNotRemovedSecondaryNumVsNCluster->Fill(nNeutralSecondariesNotRemoved,fClusterMult);
1846     fHistChargedNotRemovedSecondaryNumVsNCluster->Fill(nChargedSecondariesNotRemoved,fClusterMult);
1847     fHistNeutralRemovedSecondaryNumVsCent->Fill(nNeutralSecondariesRemoved,fCentClass);
1848     fHistChargedRemovedSecondaryNumVsCent->Fill(nChargedSecondariesRemoved,fCentClass);
1849     fHistNeutralNotRemovedSecondaryNumVsCent->Fill(nNeutralSecondariesNotRemoved,fCentClass);
1850     fHistChargedNotRemovedSecondaryNumVsCent->Fill(nChargedSecondariesNotRemoved,fCentClass);
1851     if(etGammaCrossCheckTrue>0)fHistGammaCrossCheck->Fill(fCentClass,(etGammaCrossCheck-etGammaCrossCheckTrue)/etGammaCrossCheckTrue);
1852     if(etGammaCrossCheckTrue>0)fHistGammaCrossCheckAlt->Fill(fCentClass,(etGammaCrossCheckAlt-etGammaCrossCheckTrue)/etGammaCrossCheckTrue);
1853     if(fTmCorrections->GetNeutronCorrection(fCentClass)>0)fHistNeutronCrossCheck->Fill(fCentClass,(fTmCorrections->GetNeutronCorrection(fCentClass)-etNeutronCrossCheckTrue)/fTmCorrections->GetNeutronCorrection(fCentClass));
1854     if(fTmCorrections->GetSecondaryCorrection(fCentClass)>0)fHistSecondaryCrossCheck->Fill(fCentClass,(fTmCorrections->GetSecondaryCorrection(fCentClass)-etSecondaryCrossCheckTrue)/fTmCorrections->GetSecondaryCorrection(fCentClass));
1855     if(fTmCorrections->GetHadronCorrection(fCentClass)>0)fHistHadronCrossCheck->Fill(fCentClass,(fTmCorrections->GetHadronCorrection(fCentClass)-etHadronCrossCheckTrue)/fTmCorrections->GetHadronCorrection(fCentClass));
1856     if(fTmCorrections->GetKaonCorrection(fCentClass)>0)fHistKaonCrossCheck->Fill(fCentClass,(fTmCorrections->GetKaonCorrection(fCentClass)-etKaonCrossCheckTrue)/fTmCorrections->GetKaonCorrection(fCentClass));
1857     //cout<<"Secondaries not removed: "<<fSecondaryNotRemoved<<" neutral secondaries not removed "<<nNeutralSecondariesNotRemoved<<" cluster mult "<<fClusterMult<<endl;
1858     FillHistograms();
1859     return 0;
1860 }
1861
1862 void AliAnalysisEtMonteCarlo::Init()
1863 {   // init
1864     AliAnalysisEt::Init();
1865 }
1866
1867 void AliAnalysisEtMonteCarlo::ResetEventValues()
1868 {   // reset event values
1869     AliAnalysisEt::ResetEventValues();
1870   if(!fIsMC) return;
1871
1872     fTotEtSecondary = 0;
1873     fTotEtSecondaryFromEmEtPrimary = 0;
1874     fTotEtWithSecondaryRemoved = 0;
1875
1876     // collision geometry defaults for p+p:
1877     fImpactParameter = 0;
1878     fNcoll = 1;
1879     fNpart = 2;
1880
1881     fEtNonRemovedProtons = 0;
1882     fEtNonRemovedAntiProtons = 0;
1883     fEtNonRemovedPiPlus = 0;
1884     fEtNonRemovedPiMinus = 0;
1885     fEtNonRemovedKaonPlus = 0;
1886     fEtNonRemovedKaonMinus = 0;
1887     fEtNonRemovedK0S = 0;
1888     fEtNonRemovedK0L = 0;
1889     fEtNonRemovedLambdas = 0;
1890     fEtNonRemovedElectrons = 0;
1891     fEtNonRemovedPositrons = 0;
1892     fEtNonRemovedMuPlus = 0;
1893     fEtNonRemovedMuMinus = 0;
1894     fEtNonRemovedNeutrons = 0;
1895     fEtNonRemovedAntiNeutrons = 0;
1896     fEtNonRemovedGammas = 0;
1897     fEtNonRemovedGammasFromPi0 = 0;
1898
1899     fEtRemovedProtons = 0;
1900     fEtRemovedAntiProtons = 0;
1901     fEtRemovedPiPlus = 0;
1902     fEtRemovedPiMinus = 0;
1903     fEtRemovedKaonPlus = 0;
1904     fEtRemovedKaonMinus = 0;
1905     fEtRemovedK0s = 0;
1906     fEtRemovedK0L = 0;
1907     fEtRemovedLambdas = 0;
1908     fEtRemovedElectrons = 0;
1909     fEtRemovedPositrons = 0;
1910     fEtRemovedMuPlus = 0;
1911     fEtRemovedMuMinus = 0;
1912     fEtRemovedNeutrons = 0;
1913
1914     fEtRemovedGammasFromPi0 = 0;
1915     fEtRemovedGammas = 0;
1916     fEtRemovedNeutrons = 0;
1917     fEtRemovedAntiNeutrons = 0;
1918
1919     fMultNonRemovedProtons = 0;
1920     fMultNonRemovedAntiProtons = 0;
1921     fMultNonRemovedPiPlus = 0;
1922     fMultNonRemovedPiMinus = 0;
1923     fMultNonRemovedKaonPlus = 0;
1924     fMultNonRemovedKaonMinus = 0;
1925     fMultNonRemovedK0s = 0;
1926     fMultNonRemovedK0L = 0;
1927     fMultNonRemovedLambdas = 0;
1928     fMultNonRemovedElectrons = 0;
1929     fMultNonRemovedPositrons = 0;
1930     fMultNonRemovedMuPlus = 0;
1931     fMultNonRemovedMuMinus = 0;
1932     fMultNonRemovedNeutrons = 0;
1933     fMultNonRemovedAntiNeutrons = 0;
1934     fMultNonRemovedGammas = 0;
1935
1936     fMultRemovedProtons = 0;
1937     fMultRemovedAntiProtons = 0;
1938     fMultRemovedPiPlus = 0;
1939     fMultRemovedPiMinus = 0;
1940     fMultRemovedKaonPlus = 0;
1941     fMultRemovedKaonMinus = 0;
1942     fMultRemovedK0s = 0;
1943     fMultRemovedK0L = 0;
1944     fMultRemovedLambdas = 0;
1945     fMultRemovedElectrons = 0;
1946     fMultRemovedPositrons = 0;
1947     fMultRemovedMuPlus = 0;
1948     fMultRemovedMuMinus = 0;
1949
1950     fMultRemovedGammas = 0;
1951     fMultRemovedNeutrons = 0;
1952     fMultRemovedAntiNeutrons = 0;
1953
1954     fEnergyChargedNotRemoved = 0;
1955     fEnergyChargedRemoved = 0;
1956     fEnergyNeutralNotRemoved = 0;
1957     fEnergyNeutralRemoved = 0;
1958
1959     fChargedNotRemoved = 0;
1960     fChargedRemoved = 0;
1961     fNeutralNotRemoved = 0;
1962     fNeutralRemoved = 0;
1963     fGammaRemoved = 0;
1964     fSecondaryNotRemoved = 0;
1965
1966     fTrackMultInAcc = 0;
1967
1968     fTotNeutralEtAfterMinEnergyCut = 0;
1969
1970     fSecondaryNotRemoved = 0;
1971
1972     fTotPx = 0;
1973     fTotPy = 0;
1974     fTotPz = 0;
1975
1976
1977 }
1978
1979 void AliAnalysisEtMonteCarlo::CreateHistograms()
1980 {   // histogram related additions
1981     AliAnalysisEt::CreateHistograms();
1982     if(!fIsMC) return;
1983     if (fEventSummaryTree) {
1984         fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1985         fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1986         fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1987         fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1988         fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1989         fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1990         fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
1991         fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1992         fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1993         fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1994         fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1995         fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1996         fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1997         fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1998         fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1999         fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
2000 //      fEventSummaryTree->Branch("f
2001     }
2002     Float_t scale = 1;//scale up histograms if EMCal 2011 so we have the right range
2003     if(fDataSet==2011   && !fHistogramNameSuffix.Contains("P")){
2004       scale = 2.5;
2005     }
2006
2007
2008     //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
2009     //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
2010     //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
2011     //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
2012
2013     fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
2014
2015     //I don't think any of these are actually used...  flag for later removal
2016     Int_t nbins = 15;
2017     fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", nbins, 0, 30, 10, -0.5, 9.5);
2018     fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", nbins, 0, 30, 10, -0.5, 9.5);
2019     fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", nbins, 0, 30, 10, -0.5, 9.5);
2020     fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", nbins, 0, 30, 10, -0.5, 9.5);
2021     fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", nbins, 0, 30, 10, -0.5, 9.5);
2022     fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", nbins, 0, 30, 10, -0.5, 9.5);
2023     fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", nbins, 0, 30, 10, -0.5, 9.5);
2024     fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", nbins, 0, 30, 10, -0.5, 9.5);
2025     fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", nbins, 0, 30, 10, -0.5, 9.5);
2026     fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", nbins, 0, 30, 10, -0.5, 9.5);
2027     fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", nbins, 0, 30, 10, -0.5, 9.5);
2028     fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", nbins, 0, 30, 10, -0.5, 9.5);
2029     fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", nbins, 0, 30, 10, -0.5, 9.5);
2030     fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", nbins, 0, 30, 10, -0.5, 9.5);
2031     fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", nbins, 0, 30, 10, -0.5, 9.5);
2032
2033     fHistEtNonRemovedGammas = new  TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", nbins, 0, 30, 10, -0.5, 9.5);
2034     fHistEtNonRemovedGammasFromPi0 = new  TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", nbins, 0, 30, 10, -0.5, 9.5);
2035
2036     fHistEtRemovedGammas = new  TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", nbins, 0, 30, 10, -0.5, 9.5);
2037     fHistEtRemovedNeutrons = new  TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", nbins, 0, 30, 10, -0.5, 9.5);
2038     fHistEtRemovedAntiNeutrons = new  TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", nbins, 0, 30, 10, -0.5, 9.5);
2039
2040     fHistEtRemovedCharged = new  TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", nbins, 0, 30, 10, -0.5, 9.5);
2041     fHistEtRemovedNeutrals = new  TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", nbins, 0, 30, 10, -0.5, 9.5);
2042
2043     fHistEtNonRemovedCharged = new  TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", nbins, 0, 30, 10, -0.5, 9.5);
2044     fHistEtNonRemovedNeutrals = new  TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", nbins, 0, 30, 10, -0.5, 9.5);
2045
2046     Int_t nmultbins = 10;
2047     fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2048     fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2049     fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2050     fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2051     fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2052     fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2053     fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2054     fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2055     fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2056     fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2057     fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2058     fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2059     fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2060     fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2061     fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2062
2063     fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2064
2065     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2066     fHistMultRemovedNeutrons = new  TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2067     fHistMultRemovedAntiNeutrons = new  TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", nmultbins, -0.5, 99.5, 10, -0.5, 9.5);
2068     /*
2069       fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", nbins, 0, 30, 10, -0.5, 9.5);
2070       fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", nbins, 0, 30, 10, -0.5, 9.5);
2071
2072       fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", nbins, 0, 30, 10, -0.5, 9.5);
2073       fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", nbins, 0, 30, 10, -0.5, 9.5);*/
2074
2075
2076     fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2077     fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2078
2079     fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2080     fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", nmultbins, -0.5, 99.5, nmultbins, -0.5, 99.5);
2081
2082     fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", nmultbins*10, -0.5, 999.5, nmultbins, -0.5, 99.5);
2083     fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", nmultbins*10, -0.5, 999.5, nmultbins, -0.5, 99.5);
2084     fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", nmultbins*10, -0.5, 999.5, 100, -0.5, 99.5);
2085
2086     fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", nmultbins*10, -0.5, 999.5, nmultbins, -0.5, 99.5);
2087     fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", nmultbins*10, -0.5, 999.5, nmultbins, -0.5, 99.5);
2088     fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", nmultbins*10, -0.5, 999.5, nmultbins, -0.5, 99.5);
2089
2090     Int_t ndnbins = 8;
2091     fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", ndnbins, -200, 200, ndnbins, -200, 200);
2092     fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", ndnbins, -200, 200, ndnbins, -200, 200);
2093     fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", ndnbins, -200, 200, ndnbins, -200, 200);
2094     fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", ndnbins, -200, 200, ndnbins, -200, 200);
2095
2096     if(fCalcForKaonCorrection){
2097       Int_t nEtCut = 11;
2098       Float_t etCutAxis[12] = {0.00,0.05,0.10,0.15,0.20,  0.25,0.30,0.35,0.40,0.45, 0.50,0.55};
2099       fHistK0EDepositsVsPtInAcceptance = new TH3F("fHistK0EDepositsVsPtInAcceptance","Kaon deposits with corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
2100       fHistK0EGammaVsPtInAcceptance = new TH3F("fHistK0EGammaVsPtInAcceptance","Kaon deposits without corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
2101       fHistK0EDepositsVsPtOutOfAcceptance = new TH3F("fHistK0EDepositsVsPtOutOfAcceptance","Kaon deposits with corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
2102       fHistK0EGammaVsPtOutOfAcceptance = new TH3F("fHistK0EGammaVsPtOutOfAcceptance","Kaon deposits without corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
2103       fHistK0EDepositsVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
2104       fHistK0EGammaVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
2105       fHistK0EDepositsVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
2106       fHistK0EGammaVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
2107       fHistK0EDepositsVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
2108       fHistK0EGammaVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
2109       fHistK0EDepositsVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
2110       fHistK0EGammaVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
2111       fHistK0EDepositsVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
2112       fHistK0EGammaVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
2113       fHistK0EDepositsVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
2114       fHistK0EGammaVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
2115
2116       fHistSimKaonsInAcceptance = new TH1F("fHistSimKaonsInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
2117       fHistSimK0SInAcceptance = new TH1F("fHistSimK0SInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
2118       fHistSimK0LInAcceptance = new TH1F("fHistSimK0LInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
2119       fHistSimKPlusInAcceptance = new TH1F("fHistSimKPlusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
2120       fHistSimKMinusInAcceptance = new TH1F("fHistSimKMinusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
2121       fHistSimKaonsOutOfAcceptance = new TH1F("fHistSimKaonsOutOfAcceptance","Kaons with y>0.5",fgNumOfPtBins,fgPtAxis);
2122       fHistSimKaonsInAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsInAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
2123       fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries","Secondary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
2124       fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
2125     }
2126
2127     fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", nmultbins, -0.5, 1999.5);
2128     fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", nmultbins, -0.5, 1999.5);
2129     fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", nmultbins, -0.5, 1999.5);
2130
2131     fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", nmultbins, -0.5, 1999.5);
2132     fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", nmultbins, -0.5, 1999.5);
2133     fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", nmultbins, -0.5, 1999.5);
2134
2135     if(fCuts->GetHistMakeTree())
2136     {
2137         TString treename = "fPrimaryTree" + fHistogramNameSuffix;
2138         fPrimaryTree = new TTree(treename, treename);
2139
2140         fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
2141         fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
2142         fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
2143
2144         fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
2145         fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
2146
2147         fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
2148         fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
2149
2150         fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
2151         fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
2152         fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
2153
2154         fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
2155         fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
2156         fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
2157
2158         fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
2159         fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
2160
2161
2162         fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
2163         fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
2164         fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
2165         fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
2166
2167         fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
2168         fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
2169         fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
2170
2171         fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
2172
2173
2174         fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
2175         fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
2176
2177         fPrimaryTree->Branch("fClusterMult", &fClusterMult,  "fClusterMult/I");
2178
2179
2180     }
2181
2182     fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
2183     fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
2184     fHistGammasFoundOutOfAccCent = new TH2F("fHistGammasFoundOutOfAccCent", "fHistGammasFoundOutOfAccCent",200, 0, 10,20,-0.5,19.5);
2185     fHistGammasFoundCent = new TH2F("fHistGammasFoundCent", "fHistGammasFoundCent",200, 0, 10,20,-0.5,19.5);
2186     fHistGammasFoundOutOfAccAltCent = new TH2F("fHistGammasFoundOutOfAccAltCent", "fHistGammasFoundOutOfAccCent",200, 0, 10,20,-0.5,19.5);
2187     fHistGammasFoundRecoEnergyCent = new TH2F("fHistGammasFoundRecoEnergyCent", "fHistGammasFoundRecoEnergyCent",200, 0, 10,20,-0.5,19.5);
2188     fHistAllGammasFoundRecoEnergyCent = new TH2F("fHistAllGammasFoundRecoEnergyCent", "fHistAllGammasFoundRecoEnergyCent",200, 0, 10,20,-0.5,19.5);
2189     fHistGammasFoundOutOfAccRecoEnergyCent = new TH2F("fHistGammasFoundOutOfAccRecoEnergyCent", "fHistGammasFoundOutOfAccRecoEnergyCent",200, 0, 10,20,-0.5,19.5);
2190     fHistAllGammasFoundOutOfAccRecoEnergyCent = new TH2F("fHistAllGammasFoundOutOfAccRecoEnergyCent", "fHistAllGammasFoundOutOfAccRecoEnergyCent",200, 0, 10,20,-0.5,19.5);
2191     fHistGammasFoundAltCent = new TH2F("fHistGammasFoundAltCent", "fHistGammasFoundAltCent",200, 0, 10,20,-0.5,19.5);
2192     fHistGammasGeneratedCent = new TH2F("fHistGammasGeneratedCent", "fHistGammasGeneratedCent",200, 0, 10,20,-0.5,19.5);
2193     fHistChargedTracksCut = new TH1F("fHistChargedTracksCut", "fHistChargedTracksCut",100, 0, 5);
2194     fHistChargedTracksAccepted = new TH1F("fHistChargedTracksAccepted", "fHistChargedTracksAccepted",100, 0, 5);
2195     fHistGammasCut = new TH1F("fHistGammasTracksCut", "fHistGammasTracksCut",100, 0, 5);
2196     fHistGammasAccepted = new TH1F("fHistGammasTracksAccepted", "fHistGammasTracksAccepted",100, 0, 5);
2197
2198     if(fCalcTrackMatchVsMult){
2199       fHistChargedTracksCutMult = new TH2F("fHistChargedTracksCutMult", "fHistChargedTracksCutMult",100, 0, 5,10*scale,0,100*scale);
2200       fHistChargedTracksAcceptedMult = new TH2F("fHistChargedTracksAcceptedMult", "fHistChargedTracksAcceptedMult",100, 0, 5,10*scale,0,100*scale);
2201       fHistChargedTrackDepositsAcceptedVsPt = new TH2F("fHistChargedTrackDepositsAcceptedVsPt", "fHistChargedTrackDepositsAcceptedVsPt",100, 0, 5,20,-0.5,19.5);
2202       fHistChargedTrackDepositsAllVsPt = new TH2F("fHistChargedTrackDepositsAllVsPt", "fHistChargedTrackDepositsAllVsPt",100, 0, 5,20,-0.5,19.5);
2203       fHistChargedTrackDepositsAcceptedVsPtEffCorr = new TH2F("fHistChargedTrackDepositsAcceptedVsPtEffCorr", "fHistChargedTrackDepositsAcceptedVsPtEffCorr",100, 0, 5,20,-0.5,19.5);
2204       fHistChargedTrackDepositsAllVsPtEffCorr = new TH2F("fHistChargedTrackDepositsAllVsPtEffCorr", "fHistChargedTrackDepositsAllVsPtEffCorr",100, 0, 5,20,-0.5,19.5);
2205       fHistChargedTracksAcceptedLowPtCentEffCorr = new TH2F("fHistChargedTracksAcceptedLowPtCentEffCorr", "fHistChargedTracksAcceptedLowPtCentEffCorr",100, 0, 5,20,-0.5,19.5);
2206       fHistChargedTracksAcceptedLowPtCent = new TH2F("fHistChargedTracksAcceptedLowPtCent", "fHistChargedTracksAcceptedLowPtCent",100, 0, 5,20,-0.5,19.5);
2207       fHistChargedTracksAcceptedLowPtCent500MeV = new TH2F("fHistChargedTracksAcceptedLowPtCent500MeV", "fHistChargedTracksAcceptedLowPtCent500MeV",100, 0, 5,20,-0.5,19.5);
2208       fHistChargedTracksAcceptedLowPtCentNoAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentNoAntiProtons", "fHistChargedTracksAcceptedLowPtCentNoAntiProtons",100, 0, 5,20,-0.5,19.5);
2209       fHistChargedTracksAcceptedLowPtCentAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentAntiProtons", "fHistChargedTracksAcceptedLowPtCentAntiProtons",100, 0, 5,20,-0.5,19.5);
2210       fHistGammasCutMult = new TH2F("fHistGammasTracksCutMult", "fHistGammasTracksCutMult",100, 0, 5,10*scale,0,100*scale);
2211       fHistGammasAcceptedMult = new TH2F("fHistGammasTracksAcceptedMult", "fHistGammasTracksAcceptedMult",100, 0, 5,10*scale,0,100*scale);
2212     }
2213
2214     fHistBadTrackMatches = new TH1F("fHistBadTrackMatches", "fHistBadTrackMatches",100, 0, 5);
2215     fHistMatchedTracksEvspTBkgd = new TH2F("fHistMatchedTracksEvspTBkgd", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
2216     fHistMatchedTracksEvspTSignal = new TH2F("fHistMatchedTracksEvspTSignal", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
2217     if(fCalcTrackMatchVsMult){
2218       fHistMatchedTracksEvspTBkgdPeripheral = new TH2F("fHistMatchedTracksEvspTBkgdPeripheral", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
2219       fHistMatchedTracksEvspTSignalPeripheral = new TH2F("fHistMatchedTracksEvspTSignalPeripheral", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
2220
2221       fHistMatchedTracksEvspTBkgdvsCent = new TH3F("fHistMatchedTracksEvspTBkgdvsCent", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
2222       fHistMatchedTracksEvspTSignalvsCent = new TH3F("fHistMatchedTracksEvspTSignalvsCent", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
2223       fHistMatchedTracksEvspTBkgdvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTBkgdvsCentEffCorr", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
2224       fHistMatchedTracksEvspTSignalvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTSignalvsCentEffCorr", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
2225
2226
2227       fHistChargedTracksCutPeripheral = new TH1F("fHistChargedTracksCutPeripheral", "fHistChargedTracksCut",100, 0, 5);
2228       fHistChargedTracksAcceptedPeripheral = new TH1F("fHistChargedTracksAcceptedPeripheral", "fHistChargedTracksAccepted",100, 0, 5);
2229       fHistGammasCutPeripheral = new TH1F("fHistGammasTracksCutPeripheral", "fHistGammasTracksCut",100, 0, 5);
2230       fHistGammasAcceptedPeripheral = new TH1F("fHistGammasTracksAcceptedPeripheral", "fHistGammasTracksAccepted",100, 0, 5);
2231     }
2232     fHistBadTrackMatchesdPhidEta = new TH2F("fHistBadTrackMatchesdPhidEta", "fHistBadTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
2233     fHistGoodTrackMatchesdPhidEta = new TH2F("fHistGoodTrackMatchesdPhidEta", "fHistGoodTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
2234
2235     fHistHadronDepositsAll = new TH1F("fHistHadronDepositsAll","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
2236     fHistHadronDepositsReco = new TH1F("fHistHadronDepositsReco","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
2237     //,10,0,100
2238       Int_t nMult = 50;
2239       Float_t nMultCuts[51] = { 0  , 5 ,10 ,15 ,20 , 25, 30, 35, 40, 45,
2240                                 50 ,55 ,60 ,65 ,70 , 75, 80, 85, 90, 95,
2241                                 100,105,110,115,120,125,130,135,140,145,
2242                                 150,155,160,165,170,175,180,185,190,195,
2243                                 200,205,210,215,220,225,230,235,240,245,
2244                                 250};
2245       Int_t nCent = 20;
2246       Float_t nCentCuts[21] = { 0, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
2247
2248       fHistHadronDepositsAllCent = new TH2F("fHistHadronDepositsAllCent","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2249       fHistHadronDepositsAllCent500MeV = new TH2F("fHistHadronDepositsAllCent500MeV","All Hadrons which deposited energy in calorimeter pT>500MeV",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2250       fHistHadronDepositsRecoCent = new TH2F("fHistHadronDepositsRecoCent","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2251
2252       fHistHadronDepositsAllvsECent = new TH2F("fHistHadronDepositsAllvsECent","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2253       fHistHadronDepositsRecovsECent = new TH2F("fHistHadronDepositsRecovsECent","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2254
2255       fHistHadronsAllCent = new TH2F("fHistHadronsAllCent","All Hadrons vs cluster mult",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
2256
2257       fHistMultChVsSignalVsMult = new TH3F("fHistMultChVsSignalVsMult","Charged particle Multiplicity vs Signal particle multiplicity vs Cluster Mult",nMult,nMultCuts,nMult,nMultCuts,nMult,nMultCuts);
2258       fHistNeutralRemovedSecondaryEtVsCent = new TH2F("fHistNeutralRemovedSecondaryEtVsCent","Neutral Removed Secondaries E_{T} vs centrality",100*scale,0.0,2.0*scale,20,-0.5,19.5);
2259       fHistChargedRemovedSecondaryEtVsCent = new TH2F("fHistChargedRemovedSecondaryEtVsCent","Charged Removed Secondaries E_{T} vs centrality",100*scale,0.0,2.0*scale,20,-0.5,19.5);
2260       fHistNeutralNotRemovedSecondaryEtVsCent = new TH2F("fHistNeutralNotRemovedSecondaryEtVsCent","Neutral NotRemoved Secondaries E_{T} vs centrality",100*scale,0.0,2.0*scale,20,-0.5,19.5);
2261       fHistChargedNotRemovedSecondaryEtVsCent = new TH2F("fHistChargedNotRemovedSecondaryEtVsCent","Charged NotRemoved Secondaries E_{T} vs centrality",100*scale,0.0,2.0*scale,20,-0.5,19.5);
2262       fHistNeutralRemovedSecondaryNumVsNCluster = new TH2F("fHistNeutralRemovedSecondaryNumVsNCluster","Neutral Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,250*scale,0,250*scale);
2263       fHistChargedRemovedSecondaryNumVsNCluster = new TH2F("fHistChargedRemovedSecondaryNumVsNCluster","Charged Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,250*scale,0,250*scale);
2264       fHistNeutralNotRemovedSecondaryNumVsNCluster = new TH2F("fHistNeutralNotRemovedSecondaryNumVsNCluster","Neutral NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,250*scale,0,250*scale);
2265       fHistChargedNotRemovedSecondaryNumVsNCluster = new TH2F("fHistChargedNotRemovedSecondaryNumVsNCluster","Charged NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,250*scale,0,250*scale);
2266       fHistNeutralRemovedSecondaryNumVsCent = new TH2F("fHistNeutralRemovedSecondaryNumVsCent","Neutral Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
2267       fHistChargedRemovedSecondaryNumVsCent = new TH2F("fHistChargedRemovedSecondaryNumVsCent","Charged Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
2268       fHistNeutralNotRemovedSecondaryNumVsCent = new TH2F("fHistNeutralNotRemovedSecondaryNumVsCent","Neutral NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
2269       fHistChargedNotRemovedSecondaryNumVsCent = new TH2F("fHistChargedNotRemovedSecondaryNumVsCent","Charged NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
2270       fHistNeutronsEtVsCent = new TH2F("fHistNeutronsEtVsCent","Neutrons and anti-neutrons - deposited ET vs Centrality bin",100*scale,0,4.0*scale,20,-0.5,19.5);
2271       fHistNeutronsNumVsCent = new TH2F("fHistNeutronsNumVsCent","Neutrons and anti-neutrons - number vs Centrality bin",20,-0.5,19.5,20,-0.5,19.5);
2272       fHistNotNeutronsNumVsCent = new TH2F("fHistNotNeutronsNumVsCent","Neutral particles not otherwise classified - number vs Centrality bin",20,-0.5,19.5,20,-0.5,19.5);
2273       Int_t nbinsEt = 125*scale;
2274       Float_t maxEtRange = 125*scale;
2275       Float_t maxEtRangeShort = 50*scale;
2276       Float_t minEtRange = 0;
2277       Int_t nbinsMult = 100*scale;
2278       Float_t maxMult = 3000*scale;
2279       Float_t minMult = 0;
2280       Int_t nbinsCl = 250*scale;
2281       Float_t maxCl = 500*scale;
2282       Float_t minCl = 0;
2283       fHistPiKPDepositedVsNch = new TH2F("fHistPiKPDepositedVsNch","#pi,K,p E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
2284       fHistPiKPNotTrackMatchedDepositedVsNch = new TH2F("fHistPiKPNotTrackMatchedDepositedVsNch","#pi,K,p E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
2285
2286     fHistNeutronsDepositedVsNch = new TH2F("fHistNeutronsDepositedVsNch","n deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2287     fHistAntiNeutronsDepositedVsNch = new TH2F("fHistAntiNeutronsDepositedVsNch","#bar{n} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2288     fHistProtonsDepositedVsNch = new TH2F("fHistProtonsDepositedVsNch","p deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2289     fHistAntiProtonsDepositedVsNch = new TH2F("fHistAntiProtonsDepositedVsNch","#bar{p} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2290     fHistProtonsNotTrackMatchedDepositedVsNch = new TH2F("fHistProtonsNotTrackMatchedDepositedVsNch","p not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2291     fHistAntiProtonsNotTrackMatchedDepositedVsNch = new TH2F("fHistAntiProtonsNotTrackMatchedDepositedVsNch","#bar{p} not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2292
2293
2294
2295     fHistNeutronsDepositedVsNcl = new TH2F("fHistNeutronsDepositedVsNcl","n deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2296     fHistAntiNeutronsDepositedVsNcl = new TH2F("fHistAntiNeutronsDepositedVsNcl","#bar{n} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2297     fHistProtonsDepositedVsNcl = new TH2F("fHistProtonsDepositedVsNcl","p deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2298     fHistAntiProtonsDepositedVsNcl = new TH2F("fHistAntiProtonsDepositedVsNcl","#bar{p} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2299     fHistProtonsNotTrackMatchedDepositedVsNcl = new TH2F("fHistProtonsNotTrackMatchedDepositedVsNcl","p not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2300     fHistAntiProtonsNotTrackMatchedDepositedVsNcl = new TH2F("fHistAntiProtonsNotTrackMatchedDepositedVsNcl","#bar{p} not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2301
2302     fHistSecondariesVsNch = new TH2F("fHistSecondariesVsNch","secondaries deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2303     fHistSecondariesVsNcl = new TH2F("fHistSecondariesVsNcl","secondaries deposited in calorimeter vs number of clusters",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2304     fHistSecondariesEffCorrVsNch = new TH2F("fHistSecondariesEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2305     fHistSecondariesEffCorrVsNcl = new TH2F("fHistSecondariesEffCorrVsNcl","efficiency corrected secondaries deposited in calorimeter vs number of clusters",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2306
2307     fHistSecondariesOutOfAccEffCorrVsNch = new TH2F("fHistSecondariesOutOfAccEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs number of clusters for secondary particles out of detector acceptance",nbinsEt,minEtRange,maxEtRangeShort/10.0,nbinsMult,minMult,maxMult);
2308     fHistSecondariesDetectorCoverEffCorrVsNch = new TH2F("fHistSecondariesDetectorCoverEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs number of clusters for secondaries from the detector cover",nbinsEt,minEtRange,maxEtRangeShort/10.0,nbinsMult,minMult,maxMult);
2309     //start
2310     fHistNeutronsDepositedVsNchNoEffCorr = new TH2F("fHistNeutronsDepositedVsNchNoEffCorr","n deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2311     fHistAntiNeutronsDepositedVsNchNoEffCorr = new TH2F("fHistAntiNeutronsDepositedVsNchNoEffCorr","#bar{n} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2312     fHistProtonsDepositedVsNchNoEffCorr = new TH2F("fHistProtonsDepositedVsNchNoEffCorr","p deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2313     fHistAntiProtonsDepositedVsNchNoEffCorr = new TH2F("fHistAntiProtonsDepositedVsNchNoEffCorr","#bar{p} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2314     fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr = new TH2F("fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr","p not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2315     fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr = new TH2F("fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr","#bar{p} not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
2316
2317
2318
2319     fHistNeutronsDepositedVsNclNoEffCorr = new TH2F("fHistNeutronsDepositedVsNclNoEffCorr","n deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2320     fHistAntiNeutronsDepositedVsNclNoEffCorr = new TH2F("fHistAntiNeutronsDepositedVsNclNoEffCorr","#bar{n} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2321     fHistProtonsDepositedVsNclNoEffCorr = new TH2F("fHistProtonsDepositedVsNclNoEffCorr","p deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2322     fHistAntiProtonsDepositedVsNclNoEffCorr = new TH2F("fHistAntiProtonsDepositedVsNclNoEffCorr","#bar{p} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2323     fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr = new TH2F("fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr","p not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2324     fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr = new TH2F("fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr","#bar{p} not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
2325
2326     //end
2327
2328     fHistCentVsNchVsNcl = new TH3F("fHistCentVsNchVsNcl","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
2329     float maxpos = 500*scale;
2330      int nbinspos = 50*scale;
2331     fHistSecondaryPositionInDetector = new TH3F("fHistSecondaryPositionInDetector","Position of secondaries",nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos);
2332     fHistSecondaryPositionInDetector->GetXaxis()->SetTitle("X");
2333     fHistSecondaryPositionInDetector->GetYaxis()->SetTitle("Y");
2334     fHistSecondaryPositionInDetector->GetZaxis()->SetTitle("Z");
2335 //     fHistSecondaryPositionInDetectorMultiple = new TH3F("fHistSecondaryPositionInDetectorMultiple","Position of secondaries",nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos);
2336 //     fHistSecondaryPositionInDetectorMultiple->GetXaxis()->SetTitle("X");
2337 //     fHistSecondaryPositionInDetectorMultiple->GetYaxis()->SetTitle("Y");
2338 //     fHistSecondaryPositionInDetectorMultiple->GetZaxis()->SetTitle("Z");
2339     fClusterPositionWeird =  new TH2F("fClusterPositionWeird", "Position of weird secondary clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
2340
2341    fSecondaryClusterEnergy = new TH1F("fSecondaryClusterEnergy","fSecondaryClusterEnergy", 100, 0, 5);
2342    fHistGammaCrossCheck =  new TH2F("fHistGammaCrossCheck", "(E_{T}^{#gamma,rec}-E_{T}^{#gamma,sim})/E_{T}^{#gamma,rec}",20,-0.5,19.5,100,-1.5,4);
2343    fHistGammaCrossCheckAlt =  new TH2F("fHistGammaCrossCheckAlt", "(E_{T}^{#gamma,rec}-E_{T}^{#gamma,sim})/E_{T}^{#gamma,rec}",20,-0.5,19.5,100,-1.5,4);
2344    fHistGammaEnergyCrossCheck =  new TH2F("fHistGammaEnergyCrossCheck", "(E_{T}^{#gamma,rec}-E_{T}^{#gamma,sim})/E_{T}^{#gamma,rec}",100,0.0,10.0,100,-2,2);
2345    fHistGammaEnergyCrossCheckAlt =  new TH2F("fHistGammaEnergyCrossCheckAlt", "(E_{T}^{#gamma,rec}-E_{T}^{#gamma,sim})/E_{T}^{#gamma,rec}",100,0.0,10.0,100,-2,2);
2346    fHistNeutronCrossCheck =  new TH2F("fHistNeutronCrossCheck", "(E_{T}^{n,rec}-E_{T}^{n,sim})/E_{T}^{n,rec}",20,-0.5,19.5,100,-14.8,1.2);
2347    fHistSecondaryCrossCheck =  new TH2F("fHistSecondaryCrossCheck", "(E_{T}^{sec,rec}-E_{T}^{sec,sim})/E_{T}^{sec,rec}",20,-0.5,19.5,100,-14.8,1.2);
2348    fHistHadronCrossCheck =  new TH2F("fHistHadronCrossCheck", "(E_{T}^{h,rec}-E_{T}^{h,sim})/E_{T}^{h,rec}",20,-0.5,19.5,100,-7.8,1.2);
2349    fHistKaonCrossCheck =  new TH2F("fHistKaonCrossCheck", "(E_{T}^{K,rec}-E_{T}^{K,sim})/E_{T}^{K,rec}",20,-0.5,19.5,100,-7.8,1.2);
2350
2351    fHistAllEnergy = new TH1F("fHistAllEnergy","fHistAllEnergy",20,-0.5,19.5);
2352    fHistSignalEnergy = new TH1F("fHistSignalEnergy","fHistSignalEnergy",20,-0.5,19.5);
2353    fHistNeutronEnergy = new TH1F("fHistNeutronEnergy","fHistNeutronEnergy",20,-0.5,19.5);
2354    fHistKaonEnergy = new TH1F("fHistKaonEnergy","fHistKaonEnergy",20,-0.5,19.5);
2355    fHistHadronEnergy = new TH1F("fHistHadronEnergy","fHistHadronEnergy",20,-0.5,19.5);
2356    fHistSecondaryEnergy = new TH1F("fHistSecondaryEnergy","fHistSecondaryEnergy",20,-0.5,19.5);
2357    fHistSecondaryChargedEnergy = new TH1F("fHistSecondaryChargedEnergy","fHistSecondaryChargedEnergy",20,-0.5,19.5);
2358    fHistSecondaryNeutronEnergy = new TH1F("fHistSecondaryNeutronEnergy","fHistSecondaryNeutronEnergy",20,-0.5,19.5);
2359    fHistSecondaryGammaEnergy = new TH1F("fHistSecondaryGammaEnergy","fHistSecondaryGammaEnergy",20,-0.5,19.5);
2360    fHistSecondaryElectronEnergy = new TH1F("fHistSecondaryElectronEnergy","fHistSecondaryElectronEnergy",20,-0.5,19.5);
2361    fHistSecondaryOtherEnergy = new TH1F("fHistSecondaryOtherEnergy","fHistSecondaryOtherEnergy",20,-0.5,19.5);
2362    fHistSimulatedGammaEnergy = new TH1F("fHistSimulatedGammaEnergy","fHistSimulatedGammaEnergy",20,-0.5,19.5);
2363    fHistReconstructedGammaEnergy = new TH1F("fHistReconstructedGammaEnergy","fHistReconstructedGammaEnergy",20,-0.5,19.5);
2364    fHistSimulatedGammaEnergyAboveThreshold = new TH1F("fHistSimulatedGammaEnergyAboveThreshold","fHistSimulatedGammaEnergyAboveThreshold",20,-0.5,19.5);
2365    fHistReconstructedSignalEnergy = new TH1F("fHistReconstructedSignalEnergy","fHistReconstructedSignalEnergy",20,-0.5,19.5);
2366
2367    Float_t fracMin = 0.0;
2368    Float_t fracMax = 1.0;
2369    Float_t fracMaxLow = 0.2;
2370    Float_t nBinsFrac = 50;
2371    fHistFracSignalVsNClusters = new TH2F("fHistFracSignalVsNClusters","fHistFracSignalVsNClusters",nbinsCl,minCl,maxCl,nBinsFrac,fracMin,fracMax);
2372    fHistFracHadronsVsNClusters = new TH2F("fHistFracHadronsVsNClusters","fHistFracHadronsVsNClusters",nbinsCl,minCl,maxCl,nBinsFrac,fracMin,fracMax);
2373    fHistFracNeutronsVsNClusters = new TH2F("fHistFracNeutronsVsNClusters","fHistFracNeutronsVsNClusters",nbinsCl,minCl,maxCl,nBinsFrac,fracMin,fracMaxLow);
2374    fHistFracKaonsVsNClusters = new TH2F("fHistFracKaonsVsNClusters","fHistFracKaonsVsNClusters",nbinsCl,minCl,maxCl,nBinsFrac,fracMin,fracMaxLow);
2375    fHistFracSecondariesVsNClusters = new TH2F("fHistFracSecondariesVsNClusters","fHistFracSecondariesVsNClusters",nbinsCl,minCl,maxCl,nBinsFrac,fracMin,fracMax);
2376    fHistFracSignalVsNMultiplicity = new TH2F("fHistFracSignalVsNMultiplicity","fHistFracSignalVsNMultiplicity",nbinsMult,minMult,maxMult,nBinsFrac,fracMin,fracMax);
2377    fHistFracHadronsVsNMultiplicity = new TH2F("fHistFracHadronsVsNMultiplicity","fHistFracHadronsVsNMultiplicity",nbinsMult,minMult,maxMult,nBinsFrac,fracMin,fracMax);
2378    fHistFracNeutronsVsNMultiplicity = new TH2F("fHistFracNeutronsVsNMultiplicity","fHistFracNeutronsVsNMultiplicity",nbinsMult,minMult,maxMult,nBinsFrac,fracMin,fracMaxLow);
2379    fHistFracKaonsVsNMultiplicity = new TH2F("fHistFracKaonsVsNMultiplicity","fHistFracKaonsVsNMultiplicity",nbinsMult,minMult,maxMult,nBinsFrac,fracMin,fracMaxLow);
2380    fHistFracSecondariesVsNMultiplicity = new TH2F("fHistFracSecondariesVsNMultiplicity","fHistFracSecondariesVsNMultiplicity",nbinsMult,minMult,maxMult,nBinsFrac,fracMin,fracMax);
2381    Int_t nBinsMatchedTracks = 100*scale;
2382    Float_t lowMatchedTracks = 0;
2383    Float_t highMatchedTracks = 100*scale;
2384    fHistFracSignalVsNMatchedTracks = new TH2F("fHistFracSignalVsNMatchedTracks","fHistFracSignalVsNMatchedTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2385    fHistFracHadronsVsNMatchedTracks = new TH2F("fHistFracHadronsVsNMatchedTracks","fHistFracHadronsVsNMatchedTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2386    fHistFracNeutronsVsNMatchedTracks = new TH2F("fHistFracNeutronsVsNMatchedTracks","fHistFracNeutronsVsNMatchedTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMaxLow);
2387    fHistFracKaonsVsNMatchedTracks = new TH2F("fHistFracKaonsVsNMatchedTracks","fHistFracKaonsVsNMatchedTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMaxLow);
2388    fHistFracSecondariesVsNMatchedTracks = new TH2F("fHistFracSecondariesVsNMatchedTracks","fHistFracSecondariesVsNMatchedTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2389    fHistFracSignalVsNTotalTracks = new TH2F("fHistFracSignalVsNTotalTracks","fHistFracSignalVsNTotalTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2390    fHistFracHadronsVsNTotalTracks = new TH2F("fHistFracHadronsVsNTotalTracks","fHistFracHadronsVsNTotalTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2391    fHistFracNeutronsVsNTotalTracks = new TH2F("fHistFracNeutronsVsNTotalTracks","fHistFracNeutronsVsNTotalTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMaxLow);
2392    fHistFracKaonsVsNTotalTracks = new TH2F("fHistFracKaonsVsNTotalTracks","fHistFracKaonsVsNTotalTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMaxLow);
2393    fHistFracSecondariesVsNTotalTracks = new TH2F("fHistFracSecondariesVsNTotalTracks","fHistFracSecondariesVsNTotalTracks",nBinsMatchedTracks,lowMatchedTracks,highMatchedTracks,nBinsFrac,fracMin,fracMax);
2394    fHistRCorrVsPtVsCent = new TH3F("fHistRCorrVsPtVsCent","fHistRCorrVsPtVsCent",72,0,2,50,0,10,20,-0.5,19.5);
2395    //NClusters
2396 }
2397
2398 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
2399 {   //fill the output list
2400     AliAnalysisEt::FillOutputList(list);
2401
2402   if(!fIsMC) return;
2403     if(fCuts->GetHistMakeTree())
2404     {
2405         list->Add(fPrimaryTree);
2406     }
2407
2408     list->Add(fHistRemovedOrNot);
2409
2410     list->Add(fHistEtNonRemovedProtons);
2411     list->Add(fHistEtNonRemovedAntiProtons);
2412     list->Add(fHistEtNonRemovedPiPlus);
2413     list->Add(fHistEtNonRemovedPiMinus);
2414     list->Add(fHistEtNonRemovedKaonPlus);
2415     list->Add(fHistEtNonRemovedKaonMinus);
2416     list->Add(fHistEtNonRemovedK0s);
2417     list->Add(fHistEtNonRemovedK0L);
2418     list->Add(fHistEtNonRemovedLambdas);
2419     list->Add(fHistEtNonRemovedElectrons);
2420     list->Add(fHistEtNonRemovedPositrons);
2421     list->Add(fHistEtNonRemovedMuPlus);
2422     list->Add(fHistEtNonRemovedMuMinus);
2423     list->Add(fHistEtNonRemovedNeutrons);
2424     list->Add(fHistEtNonRemovedAntiNeutrons);
2425     list->Add(fHistEtNonRemovedGammas);
2426     list->Add(fHistEtNonRemovedGammasFromPi0);
2427
2428     list->Add(fHistEtRemovedGammas);
2429     list->Add(fHistEtRemovedNeutrons);
2430     list->Add(fHistEtRemovedAntiNeutrons);
2431
2432     list->Add(fHistEtRemovedCharged);
2433     list->Add(fHistEtRemovedNeutrals);
2434
2435     list->Add(fHistEtNonRemovedCharged);
2436     list->Add(fHistEtNonRemovedNeutrals);
2437
2438     list->Add(fHistMultNonRemovedProtons);
2439     list->Add(fHistMultNonRemovedAntiProtons);
2440     list->Add(fHistMultNonRemovedPiPlus);
2441     list->Add(fHistMultNonRemovedPiMinus);
2442     list->Add(fHistMultNonRemovedKaonPlus);
2443     list->Add(fHistMultNonRemovedKaonMinus);
2444     list->Add(fHistMultNonRemovedK0s);
2445     list->Add(fHistMultNonRemovedK0L);
2446     list->Add(fHistMultNonRemovedLambdas);
2447     list->Add(fHistMultNonRemovedElectrons);
2448     list->Add(fHistMultNonRemovedPositrons);
2449     list->Add(fHistMultNonRemovedMuPlus);
2450     list->Add(fHistMultNonRemovedMuMinus);
2451     list->Add(fHistMultNonRemovedNeutrons);
2452     list->Add(fHistMultNonRemovedAntiNeutrons);
2453     list->Add(fHistMultNonRemovedGammas);
2454
2455     list->Add(fHistMultRemovedGammas);
2456     list->Add(fHistMultRemovedNeutrons);
2457     list->Add(fHistMultRemovedAntiNeutrons);
2458
2459     list->Add(fHistMultRemovedCharged);
2460     list->Add(fHistMultRemovedNeutrals);
2461
2462     list->Add(fHistMultNonRemovedCharged);
2463     list->Add(fHistMultNonRemovedNeutrals);
2464
2465     list->Add(fHistTrackMultvsNonRemovedCharged);
2466     list->Add(fHistTrackMultvsNonRemovedNeutral);
2467     list->Add(fHistTrackMultvsRemovedGamma);
2468
2469     list->Add(fHistClusterMultvsNonRemovedCharged);
2470     list->Add(fHistClusterMultvsNonRemovedNeutral);
2471     list->Add(fHistClusterMultvsRemovedGamma);
2472
2473     //list->Add(fHistDecayVertexNonRemovedCharged);
2474     //list->Add(fHistDecayVertexNonRemovedNeutral);
2475     //list->Add(fHistDecayVertexRemovedCharged);
2476     //list->Add(fHistDecayVertexRemovedNeutral);
2477
2478     list->Add(fHistDxDzNonRemovedCharged);
2479     list->Add(fHistDxDzRemovedCharged);
2480     list->Add(fHistDxDzNonRemovedNeutral);
2481     list->Add(fHistDxDzRemovedNeutral);
2482
2483     if(fCalcForKaonCorrection){
2484       list->Add(fHistK0EDepositsVsPtInAcceptance);
2485       list->Add(fHistK0EGammaVsPtInAcceptance);
2486       list->Add(fHistK0EDepositsVsPtOutOfAcceptance);
2487       list->Add(fHistK0EGammaVsPtOutOfAcceptance);
2488       list->Add(fHistSimKaonsInAcceptance);
2489       list->Add(fHistSimK0SInAcceptance);
2490       list->Add(fHistSimK0LInAcceptance);
2491       list->Add(fHistSimKPlusInAcceptance);
2492       list->Add(fHistSimKMinusInAcceptance);
2493       list->Add(fHistSimKaonsOutOfAcceptance);
2494       list->Add(fHistSimKaonsInAcceptanceWithDepositsPrimaries);
2495       list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries);
2496       list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries);
2497     }
2498
2499     list->Add(fHistPiPlusMult);
2500     list->Add(fHistPiMinusMult);
2501     list->Add(fHistPiZeroMult);
2502     list->Add(fHistPiPlusMultAcc);
2503     list->Add(fHistPiMinusMultAcc);
2504     list->Add(fHistPiZeroMultAcc);
2505
2506     list->Add(fHistGammasFound);
2507     list->Add(fHistGammasGenerated);
2508     list->Add(fHistGammasFoundOutOfAccCent);
2509     list->Add(fHistGammasFoundCent);
2510     list->Add(fHistGammasFoundOutOfAccAltCent);
2511     list->Add(fHistGammasFoundAltCent);
2512     list->Add(fHistGammasGeneratedCent);
2513     list->Add(fHistGammasFoundRecoEnergyCent);
2514     list->Add(fHistAllGammasFoundRecoEnergyCent);
2515     list->Add(fHistGammasFoundOutOfAccRecoEnergyCent);
2516     list->Add(fHistAllGammasFoundOutOfAccRecoEnergyCent);
2517     list->Add(fHistChargedTracksCut);
2518     list->Add(fHistChargedTracksAccepted);
2519     list->Add(fHistGammasCut);
2520     list->Add(fHistGammasAccepted);
2521     if(fCalcTrackMatchVsMult){
2522       list->Add(fHistChargedTracksCutMult);
2523       list->Add(fHistChargedTrackDepositsAcceptedVsPt);
2524       list->Add(fHistChargedTrackDepositsAllVsPt);
2525       list->Add(fHistChargedTrackDepositsAcceptedVsPtEffCorr);
2526       list->Add(fHistChargedTrackDepositsAllVsPtEffCorr);
2527       list->Add(fHistChargedTracksAcceptedMult);
2528       list->Add(fHistChargedTracksAcceptedLowPtCentEffCorr);
2529       list->Add(fHistChargedTracksAcceptedLowPtCent);
2530       list->Add(fHistChargedTracksAcceptedLowPtCent500MeV);
2531       list->Add(fHistChargedTracksAcceptedLowPtCentNoAntiProtons);
2532       list->Add(fHistChargedTracksAcceptedLowPtCentAntiProtons);
2533       list->Add(fHistGammasCutMult);
2534       list->Add(fHistGammasAcceptedMult);
2535     }
2536     list->Add(fHistBadTrackMatches);
2537     list->Add(fHistMatchedTracksEvspTBkgd);
2538     list->Add(fHistMatchedTracksEvspTSignal);
2539     if(fCalcTrackMatchVsMult){
2540       list->Add(fHistMatchedTracksEvspTBkgdPeripheral);
2541       list->Add(fHistMatchedTracksEvspTSignalPeripheral);
2542       list->Add(fHistMatchedTracksEvspTBkgdvsCent);
2543       list->Add(fHistMatchedTracksEvspTSignalvsCent);
2544       list->Add(fHistMatchedTracksEvspTBkgdvsCentEffCorr);
2545       list->Add(fHistMatchedTracksEvspTSignalvsCentEffCorr);
2546       list->Add(fHistChargedTracksCutPeripheral);
2547       list->Add(fHistChargedTracksAcceptedPeripheral);
2548       list->Add(fHistGammasCutPeripheral);
2549       list->Add(fHistGammasAcceptedPeripheral);
2550     }
2551     list->Add(fHistBadTrackMatchesdPhidEta);
2552     list->Add(fHistGoodTrackMatchesdPhidEta);
2553     list->Add(fHistHadronDepositsAll);
2554     list->Add(fHistHadronDepositsReco);
2555     list->Add(fHistHadronDepositsAllCent);
2556     list->Add(fHistHadronDepositsAllCent500MeV);
2557     list->Add(fHistHadronDepositsRecoCent);
2558     list->Add(fHistHadronDepositsAllvsECent);
2559     list->Add(fHistHadronDepositsRecovsECent);
2560     list->Add(fHistHadronsAllCent);
2561     list->Add(fHistMultChVsSignalVsMult);
2562     list->Add(fHistNeutralRemovedSecondaryEtVsCent);
2563     list->Add(fHistChargedRemovedSecondaryEtVsCent);
2564     list->Add(fHistNeutralNotRemovedSecondaryEtVsCent);
2565     list->Add(fHistChargedNotRemovedSecondaryEtVsCent);
2566     list->Add(fHistNeutralRemovedSecondaryNumVsNCluster);
2567     list->Add(fHistChargedRemovedSecondaryNumVsNCluster);
2568     list->Add(fHistNeutralNotRemovedSecondaryNumVsNCluster);
2569     list->Add(fHistChargedNotRemovedSecondaryNumVsNCluster);
2570     list->Add(fHistNeutralRemovedSecondaryNumVsCent);
2571     list->Add(fHistChargedRemovedSecondaryNumVsCent);
2572     list->Add(fHistNeutralNotRemovedSecondaryNumVsCent);
2573     list->Add(fHistChargedNotRemovedSecondaryNumVsCent);
2574     list->Add(fHistNeutronsEtVsCent);
2575     list->Add(fHistNeutronsNumVsCent);
2576     list->Add(fHistNotNeutronsNumVsCent);
2577     list->Add(fHistPiKPDepositedVsNch);
2578     list->Add(fHistPiKPNotTrackMatchedDepositedVsNch);
2579
2580     list->Add(fHistNeutronsDepositedVsNch);
2581     list->Add(fHistAntiNeutronsDepositedVsNch);
2582     list->Add(fHistProtonsDepositedVsNch);
2583     list->Add(fHistAntiProtonsDepositedVsNch);
2584     list->Add(fHistProtonsNotTrackMatchedDepositedVsNch);
2585     list->Add(fHistAntiProtonsNotTrackMatchedDepositedVsNch);
2586     list->Add(fHistNeutronsDepositedVsNcl);
2587     list->Add(fHistAntiNeutronsDepositedVsNcl);
2588     list->Add(fHistProtonsDepositedVsNcl);
2589     list->Add(fHistAntiProtonsDepositedVsNcl);
2590     list->Add(fHistProtonsNotTrackMatchedDepositedVsNcl);
2591     list->Add(fHistAntiProtonsNotTrackMatchedDepositedVsNcl);
2592     list->Add(fHistSecondariesVsNch);
2593     list->Add(fHistSecondariesVsNcl);
2594     list->Add(fHistSecondariesEffCorrVsNch);
2595     list->Add(fHistSecondariesEffCorrVsNcl);
2596     list->Add(fHistSecondariesOutOfAccEffCorrVsNch);
2597     list->Add(fHistSecondariesDetectorCoverEffCorrVsNch);
2598     //start
2599
2600     list->Add(fHistNeutronsDepositedVsNchNoEffCorr);
2601     list->Add(fHistAntiNeutronsDepositedVsNchNoEffCorr);
2602     list->Add(fHistProtonsDepositedVsNchNoEffCorr);
2603     list->Add(fHistAntiProtonsDepositedVsNchNoEffCorr);
2604     list->Add(fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr);
2605     list->Add(fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr);
2606     list->Add(fHistNeutronsDepositedVsNclNoEffCorr);
2607     list->Add(fHistAntiNeutronsDepositedVsNclNoEffCorr);
2608     list->Add(fHistProtonsDepositedVsNclNoEffCorr);
2609     list->Add(fHistAntiProtonsDepositedVsNclNoEffCorr);
2610     list->Add(fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr);
2611     list->Add(fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr);
2612
2613     //end
2614
2615     list->Add(fHistCentVsNchVsNcl);
2616     list->Add(fHistSecondaryPositionInDetector);
2617     list->Add(fClusterPositionWeird);
2618     //list->Add(fHistSecondaryPositionInDetectorMultiple);
2619     list->Add(fSecondaryClusterEnergy);
2620     list->Add(fHistGammaCrossCheck);
2621     list->Add(fHistGammaCrossCheckAlt);
2622     list->Add(fHistGammaEnergyCrossCheck);
2623     list->Add(fHistGammaEnergyCrossCheckAlt);
2624     list->Add(fHistNeutronCrossCheck);
2625     list->Add(fHistSecondaryCrossCheck);
2626     list->Add(fHistHadronCrossCheck);
2627     list->Add(fHistKaonCrossCheck);
2628
2629
2630     list->Add(fHistAllEnergy);
2631     list->Add(fHistSignalEnergy);
2632     list->Add(fHistNeutronEnergy);
2633     list->Add(fHistKaonEnergy);
2634     list->Add(fHistHadronEnergy);
2635     list->Add(fHistSecondaryEnergy);
2636     list->Add(fHistSecondaryChargedEnergy);
2637     list->Add(fHistSecondaryNeutronEnergy);
2638     list->Add(fHistSecondaryGammaEnergy);
2639     list->Add(fHistSecondaryElectronEnergy);
2640     list->Add(fHistSecondaryOtherEnergy);
2641     list->Add(fHistSimulatedGammaEnergy);
2642     list->Add(fHistReconstructedGammaEnergy);
2643     list->Add(fHistSimulatedGammaEnergyAboveThreshold);
2644     list->Add(fHistReconstructedSignalEnergy);
2645
2646
2647     list->Add(fHistFracSignalVsNClusters);
2648     list->Add(fHistFracHadronsVsNClusters);
2649     list->Add(fHistFracNeutronsVsNClusters);
2650     list->Add(fHistFracKaonsVsNClusters);
2651     list->Add(fHistFracSecondariesVsNClusters);
2652     list->Add(fHistFracSignalVsNMultiplicity);
2653     list->Add(fHistFracHadronsVsNMultiplicity);
2654     list->Add(fHistFracNeutronsVsNMultiplicity);
2655     list->Add(fHistFracKaonsVsNMultiplicity);
2656     list->Add(fHistFracSecondariesVsNMultiplicity);
2657     list->Add(fHistFracSignalVsNMatchedTracks);
2658     list->Add(fHistFracHadronsVsNMatchedTracks);
2659     list->Add(fHistFracNeutronsVsNMatchedTracks);
2660     list->Add(fHistFracKaonsVsNMatchedTracks);
2661     list->Add(fHistFracSecondariesVsNMatchedTracks);
2662     list->Add(fHistFracSignalVsNTotalTracks);
2663     list->Add(fHistFracHadronsVsNTotalTracks);
2664     list->Add(fHistFracNeutronsVsNTotalTracks);
2665     list->Add(fHistFracKaonsVsNTotalTracks);
2666     list->Add(fHistFracSecondariesVsNTotalTracks);
2667     list->Add(fHistRCorrVsPtVsCent);
2668
2669 }
2670
2671
2672 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
2673 {
2674     //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
2675     AliESDtrack *esdTrack = new AliESDtrack(part);
2676     // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2677
2678     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
2679
2680     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2681
2682     bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
2683     delete esdTrack;
2684
2685     return status;
2686 }
2687
2688 void AliAnalysisEtMonteCarlo::FillHistograms()
2689 {   // let base class fill its histograms, and us fill the local ones
2690     AliAnalysisEt::FillHistograms();
2691   if(!fIsMC) return;
2692     //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
2693
2694     fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
2695     fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
2696     fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
2697     fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
2698     fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
2699     fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
2700     fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
2701     fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
2702     fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
2703     fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
2704     fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
2705     fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
2706     fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
2707     fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
2708     fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
2709     fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
2710     fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
2711
2712     fHistEtRemovedGammas->Fill(fEtRemovedGammas, fClusterMult);
2713     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
2714     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
2715
2716 //     fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
2717 //                                             +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
2718 //                                             +fEtRemovedProtons.
2719 //                              fCentClass);
2720 //     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
2721 //
2722 //     fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
2723 //                                             +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
2724 //                                             +fEtNonRemovedProtons,
2725 //                              fCentClass);
2726 //     fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
2727
2728     fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fClusterMult);
2729     fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fClusterMult);
2730     fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fClusterMult);
2731     fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fClusterMult);
2732
2733     fHistMultRemovedCharged->Fill(fChargedRemoved, fClusterMult);
2734     fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fClusterMult);
2735     fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fClusterMult);
2736     fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fClusterMult);
2737
2738
2739     fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
2740     fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
2741     fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
2742     fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
2743     fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
2744     fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
2745     fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
2746     fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
2747     fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
2748     fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
2749     fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
2750     fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
2751     fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
2752     fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
2753     fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
2754     fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
2755
2756     fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
2757     fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
2758     fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
2759
2760     fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
2761                                             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
2762                                             +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
2763                                             +fMultNonRemovedProtons);
2764
2765     fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
2766                                             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
2767
2768     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
2769                                        fMultRemovedGammas);
2770
2771     fHistClusterMultvsNonRemovedCharged->Fill(fClusterMult,
2772             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
2773             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
2774             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
2775
2776     fHistClusterMultvsNonRemovedNeutral->Fill(fClusterMult,
2777             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
2778
2779     fHistClusterMultvsRemovedGamma->Fill(fClusterMult,
2780                                          fMultRemovedGammas);
2781
2782 }
2783
2784
2785
2786
2787 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
2788 { // print family tree
2789     TParticle *part = stack->Particle(partIdx);
2790 //     if(part->GetPdgCode() == fgK0SCode)
2791     {
2792       std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<<" is it from secondary interaction "<<fSelector->FromSecondaryInteraction(partIdx, *stack)<< std::endl;
2793         std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
2794         std::cout << "Energy: " << part->Energy() << std::endl;
2795         Float_t vtx = TMath::Sqrt( TMath::Power(part->Vx(),2) + TMath::Power(part->Vy(),2) + TMath::Power(part->Vz(),2) );
2796         Float_t vtxy = TMath::Sqrt( TMath::Power(part->Vx(),2) + TMath::Power(part->Vy(),2)  );
2797         std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() <<" |Vtx| "<<vtx <<" |Vtxy| "<<vtxy << std::endl;
2798     }
2799     return PrintMothers(partIdx, stack, 1);
2800 }
2801
2802 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
2803 { // print mothers
2804     char *tabs = new char[gen+1];
2805     for(Int_t i = 0; i < gen; ++i)
2806     {
2807         //std::cout << i << std::endl;
2808         tabs[i] = '\t';
2809     }
2810     tabs[gen] = '\0';
2811     Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2812     if(mothIdx < 0)
2813     {
2814       delete [] tabs;
2815       return 0;
2816     }
2817     TParticle *mother = stack->Particle(mothIdx);
2818 //     if(mother->GetPdgCode() == fgK0SCode)
2819     {
2820         //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
2821         std::cout << tabs << "Index: " << mothIdx << std::endl;
2822         std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx)<<" is it from secondary interaction "<<fSelector->FromSecondaryInteraction(partIdx, *stack) << std::endl;
2823         std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
2824         std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
2825         if(mother->GetFirstMother() >= 0)
2826         {
2827           std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
2828           if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
2829           std::cout << std::endl;
2830         }
2831         Float_t vtx = TMath::Sqrt( TMath::Power(mother->Vx(),2) + TMath::Power(mother->Vy(),2) + TMath::Power(mother->Vz(),2) );
2832         Float_t vtxy = TMath::Sqrt( TMath::Power(mother->Vx(),2) + TMath::Power(mother->Vy(),2)  );
2833         std::cout<<tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() <<"|Vtx| "<<vtx<<" |Vtxy| "<<vtxy << std::endl;
2834     }
2835     if(mother->GetPdgCode() == fgK0SCode)
2836     {
2837 //      std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
2838     }
2839 //  std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
2840 //  std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
2841 //  std::cout << "Energy: " << mother->Energy() << std::endl;
2842 //  std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
2843
2844     delete [] tabs;
2845     return PrintMothers(mothIdx, stack, gen+1) + 1;
2846 }
2847
2848 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
2849 { // get primary mother
2850     if(partIdx >= 0)
2851     {
2852         //return stack->GetPrimary(partIdx);
2853
2854         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2855         if(mothIdx < 0) return -1;
2856         TParticle *mother = stack->Particle(mothIdx);
2857         if(mother)
2858         {
2859             if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
2860             else return GetPrimMother(mothIdx, stack);
2861         }
2862         else
2863         {
2864             return -1;
2865         }
2866     }
2867     return -1;
2868 }
2869
2870 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
2871 { // get K0 in family
2872     if(partIdx >= 0)
2873     {
2874         if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
2875         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2876         if(mothIdx < 0) return -1;
2877         TParticle *mother = stack->Particle(mothIdx);
2878         if(mother)
2879         {
2880           if(mother->GetPdgCode() == fgK0SCode)
2881             {
2882                 return mothIdx;
2883             }
2884             return GetK0InFamily(mothIdx, stack);
2885         }
2886         else
2887         {
2888             return -1;
2889         }
2890     }
2891     return -1;
2892 }
2893 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTreeShort(Int_t partIdx, AliStack* stack)
2894 { // print family tree
2895     TParticle *part = stack->Particle(partIdx);
2896     if(part){
2897         Float_t vtxy = TMath::Sqrt( TMath::Power(part->Vx(),2) + TMath::Power(part->Vy(),2)  );
2898         cout<<part->GetName()<<"( is scondary "<<fSelector->FromSecondaryInteraction(partIdx, *stack)<<", vtx "<<vtxy <<", index "<<partIdx<<")";
2899         //cout<<"<-"<<part->GetName()<<"("<<fSelector->FromSecondaryInteraction(partIdx, *stack)<<","<<vtxy <<")";
2900     }
2901     else{return 0;}
2902     Int_t value = PrintMothersShort(partIdx, stack, 1);
2903     cout<<endl;
2904     return value;
2905 }
2906
2907 Int_t AliAnalysisEtMonteCarlo::PrintMothersShort(Int_t partIdx, AliStack* stack, Int_t gen)
2908 { // print mothers
2909     Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2910     if(mothIdx < 0)
2911     {
2912       return 0;
2913     }
2914     TParticle *mother = stack->Particle(mothIdx);
2915     if(mother){
2916       //Float_t vtx = TMath::Sqrt( TMath::Power(mother->Vx(),2) + TMath::Power(mother->Vy(),2) + TMath::Power(mother->Vz(),2) );
2917       Float_t vtxy = TMath::Sqrt( TMath::Power(mother->Vx(),2) + TMath::Power(mother->Vy(),2)  );
2918       cout<<"<-"<<mother->GetName()<<"( is scondary "<<fSelector->FromSecondaryInteraction(mothIdx, *stack)<<", vtx "<<vtxy <<", index "<<mothIdx<<")";
2919       //std::cout<<tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() <<"|Vtx| "<<vtx<<" |Vtxy| "<<vtxy << std::endl;
2920     }
2921     else{return 0;}
2922     return PrintMothersShort(mothIdx, stack, gen+1) + 1;
2923 }
2924
2925
2926 void AliAnalysisEtMonteCarlo::SetGeneratorMinMaxParticles(AliMCEvent *eventMC){
2927   // In case of access only to hijing particles in cocktail
2928   // get the min and max labels
2929   // TODO: Check when generator is not the first one ...
2930
2931   fNMCProducedMin = 0;
2932   fNMCProducedMax = 0;
2933
2934   AliGenEventHeader * eventHeader = eventMC->GenEventHeader();
2935
2936   AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(eventHeader);
2937
2938   if(!cocktail) return ;
2939
2940   TList *genHeaders = cocktail->GetHeaders();
2941
2942   Int_t nGenerators = genHeaders->GetEntries();
2943   //printf("N generators %d \n", nGenerators);
2944
2945   for(Int_t igen = 0; igen < nGenerators; igen++)
2946     {
2947       AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
2948       TString name = eventHeader2->GetName();
2949
2950       //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle());
2951
2952       fNMCProducedMin = fNMCProducedMax;
2953       fNMCProducedMax+= eventHeader2->NProduced();
2954
2955       if(name.Contains("Hijing",TString::kIgnoreCase)){
2956         //cout<<"Found HIJING event and set range "<<fNMCProducedMin<<"-"<<fNMCProducedMax<<endl;
2957         return ;
2958       }
2959     }
2960
2961 }
2962 AliGenEventHeader* AliAnalysisEtMonteCarlo::GetGenEventHeader(AliMCEvent *eventMC) const
2963 {
2964   // Return pointer to Generated event header
2965   // If requested and cocktail, search for the hijing generator
2966   AliGenEventHeader * eventHeader = eventMC->GenEventHeader();
2967   AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(eventHeader);
2968
2969   if(!cocktail) return 0x0 ;
2970
2971   TList *genHeaders = cocktail->GetHeaders();
2972
2973   Int_t nGenerators = genHeaders->GetEntries();
2974   //printf("N generators %d \n", nGenerators);
2975
2976   for(Int_t igen = 0; igen < nGenerators; igen++)
2977     {
2978       AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
2979       TString name = eventHeader2->GetName();
2980
2981       //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle());
2982
2983       if(name.Contains("Hijing",TString::kIgnoreCase)) return eventHeader2 ;
2984     }
2985
2986   return 0x0;
2987
2988 }
2989 Bool_t AliAnalysisEtMonteCarlo::IsHIJINGLabel(Int_t label,AliMCEvent *eventMC,AliStack *stack)
2990 {
2991
2992   // Find if cluster/track was generated by HIJING
2993
2994   AliGenHijingEventHeader*  hijingHeader =  dynamic_cast<AliGenHijingEventHeader *> (GetGenEventHeader(eventMC));
2995
2996   //printf("header %p, label %d\n",hijingHeader,label);
2997
2998   if(!hijingHeader || label < 0 ) return kFALSE;
2999
3000
3001   //printf("pass a), N produced %d\n",nproduced);
3002
3003   if(label >= fNMCProducedMin && label < fNMCProducedMax)
3004   {
3005     //printf(" accept!, label is smaller than produced, N %d\n",nproduced);
3006
3007     return kTRUE;
3008   }
3009
3010   if(!stack) return kFALSE;
3011
3012   Int_t nprimaries = stack->GetNtrack();
3013
3014   if(label > nprimaries) return kFALSE;
3015
3016   TParticle * mom = stack->Particle(label);
3017
3018   Int_t iMom = label;
3019   Int_t iParent = mom->GetFirstMother();
3020   while(iParent!=-1){
3021     if(iParent >= fNMCProducedMin && iParent < fNMCProducedMax){
3022       return kTRUE;
3023     }
3024
3025     iMom = iParent;
3026     mom = stack->Particle(iMom);
3027     iParent = mom->GetFirstMother();
3028   }
3029
3030   return kFALSE ;
3031
3032 }
3033
3034
3035
3036