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