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