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