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                                                   ,fIsMC(kTRUE)
  41     ,fImpactParameter(0)
  42     ,fNcoll(0)
  43     ,fNpart(0)
  44     ,fPrimaryTree(0)
  45     ,fTotEtWithSecondaryRemoved(0)
  46     ,fTotEtSecondaryFromEmEtPrimary(0)
  47     ,fTotEtSecondary(0)
  48     ,fPrimaryCode(0)
  49     ,fPrimaryCharge(0)
  50     ,fPrimaryE(0)
  51     ,fPrimaryEt(0)
  52     ,fPrimaryPx(0)
  53     ,fPrimaryPy(0)
  54     ,fPrimaryPz(0)
  55     ,fPrimaryVx(0)
  56     ,fPrimaryVy(0)
  57     ,fPrimaryVz(0)
  58     ,fPrimaryAccepted(0)
  59     ,fPrimaryMatched(0)
  60     ,fDepositedCode(0)
  61     ,fDepositedE(0)
  62     ,fDepositedEt(0)
  63     ,fDepositedCharge(0)
  64     ,fDepositedVx(0)
  65     ,fDepositedVy(0)
  66     ,fDepositedVz(0)
  67     ,fSecondary(kFALSE)
  68     ,fReconstructedE(0)
  69     ,fReconstructedEt(0)
  70     ,fTotPx(0)
  71     ,fTotPy(0)
  72     ,fTotPz(0)
  73     ,fClusterMult(0)
  74     ,fHistDecayVertexNonRemovedCharged(0)
  75     ,fHistDecayVertexRemovedCharged(0)
  76     ,fHistDecayVertexNonRemovedNeutral(0)
  77     ,fHistDecayVertexRemovedNeutral(0)
  78
  79     ,fHistRemovedOrNot(0)
  80     ,fHistEtNonRemovedProtons(0)
  81     ,fHistEtNonRemovedAntiProtons(0)
  82     ,fHistEtNonRemovedPiPlus(0)
  83     ,fHistEtNonRemovedPiMinus(0)
  84     ,fHistEtNonRemovedKaonPlus(0)
  85     ,fHistEtNonRemovedKaonMinus(0)
  86     ,fHistEtNonRemovedK0s(0)
  87     ,fHistEtNonRemovedK0L(0)
  88     ,fHistEtNonRemovedLambdas(0)
  89     ,fHistEtNonRemovedElectrons(0)
  90     ,fHistEtNonRemovedPositrons(0)
  91     ,fHistEtNonRemovedMuPlus(0)
  92     ,fHistEtNonRemovedMuMinus(0)
  93     ,fHistEtNonRemovedNeutrons(0)
  94     ,fHistEtNonRemovedAntiNeutrons(0)
  95     ,fHistEtNonRemovedGammas(0)
  96     ,fHistEtNonRemovedGammasFromPi0(0)
  97     ,fHistEtRemovedGammas(0)
  98     ,fHistEtRemovedNeutrons(0)
  99     ,fHistEtRemovedAntiNeutrons(0)
 100     ,fHistEtRemovedCharged(0)
 101     ,fHistEtRemovedNeutrals(0)
 102     ,fHistEtNonRemovedCharged(0)
 103     ,fHistEtNonRemovedNeutrals(0)
 104     ,fHistMultNonRemovedProtons(0)
 105     ,fHistMultNonRemovedAntiProtons(0)
 106     ,fHistMultNonRemovedPiPlus(0)
 107     ,fHistMultNonRemovedPiMinus(0)
 108     ,fHistMultNonRemovedKaonPlus(0)
 109     ,fHistMultNonRemovedKaonMinus(0)
 110     ,fHistMultNonRemovedK0s(0)
 111     ,fHistMultNonRemovedK0L(0)
 112     ,fHistMultNonRemovedLambdas(0)
 113     ,fHistMultNonRemovedElectrons(0)
 114     ,fHistMultNonRemovedPositrons(0)
 115     ,fHistMultNonRemovedMuPlus(0)
 116     ,fHistMultNonRemovedMuMinus(0)
 117     ,fHistMultNonRemovedNeutrons(0)
 118     ,fHistMultNonRemovedAntiNeutrons(0)
 119     ,fHistMultNonRemovedGammas(0)
 120     ,fHistMultRemovedGammas(0)
 121     ,fHistMultRemovedNeutrons(0)
 122     ,fHistMultRemovedAntiNeutrons(0)
 123     ,fHistMultRemovedCharged(0)
 124     ,fHistMultRemovedNeutrals(0)
 125     ,fHistMultNonRemovedCharged(0)
 126     ,fHistMultNonRemovedNeutrals(0)
 127     ,fHistTrackMultvsNonRemovedCharged(0)
 128     ,fHistTrackMultvsNonRemovedNeutral(0)
 129     ,fHistTrackMultvsRemovedGamma(0)
 130     ,fHistClusterMultvsNonRemovedCharged(0)
 131     ,fHistClusterMultvsNonRemovedNeutral(0)
 132     ,fHistClusterMultvsRemovedGamma(0)
 133     ,fHistMultvsNonRemovedChargedE(0)
 134     ,fHistMultvsNonRemovedNeutralE(0)
 135     ,fHistMultvsRemovedGammaE(0)
 136                                                   ,fCalcForKaonCorrection(kFALSE)
 137                                                   ,fHistK0EDepositsVsPtInAcceptance(0)
 138                                                   ,fHistK0EGammaVsPtInAcceptance(0)
 139                                                   ,fHistK0EDepositsVsPtOutOfAcceptance(0)
 140                                                   ,fHistK0EGammaVsPtOutOfAcceptance(0)
 141                                                   ,fHistSimKaonsInAcceptance(0)
 142                                                   ,fHistSimK0SInAcceptance(0)
 143                                                   ,fHistSimKPlusInAcceptance(0)
 144                                                   ,fHistSimKMinusInAcceptance(0)
 145                                                   ,fHistSimK0LInAcceptance(0)
 146                                                   ,fHistSimKaonsOutOfAcceptance(0)
 147                                                   ,fHistSimKaonsInAcceptanceWithDepositsPrimaries(0)
 148                                                   ,fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries(0)
 149                                                   ,fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries(0)
 150     ,fEtNonRemovedProtons(0)
 151     ,fEtNonRemovedAntiProtons(0)
 152     ,fEtNonRemovedPiPlus(0)
 153     ,fEtNonRemovedPiMinus(0)
 154     ,fEtNonRemovedKaonPlus(0)
 155     ,fEtNonRemovedKaonMinus(0)
 156     ,fEtNonRemovedK0S(0)
 157     ,fEtNonRemovedK0L(0)
 158     ,fEtNonRemovedLambdas(0)
 159     ,fEtNonRemovedElectrons(0)
 160     ,fEtNonRemovedPositrons(0)
 161     ,fEtNonRemovedMuMinus(0)
 162     ,fEtNonRemovedMuPlus(0)
 163     ,fEtNonRemovedGammas(0)
 164     ,fEtNonRemovedGammasFromPi0(0)
 165     ,fEtNonRemovedNeutrons(0)
 166     ,fEtNonRemovedAntiNeutrons(0)
 167     ,fEtRemovedProtons(0)
 168     ,fEtRemovedAntiProtons(0)
 169     ,fEtRemovedPiPlus(0)
 170     ,fEtRemovedPiMinus(0)
 171     ,fEtRemovedKaonPlus(0)
 172     ,fEtRemovedKaonMinus(0)
 173     ,fEtRemovedK0s(0)
 174     ,fEtRemovedK0L(0)
 175     ,fEtRemovedLambdas(0)
 176     ,fEtRemovedElectrons(0)
 177     ,fEtRemovedPositrons(0)
 178     ,fEtRemovedMuMinus(0)
 179     ,fEtRemovedMuPlus(0)
 180     ,fEtRemovedGammasFromPi0(0)
 181     ,fEtRemovedGammas(0)
 182     ,fEtRemovedNeutrons(0)
 183     ,fEtRemovedAntiNeutrons(0)
 184     ,fMultNonRemovedProtons(0)
 185     ,fMultNonRemovedAntiProtons(0)
 186     ,fMultNonRemovedPiPlus(0)
 187     ,fMultNonRemovedPiMinus(0)
 188     ,fMultNonRemovedKaonPlus(0)
 189     ,fMultNonRemovedKaonMinus(0)
 190     ,fMultNonRemovedK0s(0)
 191     ,fMultNonRemovedK0L(0)
 192     ,fMultNonRemovedLambdas(0)
 193     ,fMultNonRemovedElectrons(0)
 194     ,fMultNonRemovedPositrons(0)
 195     ,fMultNonRemovedMuMinus(0)
 196     ,fMultNonRemovedMuPlus(0)
 197     ,fMultNonRemovedGammas(0)
 198     ,fMultNonRemovedNeutrons(0)
 199     ,fMultNonRemovedAntiNeutrons(0)
 200     ,fMultRemovedProtons(0)
 201     ,fMultRemovedAntiProtons(0)
 202     ,fMultRemovedPiPlus(0)
 203     ,fMultRemovedPiMinus(0)
 204     ,fMultRemovedKaonPlus(0)
 205     ,fMultRemovedKaonMinus(0)
 206     ,fMultRemovedK0s(0)
 207     ,fMultRemovedK0L(0)
 208     ,fMultRemovedLambdas(0)
 209     ,fMultRemovedElectrons(0)
 210     ,fMultRemovedPositrons(0)
 211     ,fMultRemovedMuMinus(0)
 212     ,fMultRemovedMuPlus(0)
 213     ,fMultRemovedGammas(0)
 214     ,fMultRemovedNeutrons(0)
 215     ,fMultRemovedAntiNeutrons(0)
 216     ,fTrackMultInAcc(0)
 217     ,fHistDxDzNonRemovedCharged(0)
 218     ,fHistDxDzRemovedCharged(0)
 219     ,fHistDxDzNonRemovedNeutral(0)
 220     ,fHistDxDzRemovedNeutral(0)
 221     ,fHistPiPlusMult(0)
 222     ,fHistPiMinusMult(0)
 223     ,fHistPiZeroMult(0)
 224     ,fHistPiPlusMultAcc(0)
 225     ,fHistPiMinusMultAcc(0)
 226     ,fHistPiZeroMultAcc(0)
 227 //     ,fPiPlusMult(0)
 228 //     ,fPiMinusMult(0)
 229     ,fPiZeroMult(0)
 230     ,fPiPlusMultAcc(0)
 231     ,fPiMinusMultAcc(0)
 232     ,fPiZeroMultAcc(0)
 233     ,fNeutralRemoved(0)
 234     ,fChargedRemoved(0)
 235     ,fChargedNotRemoved(0)
 236     ,fNeutralNotRemoved(0)
 237     ,fGammaRemoved(0)
 238     ,fSecondaryNotRemoved(0)
 239     ,fEnergyNeutralRemoved(0)
 240     ,fEnergyChargedRemoved(0)
 241     ,fEnergyChargedNotRemoved(0)
 242     ,fEnergyNeutralNotRemoved(0)
 243     ,fEnergyGammaRemoved(0)
 244     ,fNClusters(0)
 245     ,fTotNeutralEtAfterMinEnergyCut(0)
 246                                                   ,fCalcTrackMatchVsMult(kFALSE)
 247                                                   ,fHistGammasFound(0)
 248                                                   ,fHistGammasGenerated(0)
 249                                                   ,fHistGammasFoundCent(0)
 250                                                   ,fHistGammasFoundOutOfAccCent(0)
 251                                                   ,fHistGammasGeneratedCent(0)
 252                                                   ,fHistChargedTracksCut(0)
 253                                                   ,fHistChargedTracksAccepted(0)
 254                                                   ,fHistGammasCut(0)
 255                                                   ,fHistGammasAccepted(0)
 256                                                   ,fHistChargedTracksCutMult(0)
 257                                                   ,fHistChargedTracksAcceptedMult(0)
 258                                                   ,fHistChargedTracksAcceptedLowPtCent(0)
 259                                                   ,fHistChargedTracksAcceptedLowPtCent500MeV(0)
 260                                                   ,fHistChargedTracksAcceptedLowPtCentNoAntiProtons(0)
 261                                                   ,fHistChargedTracksAcceptedLowPtCentAntiProtons(0)
 262                                                   ,fHistGammasCutMult(0)
 263                                                   ,fHistGammasAcceptedMult(0)
 264                                                   ,fHistBadTrackMatches(0)
 265                                                   ,fHistMatchedTracksEvspTBkgd(0)
 266                                                   ,fHistMatchedTracksEvspTSignal(0)
 267                                                   ,fHistMatchedTracksEvspTBkgdPeripheral(0)
 268                                                   ,fHistMatchedTracksEvspTSignalPeripheral(0)
 269                                                   ,fHistMatchedTracksEvspTBkgdvsCent(0)
 270                                                   ,fHistMatchedTracksEvspTSignalvsCent(0)
 271                                                   ,fHistMatchedTracksEvspTBkgdvsCentEffCorr(0)
 272                                                   ,fHistMatchedTracksEvspTSignalvsCentEffCorr(0)
 273
 274                                                   ,fHistChargedTracksCutPeripheral(0)
 275                                                   ,fHistChargedTracksAcceptedPeripheral(0)
 276                                                   ,fHistGammasCutPeripheral(0)
 277                                                   ,fHistGammasAcceptedPeripheral(0)
 278                                                   ,fHistBadTrackMatchesdPhidEta(0)
 279                                                   ,fHistGoodTrackMatchesdPhidEta(0)
 280                                                   ,fHistHadronDepositsAll(0)
 281                                                   ,fHistHadronDepositsReco(0)
 282                                                   ,fHistHadronDepositsAllCent(0)
 283                                                   ,fHistHadronDepositsAllCent500MeV(0)
 284                                                   ,fHistHadronDepositsRecoCent(0)
 285                                                   ,fHistHadronsAllCent(0)
 286                                                   ,fHistMultChVsSignalVsMult(0)
 287                                                   ,fHistNeutralRemovedSecondaryEtVsCent(0)
 288                                                   ,fHistChargedRemovedSecondaryEtVsCent(0)
 289                                                   ,fHistNeutralNotRemovedSecondaryEtVsCent(0)
 290                                                   ,fHistChargedNotRemovedSecondaryEtVsCent(0)
 291                                                   ,fHistNeutralRemovedSecondaryNumVsNCluster(0)
 292                                                   ,fHistChargedRemovedSecondaryNumVsNCluster(0)
 293                                                   ,fHistNeutralNotRemovedSecondaryNumVsNCluster(0)
 294                                                   ,fHistChargedNotRemovedSecondaryNumVsNCluster(0)
 295                                                   ,fHistNeutralRemovedSecondaryNumVsCent(0)
 296                                                   ,fHistChargedRemovedSecondaryNumVsCent(0)
 297                                                   ,fHistNeutralNotRemovedSecondaryNumVsCent(0)
 298                                                   ,fHistChargedNotRemovedSecondaryNumVsCent(0)
 299                                                   ,fHistNeutronsEtVsCent(0)
 300                                                   ,fHistNeutronsNumVsCent(0)
 301                                                   ,fHistNotNeutronsNumVsCent(0)
 302                                                   ,fHistPiKPDepositedVsNch(0)
 303                                                   ,fHistPiKPNotTrackMatchedDepositedVsNch(0)
 304                                                   ,fHistNeutronsDepositedVsNch(0)
 305                                                   ,fHistAntiNeutronsDepositedVsNch(0)
 306                                                   ,fHistProtonsDepositedVsNch(0)
 307                                                   ,fHistAntiProtonsDepositedVsNch(0)
 308                                                   ,fHistProtonsNotTrackMatchedDepositedVsNch(0)
 309                                                   ,fHistAntiProtonsNotTrackMatchedDepositedVsNch(0)
 310                                                   ,fHistSecondariesVsNch(0)
 311                                                   ,fHistSecondariesVsNcl(0)
 312                                                   ,fHistSecondariesEffCorrVsNch(0)
 313                                                   ,fHistSecondariesEffCorrVsNcl(0)
 314                                                   ,fHistSecondariesOutOfAccEffCorrVsNch(0)
 315                                                   ,fHistSecondariesDetectorCoverEffCorrVsNch(0)
 316                                                   ,fHistCentVsNchVsNcl(0)
 317                                                 //,fHistSecondaryPositionInDetector(0)
 318                                                   ,fClusterPositionWeird(0)
 319                                                 //,fHistSecondaryPositionInDetectorMultiple(0)
 320                                                   ,fSecondaryClusterEnergy(0)
 321 {
 322 }
 323
 324 // dtor
 325 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
 326 {   //Destructor
 327
 328   if(fPrimaryTree){
 329     fPrimaryTree->Clear();
 330     delete fPrimaryTree;
 331   }
 332     delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
 333     delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
 334     delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
 335     delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
 336
 337     delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
 338
 339     delete fHistEtNonRemovedProtons; // enter comment here
 340     delete fHistEtNonRemovedAntiProtons; // enter comment here
 341     delete fHistEtNonRemovedPiPlus; // enter comment here
 342     delete fHistEtNonRemovedPiMinus; // enter comment here
 343     delete fHistEtNonRemovedKaonPlus; // enter comment here
 344     delete fHistEtNonRemovedKaonMinus; // enter comment here
 345     delete fHistEtNonRemovedK0s; // enter comment here
 346     delete fHistEtNonRemovedK0L; // enter comment here
 347     delete fHistEtNonRemovedLambdas; // enter comment here
 348     delete fHistEtNonRemovedElectrons; // enter comment here
 349     delete fHistEtNonRemovedPositrons; // enter comment here
 350     delete fHistEtNonRemovedMuPlus; // enter comment here
 351     delete fHistEtNonRemovedMuMinus; // enter comment here
 352     delete fHistEtNonRemovedNeutrons; // enter comment here
 353     delete fHistEtNonRemovedAntiNeutrons; // enter comment here
 354     delete fHistEtNonRemovedGammas; // enter comment here
 355     delete fHistEtNonRemovedGammasFromPi0; // enter comment here
 356
 357     delete fHistEtRemovedGammas; // enter comment here
 358     delete fHistEtRemovedNeutrons; // enter comment here
 359     delete fHistEtRemovedAntiNeutrons; // enter comment here
 360
 361
 362     delete fHistMultNonRemovedProtons; // enter comment here
 363     delete fHistMultNonRemovedAntiProtons; // enter comment here
 364     delete fHistMultNonRemovedPiPlus; // enter comment here
 365     delete fHistMultNonRemovedPiMinus; // enter comment here
 366     delete fHistMultNonRemovedKaonPlus; // enter comment here
 367     delete fHistMultNonRemovedKaonMinus; // enter comment here
 368     delete fHistMultNonRemovedK0s; // enter comment here
 369     delete fHistMultNonRemovedK0L; // enter comment here
 370     delete fHistMultNonRemovedLambdas; // enter comment here
 371     delete fHistMultNonRemovedElectrons; // enter comment here
 372     delete fHistMultNonRemovedPositrons; // enter comment here
 373     delete fHistMultNonRemovedMuPlus; // enter comment here
 374     delete fHistMultNonRemovedMuMinus; // enter comment here
 375     delete fHistMultNonRemovedNeutrons; // enter comment here
 376     delete fHistMultNonRemovedAntiNeutrons; // enter comment here
 377     delete fHistMultNonRemovedGammas; // enter comment here
 378
 379     delete fHistMultRemovedGammas; // enter comment here
 380     delete fHistMultRemovedNeutrons; // enter comment here
 381     delete fHistMultRemovedAntiNeutrons; // enter comment here
 382
 383     delete fHistTrackMultvsNonRemovedCharged; // enter comment here
 384     delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
 385     delete fHistTrackMultvsRemovedGamma; // enter comment here
 386
 387     delete fHistClusterMultvsNonRemovedCharged; // enter comment here
 388     delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
 389     delete fHistClusterMultvsRemovedGamma; // enter comment here
 390
 391     delete fHistMultvsNonRemovedChargedE; // enter comment here
 392     delete fHistMultvsNonRemovedNeutralE; // enter comment here
 393     delete fHistMultvsRemovedGammaE; // enter comment here
 394     delete fHistK0EDepositsVsPtInAcceptance;//enter comment here
 395     delete fHistK0EGammaVsPtInAcceptance;//enter comment here
 396     delete fHistK0EDepositsVsPtOutOfAcceptance;
 397     delete fHistK0EGammaVsPtOutOfAcceptance;
 398
 399     delete fHistSimKaonsInAcceptance;
 400     delete fHistSimK0SInAcceptance;
 401     delete fHistSimKPlusInAcceptance;
 402     delete fHistSimKMinusInAcceptance;
 403     delete fHistSimK0LInAcceptance;
 404     delete fHistSimKaonsOutOfAcceptance;
 405     delete fHistSimKaonsInAcceptanceWithDepositsPrimaries;
 406     delete fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries;
 407     delete fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries;
 408
 409     delete fHistDxDzNonRemovedCharged; // enter comment here
 410     delete fHistDxDzRemovedCharged; // enter comment here
 411     delete fHistDxDzNonRemovedNeutral; // enter comment here
 412     delete fHistDxDzRemovedNeutral; // enter comment here
 413
 414     delete fHistPiPlusMult; // enter comment here
 415     delete fHistPiMinusMult; // enter comment here
 416     delete fHistPiZeroMult; // enter comment here
 417
 418     delete fHistPiPlusMultAcc; // enter comment here
 419     delete fHistPiMinusMultAcc; // enter comment here
 420     delete fHistPiZeroMultAcc; // enter comment here
 421     delete fHistGammasFound; // enter comment here
 422     delete fHistGammasGenerated; // enter comment here
 423     delete fHistGammasFoundOutOfAccCent; // enter comment here
 424     delete fHistGammasFoundCent; // enter comment here
 425     delete fHistGammasGeneratedCent; // enter comment here
 426     delete fHistChargedTracksCut;
 427     delete fHistChargedTracksAccepted;
 428     delete fHistGammasCut;
 429     delete fHistGammasAccepted;
 430     delete fHistChargedTracksCutMult;
 431     delete fHistChargedTracksAcceptedMult;
 432     delete fHistChargedTracksAcceptedLowPtCent;
 433     delete fHistChargedTracksAcceptedLowPtCent500MeV;
 434     delete fHistChargedTracksAcceptedLowPtCentNoAntiProtons;
 435     delete fHistChargedTracksAcceptedLowPtCentAntiProtons;
 436     delete fHistGammasCutMult;
 437     delete fHistGammasAcceptedMult;
 438     delete fHistBadTrackMatches;
 439     delete fHistMatchedTracksEvspTBkgd;
 440     delete fHistMatchedTracksEvspTSignal;
 441     delete fHistMatchedTracksEvspTBkgdPeripheral;
 442     delete fHistMatchedTracksEvspTSignalPeripheral;
 443     delete fHistMatchedTracksEvspTBkgdvsCent;
 444     delete fHistMatchedTracksEvspTSignalvsCent;
 445     delete fHistMatchedTracksEvspTBkgdvsCentEffCorr;
 446     delete fHistMatchedTracksEvspTSignalvsCentEffCorr;
 447     delete fHistChargedTracksCutPeripheral;
 448     delete fHistChargedTracksAcceptedPeripheral;
 449     delete fHistGammasCutPeripheral;
 450     delete fHistGammasAcceptedPeripheral;
 451     delete fHistBadTrackMatchesdPhidEta;
 452     delete fHistGoodTrackMatchesdPhidEta;
 453     delete fHistHadronDepositsAll;
 454     delete fHistHadronDepositsReco;
 455     delete fHistHadronDepositsAllCent;
 456     delete fHistHadronDepositsAllCent500MeV;
 457     delete fHistHadronDepositsRecoCent;
 458     delete fHistHadronsAllCent;
 459     delete fHistMultChVsSignalVsMult;
 460     delete fHistNeutralRemovedSecondaryEtVsCent;
 461     delete fHistChargedRemovedSecondaryEtVsCent;
 462     delete fHistNeutralNotRemovedSecondaryEtVsCent;
 463     delete fHistChargedNotRemovedSecondaryEtVsCent;
 464     delete fHistNeutralRemovedSecondaryNumVsNCluster;
 465     delete fHistChargedRemovedSecondaryNumVsNCluster;
 466     delete fHistNeutralNotRemovedSecondaryNumVsNCluster;
 467     delete fHistChargedNotRemovedSecondaryNumVsNCluster;
 468     delete fHistNeutralRemovedSecondaryNumVsCent;
 469     delete fHistChargedRemovedSecondaryNumVsCent;
 470     delete fHistNeutralNotRemovedSecondaryNumVsCent;
 471     delete fHistChargedNotRemovedSecondaryNumVsCent;
 472     delete fHistNeutronsEtVsCent;
 473     delete fHistNeutronsNumVsCent;
 474     delete fHistNotNeutronsNumVsCent;
 475     delete fHistPiKPDepositedVsNch;
 476     delete fHistPiKPNotTrackMatchedDepositedVsNch;
 477     delete fHistNeutronsDepositedVsNch;
 478     delete fHistAntiNeutronsDepositedVsNch;
 479     delete fHistProtonsDepositedVsNch;
 480     delete fHistAntiProtonsDepositedVsNch;
 481     delete fHistProtonsNotTrackMatchedDepositedVsNch;
 482     delete fHistAntiProtonsNotTrackMatchedDepositedVsNch;
 483     delete fHistSecondariesVsNch;
 484     delete fHistSecondariesVsNcl;
 485     delete fHistSecondariesEffCorrVsNch;
 486     delete fHistSecondariesEffCorrVsNcl;
 487     delete fHistSecondariesOutOfAccEffCorrVsNch;
 488     delete fHistSecondariesDetectorCoverEffCorrVsNch;
 489     delete fHistCentVsNchVsNcl;
 490     //delete fHistSecondaryPositionInDetector;
 491     delete fClusterPositionWeird;
 492     //delete fHistSecondaryPositionInDetectorMultiple;
 493     delete fSecondaryClusterEnergy;
 494 }
 495
 496 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 497 {   // analyse MC event
 498   if(!fIsMC) return 0;
 499     ResetEventValues();
 500
 501
 502     fPiPlusMult = 0;
 503     fPiMinusMult = 0;
 504     fPiZeroMult = 0;
 505     if (fCentrality)
 506     {
 507         fCentClass = fCentrality->GetCentralityClass5(fCentralityMethod);
 508
 509     }
 510
 511     // Get us an mc event
 512     if (!ev) {
 513         AliFatal("ERROR: Event does not exist");
 514         return 0;
 515     }
 516     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
 517     if (!event) {
 518         AliFatal("ERROR: MC Event does not exist");
 519         return 0;
 520     }
 521
 522     Double_t protonMass =fgProtonMass;
 523
 524     // Hijing header
 525     AliGenEventHeader* genHeader = event->GenEventHeader();
 526     if (!genHeader) {
 527         Printf("ERROR: Event generation header does not exist");
 528         return 0;
 529     }
 530     AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 531     if (hijingGenHeader) {
 532         fImpactParameter = hijingGenHeader->ImpactParameter();
 533         fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
 534         fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
 535     }
 536
 537     // Let's play with the stack!
 538     AliStack *stack = event->Stack();
 539
 540     Int_t nPrim = stack->GetNtrack();
 541
 542     Int_t partCount = 0;
 543     for (Int_t iPart = 0; iPart < nPrim; iPart++)
 544     {
 545
 546         TParticle *part = stack->Particle(iPart);
 547
 548
 549
 550         if (!part)
 551         {
 552             Printf("ERROR: Could not get particle %d", iPart);
 553             continue;
 554         }
 555         TParticlePDG *pdg = part->GetPDG(0);
 556
 557         if (!pdg)
 558         {
 559             Printf("ERROR: Could not get particle PDG %d", iPart);
 560             continue;
 561         }
 562
 563         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
 564         Int_t code = pdg->PdgCode();
 565
 566         if(stack->IsPhysicalPrimary(iPart))
 567         {
 568           fTotPx += part->Px();
 569           fTotPy += part->Py();
 570           fTotPz += part->Pz();
 571         }
 572
 573         // Check for reasonable (for now neutral and singly charged) charge on the particle
 574         if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
 575         {
 576
 577             fMultiplicity++;
 578 //             PrintFamilyTree(iPart, stack);
 579 //
 580             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 581             {
 582                 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
 583                 // calculate E_T
 584                 if (
 585                     TMath::Abs(code) == fgProtonCode ||
 586                     TMath::Abs(code) == fgNeutronCode ||
 587                     TMath::Abs(code) == fgLambdaCode ||
 588                     TMath::Abs(code) == fgXiCode ||
 589                     TMath::Abs(code) == fgXi0Code ||
 590                     TMath::Abs(code) == fgOmegaCode
 591                 )
 592                 {
 593                     if (code > 0) {
 594                         particleMassPart = - protonMass;
 595                     }
 596                     if (code < 0) {
 597                         particleMassPart = protonMass;
 598                     }
 599                 }
 600                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 601
 602
 603 //                 // Fill up total E_T counters for each particle species
 604 //                 if (code == fgProtonCode || code == fgAntiProtonCode)
 605 //                 {
 606 //                 }
 607 //                 if (code == fgPiPlusCode || code == fgPiMinusCode)
 608 //                 {
 609 //                     if (code == fgPiPlusCode)
 610 //                     {
 611 //                     }
 612 //                     else
 613 //                     {
 614 //                     }
 615 //                 }
 616 //                 if (code == fgGammaCode)
 617 //                 {
 618 //                 }
 619 //                 if (code == fgKPlusCode)
 620 //                 {
 621 //                 }
 622 //                 if(code == fgKMinusCode)
 623 //                 {
 624 //                 }
 625 //                 if (code == fgMuPlusCode || code == fgMuMinusCode)
 626 //                 {
 627 //                 }
 628 //                 if (code == fgEPlusCode || code == fgEMinusCode)
 629 //                 {
 630 //                 }
 631 //                 // some neutrals also
 632 //                 if (code == fgNeutronCode)
 633 //                 {
 634 //                 }
 635 //                 if (code == fgAntiNeutronCode)
 636 //                 {
 637 //                 }
 638 //                 if (code == fgGammaCode)
 639 //                 {
 640 //                 }
 641
 642                 // Neutral particles
 643                 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
 644
 645                 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
 646                 {
 647                 //  PrintFamilyTree(iPart,stack);
 648                     //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
 649                     //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
 650
 651                     // inside EMCal acceptance
 652
 653                     //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
 654
 655                     if(fSelector->CutGeometricalAcceptance(*part) )
 656                     {
 657                         fNeutralMultiplicity++;
 658                         fTotNeutralEt += et;
 659                         if(fSelector->PassMinEnergyCut(*part))
 660                         {
 661                             fTotNeutralEtAfterMinEnergyCut += et;
 662                         }
 663                         if (part->Energy() > 0.05) partCount++;
 664                     }
 665                 }
 666                 //Charged particles
 667                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
 668                 {
 669
 670                     // inside EMCal acceptance
 671                     if (fSelector->CutGeometricalAcceptance(*part))
 672                     {
 673
 674                         fChargedMultiplicity++;
 675
 676                         fTotChargedEt += et;
 677
 678                         if (code == fgProtonCode || code == fgAntiProtonCode)
 679                         {
 680                         }
 681                         if (code == fgPiPlusCode || code == fgPiMinusCode)
 682                         {
 683                         }
 684                         if (code == fgKPlusCode || code == fgKMinusCode)
 685                         {
 686                         }
 687                         if (code == fgMuPlusCode || code == fgMuMinusCode)
 688                         {
 689                         }
 690
 691                         if (code == fgEPlusCode || code == fgEMinusCode)
 692                         {
 693                             fTotNeutralEt += et; // calling electrons neutral
 694                             fTotChargedEt -= et;
 695                         }
 696                     } // inside EMCal acceptance
 697
 698                     if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 699                     {
 700                         if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 701                         else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 702                     }
 703                 }
 704             }
 705
 706         }
 707     }
 708    // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
 709     fTotEt = fTotChargedEt + fTotNeutralEt;
 710     //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
 711     //std::cout << "Event done! # of particles: " << partCount << std::endl;
 712     return 0;
 713 }
 714 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 715 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 716 {   // analyse MC and real event info
 717     //if(!mcEvent || !realEvent){
 718   if(!fIsMC) return 0;
 719     if (!ev || !ev2) {
 720         AliFatal("ERROR: Event does not exist");
 721         return 0;
 722     }
 723     AliAnalysisEt::AnalyseEvent(ev);
 724     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 725     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
 726     if (!mcEvent || !realEvent) {
 727         AliFatal("ERROR: mcEvent or realEvent does not exist");
 728
 729     }
 730
 731     std::vector<Int_t> foundGammas;
 732
 733     fSelector->SetEvent(realEvent);
 734
 735     AnalyseEvent(ev);
 736
 737     AliStack *stack = mcEvent->Stack();
 738
 739     // get all detector clusters
 740     //  TRefArray* caloClusters = new TRefArray();
 741
 742 //    if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
 743     //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
 744     //else {
 745     //AliFatal("Detector ID has not been specified");
 746     //return -1;
 747 //    }
 748
 749 //Note that this only returns clusters for the selected detector.  fSelector actually calls the right GetClusters... for the detector
 750 //It does not apply any cuts on these clusters
 751     TRefArray *caloClusters = fSelector->GetClusters();
 752
 753     Int_t nCluster = caloClusters->GetEntries();
 754     fClusterMult = nCluster;
 755     fNClusters = 0;
 756     Int_t fClusterMultChargedTracks = 0;
 757     Int_t fClusterMultGammas = 0;
 758     Int_t nChargedSecondariesRemoved = 0;
 759     Int_t nChargedSecondariesNotRemoved = 0;
 760     Int_t nNeutralSecondariesRemoved = 0;
 761     Int_t nNeutralSecondariesNotRemoved = 0;
 762     Int_t nNeutrons = 0;
 763     Int_t nNotNeutrons = 0;
 764     //All efficiency corrected except etSecondaries -->
 765     Float_t etPiKPDeposited = 0.0;//
 766     Float_t etPiKPDepositedNotTrackMatched = 0.0;//
 767     Float_t etProtonDepositedNotTrackMatched = 0.0;//
 768     Float_t etAntiProtonDepositedNotTrackMatched = 0.0;//
 769     //Float_t etPIDProtonDepositedNotTrackMatched = 0.0;//Still has to be filled!
 770     //1Float_t etPIDAntiProtonDepositedNotTrackMatched = 0.0;//Still has to be filled
 771     Float_t etProtonDeposited = 0.0;//
 772     Float_t etAntiProtonDeposited = 0.0;//
 773     Float_t etNeutronDeposited = 0.0;
 774     Float_t etAntiNeutronDeposited = 0.0;
 775     Float_t etSecondaries = 0.0;
 776     Float_t etSecondariesEffCorr = 0.0;
 777     Float_t etSecondariesOutOfAccEffCorr = 0.0;
 778     Float_t etSecondariesDetectorCoverEffCorr = 0.0;
 779     Float_t multiplicity = fEsdtrackCutsTPC->GetReferenceMultiplicity(realEvent,kTRUE);
 780     // loop the clusters
 781     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 782     {
 783         Int_t cf = 0;
 784         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 785         //Float_t caloE = caloCluster->E()
 786         if (!fSelector->CutGeometricalAcceptance(*caloCluster)) continue;
 787         fNClusters++;
 788         //const UInt_t iPart = (UInt_t)TMath::Abs(fSelector->GetLabel(caloCluster));//->GetLabel());
 789         const UInt_t iPart = fSelector->GetLabel(caloCluster,stack);
 790         //const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
 791         TParticle *part  =  stack->Particle(iPart);
 792
 793         if (!part)
 794         {
 795             Printf("No MC particle %d", iCluster);
 796             continue;
 797         }
 798
 799         int primIdx = iPart;
 800         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
 801         {
 802             primIdx = GetPrimMother(iPart, stack);
 803             if(primIdx != stack->GetPrimary(iPart))
 804             {
 805               //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
 806               //PrintFamilyTree(iPart, stack);
 807             }
 808             //if it is from a K0S
 809             if(primIdx < 0)
 810             {
 811               //std::cout << "What!? No primary?" << std::endl;
 812               //PrintFamilyTree(iPart, stack);
 813               //continue;
 814               //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.
 815               primIdx = iPart;
 816             }
 817
 818         } // end of primary particle check
 819         //const int primCode = stack->Particle(primIdx)->GetPdgCode();
 820         TParticlePDG *pdg = part->GetPDG();
 821         if (!pdg)
 822         {
 823             Printf("ERROR: Could not get particle PDG %d", iPart);
 824             continue;
 825         }
 826
 827         //Int_t code = pdg->PdgCode();
 828 //      if(primCode == fgGammaCode)
 829 //      {
 830
 831
 832
 833         Bool_t nottrackmatched = kTRUE;//default to no track matched
 834         Float_t matchedTrackp = 0.0;
 835         Float_t matchedTrackpt = 0.0;
 836         nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
 837         //by default ALL matched tracks are accepted, whether or not the match is good.  So we check to see if the track is good.
 838         if(!nottrackmatched){//if the track is trackmatched
 839           Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
 840           if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
 841           AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
 842           matchedTrackp = track->P();
 843           matchedTrackpt = track->Pt();
 844           //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
 845           nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));//if the track is bad, this is track matched
 846           //if(!nottrackmatched) cout<<"Matched track p: "<<matchedTrackp<<" sim "<<part->P()<<endl;
 847         }
 848
 849
 850         for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
 851         {
 852           Int_t pIdx = caloCluster->GetLabelAt(i);
 853
 854           //TParticle *p = stack->Particle(pIdx);
 855
 856           if(!stack->IsPhysicalPrimary(pIdx))
 857           {
 858 //          PrintFamilyTree(pIdx, stack);
 859             pIdx = GetPrimMother(pIdx, stack);
 860           }
 861           if(fSelector->PassDistanceToBadChannelCut(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
 862           {
 863
 864 //          std::cout << "Gamma primary: " << primIdx << std::endl;
 865 //          foundGammas.push_back(primIdx);
 866             if(nottrackmatched){
 867               foundGammas.push_back(pIdx);
 868             }
 869           }
 870         }
 871         fCutFlow->Fill(cf++);
 872         if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
 873         Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster,fCentClass);
 874 //      if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
 875         if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
 876
 877
 878         fCutFlow->Fill(cf++);
 879         Float_t pos[3];
 880         //PrintFamilyTree(
 881         caloCluster->GetPosition(pos);
 882         TVector3 cp(pos);
 883
 884         TParticle *primPart = stack->Particle(primIdx);
 885         fPrimaryCode = primPart->GetPdgCode();
 886         fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
 887
 888         fPrimaryE = primPart->Energy();
 889         fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
 890         fPrimaryPx = primPart->Px();
 891         fPrimaryPy = primPart->Py();
 892         fPrimaryPz = primPart->Pz();
 893         //cout<<"I have a cluster and it's good energy "<<caloCluster->E()<<" simulated "<<fPrimaryE<<endl;
 894         fPrimaryVx = primPart->Vx();
 895         fPrimaryVy = primPart->Vy();
 896         fPrimaryVz = primPart->Vz();
 897
 898         fPrimaryAccepted = false;
 899         fPrimaryMatched = false;
 900
 901         fDepositedCode = part->GetPdgCode();
 902         fDepositedE = part->Energy();
 903         fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
 904         fDepositedCharge = (Int_t) part->GetPDG()->Charge();
 905
 906         fDepositedVx = part->Vx();
 907         fDepositedVy = part->Vy();
 908         fDepositedVz = part->Vz();
 909         fReconstructedE = caloCluster->E();
 910         fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
 911
 912         //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
 913         fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
 914         //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
 915 //      if(fSecondary)
 916 //      {
 917 //        //std::cout << "Have secondary!" << std::endl;
 918 //        //PrintFamilyTree(iPart, stack);
 919 //      }
 920
 921         pdg = primPart->GetPDG(0);
 922         //Int_t code = primPart->GetPdgCode();
 923
 924         Bool_t written = kFALSE;
 925
 926 //      Bool_t nottrackmatched = kTRUE;//default to no track matched
 927 //      nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
 928 //      //by default ALL matched tracks are accepted, whether or not the match is good.  So we check to see if the track is good.
 929 //      if(!nottrackmatched){
 930 //        Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
 931 //        if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
 932 //        AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
 933 //        //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
 934 //        nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));
 935 //      }
 936
 937         if(fSecondary){//all particles from secondary interactions
 938           written = kTRUE;
 939           if(nottrackmatched){//secondaries not removed
 940 //          Float_t vtx = TMath::Sqrt( TMath::Power(part->Vx(),2) + TMath::Power(part->Vy(),2) + TMath::Power(part->Vz(),2) );
 941
 942
 943
 944             if(!fSelector->CutGeometricalAcceptance(*part)){
 945               if(TMath::Sqrt(part->Vx() * part->Vx() + part->Vy() * part->Vy() + part->Vz()* part->Vz())>100){
 946                 //            fClusterPositionWeird->Fill(cp.Phi(), cp.PseudoRapidity());
 947                 //            //Principal arc tangent of x, in the interval [-pi/2,+pi/2] radians
 948                 //            float phiVertex = TMath::ATan(part->Vy()/part->Vx());
 949                 //            //PHOS acceptance is -0.222 pi to -0.555 pi
 950                 //            //but tan(phi) = tan(phi+pi)
 951                 //            if(part->Vx()<0) phiVertex + TMath::Pi();
 952                 //            //maximum angular distance from PHOS
 953                 //            float dphi1 = phiVertex +40.0/180.0*TMath::Pi();
 954                 //            float dphi2 = phiVertex +100.0/180.0*TMath::Pi();
 955                 //            float dphi = dphi1;
 956                 //            if(TMath::Abs(dphi1)>TMath::Abs(dphi2)) dphi= dphi2;
 957                 //            cout<<"DPhi:  "<<dphi<<" phi "<<phiVertex<<" dphi1 "<<dphi1<<" dphi2 "<<dphi2<<endl;
 958                 //            //cout<<"dphi "<<dphi1/TMath::Pi()<<"pi "<<dphi2/TMath::Pi()<<"pi phi "<<phiVertex/TMath::Pi()<<endl;
 959                 //            PrintFamilyTree(iPart, stack);
 960                 //out of acceptance clusters
 961
 962                 etSecondariesOutOfAccEffCorr += clEt;
 963               }
 964             }
 965             else{
 966               if(TMath::Sqrt(part->Vx() * part->Vx() + part->Vy() * part->Vy() + part->Vz()* part->Vz())>430){
 967                 //clusters which are in the cover of the PHOS/EMCal
 968                 etSecondariesDetectorCoverEffCorr += clEt;
 969               }
 970             }
 971
 972
 973             fSecondaryClusterEnergy->Fill(fReconstructedEt);
 974             //if(fReconstructedEt>0.3){//patch to do quick cross check THIS SHOULD NOT GET COMMITTED!!!  IF IT DID YELL AT CHRISTINE ASAP
 975 //            fHistSecondaryPositionInDetector->Fill(part->Vx(),part->Vy(),part->Vz());
 976 //            if(caloCluster->GetNLabels()>1){
 977 //              fHistSecondaryPositionInDetectorMultiple->Fill(part->Vx(),part->Vy(),part->Vz());
 978 //            }
 979               //}
 980 //          if(vtx>300){
 981 //            //cout<<"Vtx "<<vtx<<endl;
 982 //            if(fPrimaryCode==fgProtonCode ||  fPrimaryCode==fgAntiProtonCode || fPrimaryCode==fgPiPlusCode || fPrimaryCode==fgPiMinusCode || fPrimaryCode==fgKPlusCode || fPrimaryCode==fgKMinusCode){
 983 //              cout<<"I think I am really a charged hadron!"<<endl;
 984 //            }
 985 //            //PrintFamilyTree(iPart, stack);
 986 //            }
 987             fSecondaryNotRemoved++;
 988             etSecondaries += fReconstructedEt;
 989             etSecondariesEffCorr += clEt;
 990             if (fDepositedCharge != 0){//charged track not removed
 991               fChargedNotRemoved++;
 992               fEnergyChargedNotRemoved += clEt;
 993               fHistRemovedOrNot->Fill(2.0, fCentClass);
 994               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 995               fHistChargedNotRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
 996               nChargedSecondariesNotRemoved++;
 997             }
 998             else{
 999               fHistNeutralNotRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1000               nNeutralSecondariesNotRemoved++;
1001             }
1002           }
1003           else{//secondaries removed
1004             if (fDepositedCharge != 0){
1005               fChargedRemoved++;
1006               fEnergyChargedRemoved += clEt;
1007               fHistRemovedOrNot->Fill(0.0, fCentClass);
1008               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1009               fHistChargedRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1010               nChargedSecondariesRemoved++;
1011               fHistChargedTracksCut->Fill(fDepositedEt);
1012               if(fCalcTrackMatchVsMult){
1013                 fHistChargedTracksCutMult->Fill(fDepositedEt,fClusterMult);
1014                 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
1015               }
1016               fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
1017               if(fCalcTrackMatchVsMult){
1018                 if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1019                 fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1020                 fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//Fill with the efficiency corrected energy
1021               }
1022               //Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
1023               UInt_t trackindex = fSelector->GetLabel(caloCluster,stack);//(caloCluster->GetLabelsArray())->At(1);
1024               if(((UInt_t)caloCluster->GetLabel())!=trackindex){
1025               //if(fSelector->GetLabel(caloCluster,stack) !=trackindex){
1026                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1027                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1028                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1029               }
1030               else{
1031                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1032               }
1033             }
1034             else{//neutral energy removed
1035               fNeutralRemoved++;
1036               fEnergyNeutralRemoved += clEt;
1037               fHistRemovedOrNot->Fill(1.0, fCentClass);
1038               fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1039               fHistNeutralRemovedSecondaryEtVsCent->Fill(fReconstructedEt,fCentClass);
1040
1041               nNeutralSecondariesRemoved++;
1042             }
1043           }
1044         }
1045         else{//not a secondary
1046
1047           if (fDepositedCharge != 0 && fDepositedCode!=fgEMinusCode && fDepositedCode!=fgEPlusCode){//if the particle hitting the calorimeter is pi/k/p/mu
1048             written = kTRUE;
1049             fClusterMultChargedTracks++;
1050             etPiKPDeposited += clEt;
1051             if(fDepositedCode==fgProtonCode){
1052               etProtonDeposited += clEt;
1053             }
1054             if(fDepositedCode==fgAntiProtonCode){
1055               etAntiProtonDeposited += clEt;
1056             }
1057             if(nottrackmatched){//not removed but should be
1058               etPiKPDepositedNotTrackMatched += clEt;
1059               if(fDepositedCode==fgProtonCode){
1060                 etProtonDepositedNotTrackMatched += clEt;
1061               }
1062               if(fDepositedCode==fgAntiProtonCode){
1063                 etAntiProtonDepositedNotTrackMatched += clEt;
1064               }
1065               fHistHadronDepositsAll->Fill(part->Pt());
1066               fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
1067               if(fReconstructedEt>0.5) fHistHadronDepositsAllCent500MeV->Fill(part->Pt(), fCentClass);
1068               fChargedNotRemoved++;
1069               fEnergyChargedNotRemoved += clEt;
1070               fHistRemovedOrNot->Fill(2.0, fCentClass);
1071               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1072               fHistChargedTracksAccepted->Fill(fDepositedEt);
1073               if(fCalcTrackMatchVsMult){
1074                 if(matchedTrackpt<0.5){//if we could never have matched this because of its pt, how much energy did it deposit?
1075                   fHistChargedTracksAcceptedLowPtCent->Fill(fDepositedEt, fCentClass);
1076                   if(fDepositedEt>=0.5) fHistChargedTracksAcceptedLowPtCent500MeV->Fill(fDepositedEt, fCentClass);
1077                   if(pdg->PdgCode()!=fgAntiProtonCode){
1078                     fHistChargedTracksAcceptedLowPtCentNoAntiProtons->Fill(fDepositedEt, fCentClass);
1079                   }
1080                   else{
1081                     fHistChargedTracksAcceptedLowPtCentAntiProtons->Fill(fDepositedEt, fCentClass);
1082                   }
1083                 }
1084                 fHistChargedTracksAcceptedMult->Fill(fDepositedEt,fClusterMult);
1085                 if(fClusterMult<25){fHistChargedTracksAcceptedPeripheral->Fill(fDepositedEt);}
1086               }
1087             }
1088             else{//removed and should have been
1089               Int_t trackindex =  fSelector->GetLabel(caloCluster,stack);// (caloCluster->GetLabelsArray())->At(0);
1090               fHistHadronDepositsReco->Fill(part->Pt());
1091               fHistHadronDepositsRecoCent->Fill(part->Pt(), fCentClass);
1092               fHistHadronDepositsAll->Fill(part->Pt());
1093               fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
1094               if(fReconstructedEt>0.5) fHistHadronDepositsAllCent500MeV->Fill(part->Pt(), fCentClass);
1095               if(caloCluster->GetLabel()!=trackindex){
1096                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1097                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1098                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1099               }
1100               else{
1101                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1102               }
1103               fChargedRemoved++;
1104               fEnergyChargedRemoved += clEt;
1105               fHistRemovedOrNot->Fill(0.0, fCentClass);
1106               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1107               fHistChargedTracksCut->Fill(fDepositedEt);
1108               if(fCalcTrackMatchVsMult){
1109                 fHistChargedTracksCutMult->Fill(fDepositedEt,fClusterMult);
1110                 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
1111               }
1112               fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
1113               if(fCalcTrackMatchVsMult){
1114                 if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1115                 fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1116                 fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//fill with the efficiency corrected energy
1117               }
1118             }
1119           }
1120           //K0L and any neutral particles from the decay of K+/- or K0S
1121           if(!written && (fPrimaryCode==fgKPlusCode || fPrimaryCode==fgKMinusCode || fPrimaryCode==fgK0SCode ||fPrimaryCode==fgK0LCode)){
1122             written = kTRUE;//At this point we are not tracking them but we don't count them as neutrals accidentally removed
1123           }
1124
1125           if(!written && (fDepositedCode==fgGammaCode || fDepositedCode==fgEMinusCode || fDepositedCode ==fgEPlusCode)){//if the particle hitting the calorimeter is gamma, electron and not from a kaon
1126             fClusterMultGammas++;
1127             written = kTRUE;
1128             if(nottrackmatched){//Not removed and not supposed to be removed - signal
1129               fEtNonRemovedGammas += clEt;
1130               fMultNonRemovedGammas++;
1131               fNeutralNotRemoved--;
1132               fEnergyNeutralNotRemoved -= clEt;
1133               fHistGammasAccepted->Fill(fDepositedEt);
1134               if(fCalcTrackMatchVsMult){
1135                 fHistGammasAcceptedMult->Fill(fDepositedEt,fClusterMult);
1136                 if(fClusterMult<25){fHistGammasAcceptedPeripheral->Fill(fDepositedEt);}
1137               }
1138             }
1139             else{//removed but shouldn't have been
1140               Int_t trackindex = fSelector->GetLabel(caloCluster,stack);// (caloCluster->GetLabelsArray())->At(1);
1141               if(caloCluster->GetLabel()!=trackindex){
1142                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
1143                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
1144 //              cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
1145 //              PrintFamilyTree(trackindex, stack);
1146 //              cout<<"Cluster"<<endl;
1147               }
1148               fGammaRemoved++;
1149               fGammaRemovedEt+=clEt;
1150               fHistGammasCut->Fill(fDepositedEt);
1151               if(fCalcTrackMatchVsMult){
1152                 fHistGammasCutMult->Fill(fDepositedEt,fClusterMult);
1153                 if(fClusterMult<25){fHistGammasCutPeripheral->Fill(fDepositedEt);}
1154               }
1155               fHistMatchedTracksEvspTSignal->Fill(matchedTrackp,fReconstructedE);
1156               if(fCalcTrackMatchVsMult){
1157                 if(fClusterMult<25){fHistMatchedTracksEvspTSignalPeripheral->Fill(matchedTrackp,fReconstructedEt);}
1158                 fHistMatchedTracksEvspTSignalvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
1159                 fHistMatchedTracksEvspTSignalvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);
1160               }
1161             }
1162           }
1163           //all other cases - neutron, anti-neutron, not aware of other cases
1164           if(!written){
1165             fNeutralNotRemoved++;
1166             fEnergyNeutralNotRemoved += clEt;//this is the efficiency corrected energy
1167             fHistRemovedOrNot->Fill(3.0, fCentClass);
1168             if (fDepositedCode == fgNeutronCode || fDepositedCode == fgAntiNeutronCode){
1169               fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
1170               fHistNeutronsEtVsCent->Fill(clEt,fCentClass);
1171               nNeutrons++;
1172               if(fDepositedCode == fgNeutronCode){
1173                 etNeutronDeposited += clEt;
1174               }
1175               if(fDepositedCode == fgAntiNeutronCode){
1176                 etAntiNeutronDeposited += clEt;
1177               }
1178               //cout<<"I am a";
1179               //if(fDepositedCode == fgAntiNeutronCode) cout<<"n anti-";
1180               //else{cout<<" ";}
1181               //cout<<"neutron!! pt "<<part->Pt()<<" eta "<<part->Eta()<<" phi "<<part->Phi()<<" Et "<<clEt<<endl;
1182               //PrintFamilyTree(iPart, stack);
1183             }
1184             else{
1185               nNotNeutrons++;
1186 //            cout<<"I am a";
1187 //            if(fDepositedCode == fgAntiNeutronCode) cout<<"n anti-";
1188 //            else{cout<<" ";}
1189 //            cout<<"neutron!! pt "<<part->Pt()<<" eta "<<part->Eta()<<" phi "<<part->Phi()<<" Et "<<clEt<<endl;
1190 //            PrintFamilyTree(iPart, stack);
1191 //            cout<<"I am not a neutron!!"<<endl;
1192 //            PrintFamilyTree(iPart, stack);
1193             }
1194           }
1195         }
1196         fPrimaryTree->Fill();
1197     } // end of loop over clusters
1198     fHistPiKPNotTrackMatchedDepositedVsNch->Fill(etPiKPDepositedNotTrackMatched,multiplicity);
1199     fHistPiKPDepositedVsNch->Fill(etPiKPDeposited,multiplicity);
1200     fHistNeutronsDepositedVsNch->Fill(etNeutronDeposited,multiplicity);
1201     fHistAntiNeutronsDepositedVsNch->Fill(etAntiNeutronDeposited,multiplicity);
1202     fHistProtonsDepositedVsNch->Fill(etProtonDeposited,multiplicity);
1203     fHistAntiProtonsDepositedVsNch->Fill(etAntiProtonDeposited,multiplicity);
1204     fHistProtonsNotTrackMatchedDepositedVsNch->Fill(etProtonDepositedNotTrackMatched,multiplicity);
1205     fHistAntiProtonsNotTrackMatchedDepositedVsNch->Fill(etAntiProtonDepositedNotTrackMatched,multiplicity);
1206     fHistSecondariesVsNch->Fill(etSecondaries,multiplicity);
1207     fHistSecondariesVsNcl->Fill(etSecondaries,fNClusters);
1208     fHistSecondariesEffCorrVsNch->Fill(etSecondariesEffCorr,multiplicity);
1209     fHistSecondariesEffCorrVsNcl->Fill(etSecondariesEffCorr,fNClusters);
1210     fHistSecondariesOutOfAccEffCorrVsNch->Fill(etSecondariesOutOfAccEffCorr,multiplicity);
1211     fHistSecondariesDetectorCoverEffCorrVsNch->Fill(etSecondariesDetectorCoverEffCorr,multiplicity);
1212     fHistCentVsNchVsNcl->Fill(fCentClass,multiplicity, fNClusters);
1213
1214     fHistNeutronsNumVsCent->Fill(nNeutrons,fCentClass);
1215     fHistNotNeutronsNumVsCent->Fill(nNotNeutrons,fCentClass);
1216
1217     std::sort(foundGammas.begin(), foundGammas.end());
1218     for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
1219     {
1220
1221         if(!stack->IsPhysicalPrimary(iPart)) continue;
1222
1223         TParticle *part = stack->Particle(iPart);
1224
1225         if (!part)
1226         {
1227             Printf("ERROR: Could not get particle %d", iPart);
1228             continue;
1229         }
1230         TParticlePDG *pdg = part->GetPDG(0);
1231
1232         if (!pdg)
1233         {
1234             Printf("ERROR: Could not get particle PDG %d", iPart);
1235             continue;
1236         }
1237
1238         if(pdg->PdgCode()==fgGammaCode)// TMath::Abs(part->Eta()) < 0.12)
1239         {
1240           if(fSelector->CutGeometricalAcceptance(*part)){
1241             fHistGammasGenerated->Fill(part->Energy());
1242             fHistGammasGeneratedCent->Fill(part->Energy(),fCentClass);
1243           }
1244           if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
1245           {
1246             if(!fSelector->CutGeometricalAcceptance(*part)){
1247               //cout<<"Gamma NOT in acceptance"<<endl;
1248               fHistGammasFoundOutOfAccCent->Fill(part->Energy(),fCentClass);
1249             }
1250             else{
1251               fHistGammasFound->Fill(part->Energy());
1252               fHistGammasFoundCent->Fill(part->Energy(),fCentClass);
1253               //cout<<"Gamma IN acceptance"<<endl;
1254             }
1255           }
1256         }
1257         if(pdg->PdgCode()==fgPiPlusCode || pdg->PdgCode()==fgPiMinusCode || pdg->PdgCode()==fgProtonCode || pdg->PdgCode()==fgAntiProtonCode){//section here for all hadrons generated
1258           fHistHadronsAllCent->Fill(part->Pt(), fCentClass);
1259         }
1260
1261
1262     }
1263     if(fCalcForKaonCorrection){
1264       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};
1265       Int_t nEtCuts = 11;
1266       //loop over simulated particles in order to find K0S
1267       for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++){
1268         TParticle *part = stack->Particle(iPart);
1269         if (!part){
1270           //Printf("ERROR: Could not get particle %d", iPart);
1271           continue;
1272         }
1273         TParticlePDG *pdg = part->GetPDG(0);
1274         if (!pdg){
1275           //Printf("ERROR: Could not get particle PDG %d", iPart);
1276           continue;
1277         }
1278         //if(stack->IsPhysicalPrimary(iPart)){//if it is a K0 it might have decayed into four pions
1279         //fgK0SCode, fgGammaCode, fgPi0Code
1280         Int_t code = pdg->PdgCode();
1281         if(code == fgK0SCode || code==fgKPlusCode || code==fgKMinusCode ||code==fgK0LCode || code==fgK0Code){//this is a kaon
1282           //cout<<"I am a kaon too! "<<stack->Particle(iPart)->GetName()<<" "<<code<<endl;
1283           Float_t pTk = stack->Particle(iPart)->Pt();
1284           if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//these are particles which would be included in our spectra measurements
1285             fHistSimKaonsInAcceptance->Fill(pTk);
1286             if(code == fgK0SCode){fHistSimK0SInAcceptance->Fill(pTk);}
1287             if(code == fgK0LCode){fHistSimK0LInAcceptance->Fill(pTk);}
1288             if(code == fgKPlusCode){fHistSimKPlusInAcceptance->Fill(pTk);}
1289             if(code == fgKMinusCode){fHistSimKMinusInAcceptance->Fill(pTk);}
1290             if(code == fgK0Code){//Split K0's between the two
1291               fHistSimK0SInAcceptance->Fill(pTk,0.5);
1292               fHistSimK0LInAcceptance->Fill(pTk,0.5);
1293             }
1294           }
1295           else{
1296             fHistSimKaonsOutOfAcceptance->Fill(pTk);
1297             //      if(!stack->IsPhysicalPrimary(iPart)){
1298             //        PrintFamilyTree(iPart, stack);
1299             //      }
1300           }
1301           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};
1302           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};
1303           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...
1304             AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
1305             if (!fSelector->CutGeometricalAcceptance(*caloCluster)) continue;
1306             const Int_t myPart = fSelector->GetLabel(caloCluster,stack);
1307             //const Int_t myPart = TMath::Abs(caloCluster->GetLabel());
1308             //identify the primary particle which created this cluster
1309             int primIdx = myPart;
1310             if (!stack->IsPhysicalPrimary(myPart)){
1311               primIdx = GetPrimMother(iPart, stack);
1312             } // end of primary particle check
1313             TParticle *hitPart = stack->Particle(myPart);
1314             Bool_t hitsAsChargedKaon = kFALSE;
1315             if(hitPart->GetPdgCode()== fgKPlusCode || hitPart->GetPdgCode()== fgKPlusCode){
1316               if(myPart==primIdx){
1317                 //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!
1318                 hitsAsChargedKaon = kTRUE;
1319                 //cout<<"Found primary charged kaon cluster!"<<endl;
1320               }
1321             }
1322             if(primIdx==iPart && primIdx>0 && !hitsAsChargedKaon){//This cluster is from our primary particle and our primary particle is a kaon
1323               //cout<<"I have a particle match! prim code"<<code<<" id "<<primIdx <<endl;
1324               Float_t pos[3];
1325               caloCluster->GetPosition(pos);
1326               TVector3 cp(pos);
1327               Double_t clEt = caloCluster->E()*TMath::Sin(cp.Theta());
1328               Double_t clEtCorr = CorrectForReconstructionEfficiency(*caloCluster,fCentClass);
1329               for(int l=0;l<nEtCuts;l++){//loop over cut values
1330                 if(clEt>=etCuts[l]){
1331                   //cout<<", "<<clEt<<">="<<etCuts[l];
1332                   totalClusterEts[l] += clEtCorr;//if cluster et is above the cut off energy add it
1333                   totalGammaEts[l] += clEt;//if cluster et is above the cut off energy add it
1334                 }
1335               }
1336             }
1337           }
1338           //      cout<<"Deposits:  pT: "<<pTk;
1339           //      for(int l=0;l<nEtCuts;l++){//loop over cut values
1340           //        cout<<" "<<totalClusterEts[l];
1341           //      }
1342           //      cout<<endl;
1343           if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//within the acceptance of our spectra and is a primary particle
1344             if(totalClusterEts[0]>0.0){fHistSimKaonsInAcceptanceWithDepositsPrimaries->Fill(pTk);}
1345             //cout<<"I have a particle match! prim code"<<code<<" id "<<iPart <<endl;
1346             for(int l=0;l<nEtCuts;l++){
1347               fHistK0EDepositsVsPtInAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1348               fHistK0EGammaVsPtInAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
1349             }
1350           }
1351           else{//outside the acceptance of our spectra
1352             if(totalClusterEts[0]>0.0){
1353               if(stack->IsPhysicalPrimary(iPart)){fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries->Fill(pTk);}
1354               else{fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries->Fill(pTk);}
1355             }
1356             for(int l=0;l<nEtCuts;l++){
1357               fHistK0EDepositsVsPtOutOfAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1358               fHistK0EGammaVsPtOutOfAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
1359             }
1360           }
1361
1362         }
1363       }
1364     }
1365     fHistMultChVsSignalVsMult->Fill(fClusterMultChargedTracks,fClusterMultGammas,fNClusters);
1366     fHistNeutralRemovedSecondaryNumVsNCluster->Fill(nNeutralSecondariesRemoved,fNClusters);
1367     fHistChargedRemovedSecondaryNumVsNCluster->Fill(nChargedSecondariesRemoved,fNClusters);
1368     fHistNeutralNotRemovedSecondaryNumVsNCluster->Fill(nNeutralSecondariesNotRemoved,fNClusters);
1369     fHistChargedNotRemovedSecondaryNumVsNCluster->Fill(nChargedSecondariesNotRemoved,fNClusters);
1370     fHistNeutralRemovedSecondaryNumVsCent->Fill(nNeutralSecondariesRemoved,fCentClass);
1371     fHistChargedRemovedSecondaryNumVsCent->Fill(nChargedSecondariesRemoved,fCentClass);
1372     fHistNeutralNotRemovedSecondaryNumVsCent->Fill(nNeutralSecondariesNotRemoved,fCentClass);
1373     fHistChargedNotRemovedSecondaryNumVsCent->Fill(nChargedSecondariesNotRemoved,fCentClass);
1374     //cout<<"Secondaries not removed: "<<fSecondaryNotRemoved<<" neutral secondaries not removed "<<nNeutralSecondariesNotRemoved<<" cluster mult "<<fClusterMult<<endl;
1375     FillHistograms();
1376     return 0;
1377 }
1378
1379 void AliAnalysisEtMonteCarlo::Init()
1380 {   // init
1381     AliAnalysisEt::Init();
1382 }
1383
1384 void AliAnalysisEtMonteCarlo::ResetEventValues()
1385 {   // reset event values
1386     AliAnalysisEt::ResetEventValues();
1387   if(!fIsMC) return;
1388
1389     fTotEtSecondary = 0;
1390     fTotEtSecondaryFromEmEtPrimary = 0;
1391     fTotEtWithSecondaryRemoved = 0;
1392
1393     // collision geometry defaults for p+p:
1394     fImpactParameter = 0;
1395     fNcoll = 1;
1396     fNpart = 2;
1397
1398     fEtNonRemovedProtons = 0;
1399     fEtNonRemovedAntiProtons = 0;
1400     fEtNonRemovedPiPlus = 0;
1401     fEtNonRemovedPiMinus = 0;
1402     fEtNonRemovedKaonPlus = 0;
1403     fEtNonRemovedKaonMinus = 0;
1404     fEtNonRemovedK0S = 0;
1405     fEtNonRemovedK0L = 0;
1406     fEtNonRemovedLambdas = 0;
1407     fEtNonRemovedElectrons = 0;
1408     fEtNonRemovedPositrons = 0;
1409     fEtNonRemovedMuPlus = 0;
1410     fEtNonRemovedMuMinus = 0;
1411     fEtNonRemovedNeutrons = 0;
1412     fEtNonRemovedAntiNeutrons = 0;
1413     fEtNonRemovedGammas = 0;
1414     fEtNonRemovedGammasFromPi0 = 0;
1415
1416     fEtRemovedProtons = 0;
1417     fEtRemovedAntiProtons = 0;
1418     fEtRemovedPiPlus = 0;
1419     fEtRemovedPiMinus = 0;
1420     fEtRemovedKaonPlus = 0;
1421     fEtRemovedKaonMinus = 0;
1422     fEtRemovedK0s = 0;
1423     fEtRemovedK0L = 0;
1424     fEtRemovedLambdas = 0;
1425     fEtRemovedElectrons = 0;
1426     fEtRemovedPositrons = 0;
1427     fEtRemovedMuPlus = 0;
1428     fEtRemovedMuMinus = 0;
1429     fEtRemovedNeutrons = 0;
1430
1431     fEtRemovedGammasFromPi0 = 0;
1432     fEtRemovedGammas = 0;
1433     fEtRemovedNeutrons = 0;
1434     fEtRemovedAntiNeutrons = 0;
1435
1436     fMultNonRemovedProtons = 0;
1437     fMultNonRemovedAntiProtons = 0;
1438     fMultNonRemovedPiPlus = 0;
1439     fMultNonRemovedPiMinus = 0;
1440     fMultNonRemovedKaonPlus = 0;
1441     fMultNonRemovedKaonMinus = 0;
1442     fMultNonRemovedK0s = 0;
1443     fMultNonRemovedK0L = 0;
1444     fMultNonRemovedLambdas = 0;
1445     fMultNonRemovedElectrons = 0;
1446     fMultNonRemovedPositrons = 0;
1447     fMultNonRemovedMuPlus = 0;
1448     fMultNonRemovedMuMinus = 0;
1449     fMultNonRemovedNeutrons = 0;
1450     fMultNonRemovedAntiNeutrons = 0;
1451     fMultNonRemovedGammas = 0;
1452
1453     fMultRemovedProtons = 0;
1454     fMultRemovedAntiProtons = 0;
1455     fMultRemovedPiPlus = 0;
1456     fMultRemovedPiMinus = 0;
1457     fMultRemovedKaonPlus = 0;
1458     fMultRemovedKaonMinus = 0;
1459     fMultRemovedK0s = 0;
1460     fMultRemovedK0L = 0;
1461     fMultRemovedLambdas = 0;
1462     fMultRemovedElectrons = 0;
1463     fMultRemovedPositrons = 0;
1464     fMultRemovedMuPlus = 0;
1465     fMultRemovedMuMinus = 0;
1466
1467     fMultRemovedGammas = 0;
1468     fMultRemovedNeutrons = 0;
1469     fMultRemovedAntiNeutrons = 0;
1470
1471     fEnergyChargedNotRemoved = 0;
1472     fEnergyChargedRemoved = 0;
1473     fEnergyNeutralNotRemoved = 0;
1474     fEnergyNeutralRemoved = 0;
1475
1476     fChargedNotRemoved = 0;
1477     fChargedRemoved = 0;
1478     fNeutralNotRemoved = 0;
1479     fNeutralRemoved = 0;
1480     fGammaRemoved = 0;
1481     fSecondaryNotRemoved = 0;
1482
1483     fTrackMultInAcc = 0;
1484
1485     fTotNeutralEtAfterMinEnergyCut = 0;
1486
1487     fSecondaryNotRemoved = 0;
1488
1489     fTotPx = 0;
1490     fTotPy = 0;
1491     fTotPz = 0;
1492
1493
1494 }
1495
1496 void AliAnalysisEtMonteCarlo::CreateHistograms()
1497 {   // histogram related additions
1498     AliAnalysisEt::CreateHistograms();
1499     if(!fIsMC) return;
1500     if (fEventSummaryTree) {
1501         fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1502         fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1503         fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1504         fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1505         fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1506         fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1507         fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
1508         fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1509         fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1510         fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1511         fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1512         fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1513         fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1514         fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1515         fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1516         fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
1517 //      fEventSummaryTree->Branch("f
1518     }
1519
1520     //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1521     //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1522     //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1523     //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1524
1525     fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1526
1527     fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1528     fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1529     fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1530     fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1531     fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1532     fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1533     fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
1534     fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
1535     fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1536     fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1537     fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1538     fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1539     fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1540     fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1541     fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1542
1543     fHistEtNonRemovedGammas = new  TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1544     fHistEtNonRemovedGammasFromPi0 = new  TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1545
1546     fHistEtRemovedGammas = new  TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1547     fHistEtRemovedNeutrons = new  TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1548     fHistEtRemovedAntiNeutrons = new  TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1549
1550     fHistEtRemovedCharged = new  TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1551     fHistEtRemovedNeutrals = new  TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1552
1553     fHistEtNonRemovedCharged = new  TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1554     fHistEtNonRemovedNeutrals = new  TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1555
1556     fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1557     fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1558     fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1559     fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1560     fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1561     fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1562     fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
1563     fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
1564     fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1565     fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1566     fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1567     fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1568     fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1569     fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1570     fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1571
1572     fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1573
1574     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1575     fHistMultRemovedNeutrons = new  TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1576     fHistMultRemovedAntiNeutrons = new  TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1577     /*
1578       fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1579       fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1580
1581       fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1582       fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
1583
1584
1585     fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1586     fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1587
1588     fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1589     fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1590
1591     fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1592     fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1593     fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1594
1595     fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1596     fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1597     fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1598
1599     fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1600     fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1601     fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1602     fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1603
1604     if(fCalcForKaonCorrection){
1605       Int_t nEtCut = 11;
1606       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};
1607       fHistK0EDepositsVsPtInAcceptance = new TH3F("fHistK0EDepositsVsPtInAcceptance","Kaon deposits with corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1608       fHistK0EGammaVsPtInAcceptance = new TH3F("fHistK0EGammaVsPtInAcceptance","Kaon deposits without corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1609       fHistK0EDepositsVsPtOutOfAcceptance = new TH3F("fHistK0EDepositsVsPtOutOfAcceptance","Kaon deposits with corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1610       fHistK0EGammaVsPtOutOfAcceptance = new TH3F("fHistK0EGammaVsPtOutOfAcceptance","Kaon deposits without corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1611       fHistK0EDepositsVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1612       fHistK0EGammaVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1613       fHistK0EDepositsVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1614       fHistK0EGammaVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1615       fHistK0EDepositsVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1616       fHistK0EGammaVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1617       fHistK0EDepositsVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1618       fHistK0EGammaVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1619       fHistK0EDepositsVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1620       fHistK0EGammaVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1621       fHistK0EDepositsVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1622       fHistK0EGammaVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1623
1624       fHistSimKaonsInAcceptance = new TH1F("fHistSimKaonsInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1625       fHistSimK0SInAcceptance = new TH1F("fHistSimK0SInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1626       fHistSimK0LInAcceptance = new TH1F("fHistSimK0LInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1627       fHistSimKPlusInAcceptance = new TH1F("fHistSimKPlusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1628       fHistSimKMinusInAcceptance = new TH1F("fHistSimKMinusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1629       fHistSimKaonsOutOfAcceptance = new TH1F("fHistSimKaonsOutOfAcceptance","Kaons with y>0.5",fgNumOfPtBins,fgPtAxis);
1630       fHistSimKaonsInAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsInAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1631       fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries","Secondary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1632       fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1633     }
1634
1635     fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1636     fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1637     fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1638
1639     fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1640     fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1641     fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1642
1643     if(fCuts->GetHistMakeTree())
1644     {
1645         TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1646         fPrimaryTree = new TTree(treename, treename);
1647
1648         fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1649         fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1650         fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1651
1652         fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1653         fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1654
1655         fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1656         fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1657
1658         fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1659         fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1660         fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1661
1662         fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1663         fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1664         fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1665
1666         fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
1667         fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
1668
1669
1670         fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1671         fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
1672         fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
1673         fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1674
1675         fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1676         fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1677         fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1678
1679         fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
1680
1681
1682         fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
1683         fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
1684
1685         fPrimaryTree->Branch("fClusterMult", &fClusterMult,  "fClusterMult/I");
1686
1687
1688     }
1689
1690     fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
1691     fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
1692     fHistGammasFoundOutOfAccCent = new TH2F("fHistGammasFoundOutOfAccCent", "fHistGammasFoundOutOfAccCent",200, 0, 10,20,-0.5,19.5);
1693     fHistGammasFoundCent = new TH2F("fHistGammasFoundCent", "fHistGammasFoundCent",200, 0, 10,20,-0.5,19.5);
1694     fHistGammasGeneratedCent = new TH2F("fHistGammasGeneratedCent", "fHistGammasGeneratedCent",200, 0, 10,20,-0.5,19.5);
1695     fHistChargedTracksCut = new TH1F("fHistChargedTracksCut", "fHistChargedTracksCut",100, 0, 5);
1696     fHistChargedTracksAccepted = new TH1F("fHistChargedTracksAccepted", "fHistChargedTracksAccepted",100, 0, 5);
1697     fHistGammasCut = new TH1F("fHistGammasTracksCut", "fHistGammasTracksCut",100, 0, 5);
1698     fHistGammasAccepted = new TH1F("fHistGammasTracksAccepted", "fHistGammasTracksAccepted",100, 0, 5);
1699
1700     if(fCalcTrackMatchVsMult){
1701       fHistChargedTracksCutMult = new TH2F("fHistChargedTracksCutMult", "fHistChargedTracksCutMult",100, 0, 5,10,0,100);
1702       fHistChargedTracksAcceptedMult = new TH2F("fHistChargedTracksAcceptedMult", "fHistChargedTracksAcceptedMult",100, 0, 5,10,0,100);
1703       fHistChargedTracksAcceptedLowPtCent = new TH2F("fHistChargedTracksAcceptedLowPtCent", "fHistChargedTracksAcceptedLowPtCent",100, 0, 5,20,-0.5,19.5);
1704       fHistChargedTracksAcceptedLowPtCent500MeV = new TH2F("fHistChargedTracksAcceptedLowPtCent500MeV", "fHistChargedTracksAcceptedLowPtCent500MeV",100, 0, 5,20,-0.5,19.5);
1705       fHistChargedTracksAcceptedLowPtCentNoAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentNoAntiProtons", "fHistChargedTracksAcceptedLowPtCentNoAntiProtons",100, 0, 5,20,-0.5,19.5);
1706       fHistChargedTracksAcceptedLowPtCentAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentAntiProtons", "fHistChargedTracksAcceptedLowPtCentAntiProtons",100, 0, 5,20,-0.5,19.5);
1707       fHistGammasCutMult = new TH2F("fHistGammasTracksCutMult", "fHistGammasTracksCutMult",100, 0, 5,10,0,100);
1708       fHistGammasAcceptedMult = new TH2F("fHistGammasTracksAcceptedMult", "fHistGammasTracksAcceptedMult",100, 0, 5,10,0,100);
1709     }
1710
1711     fHistBadTrackMatches = new TH1F("fHistBadTrackMatches", "fHistBadTrackMatches",100, 0, 5);
1712     fHistMatchedTracksEvspTBkgd = new TH2F("fHistMatchedTracksEvspTBkgd", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1713     fHistMatchedTracksEvspTSignal = new TH2F("fHistMatchedTracksEvspTSignal", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
1714     if(fCalcTrackMatchVsMult){
1715       fHistMatchedTracksEvspTBkgdPeripheral = new TH2F("fHistMatchedTracksEvspTBkgdPeripheral", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1716       fHistMatchedTracksEvspTSignalPeripheral = new TH2F("fHistMatchedTracksEvspTSignalPeripheral", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
1717
1718       fHistMatchedTracksEvspTBkgdvsCent = new TH3F("fHistMatchedTracksEvspTBkgdvsCent", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
1719       fHistMatchedTracksEvspTSignalvsCent = new TH3F("fHistMatchedTracksEvspTSignalvsCent", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
1720       fHistMatchedTracksEvspTBkgdvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTBkgdvsCentEffCorr", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
1721       fHistMatchedTracksEvspTSignalvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTSignalvsCentEffCorr", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
1722
1723
1724       fHistChargedTracksCutPeripheral = new TH1F("fHistChargedTracksCutPeripheral", "fHistChargedTracksCut",100, 0, 5);
1725       fHistChargedTracksAcceptedPeripheral = new TH1F("fHistChargedTracksAcceptedPeripheral", "fHistChargedTracksAccepted",100, 0, 5);
1726       fHistGammasCutPeripheral = new TH1F("fHistGammasTracksCutPeripheral", "fHistGammasTracksCut",100, 0, 5);
1727       fHistGammasAcceptedPeripheral = new TH1F("fHistGammasTracksAcceptedPeripheral", "fHistGammasTracksAccepted",100, 0, 5);
1728     }
1729     fHistBadTrackMatchesdPhidEta = new TH2F("fHistBadTrackMatchesdPhidEta", "fHistBadTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
1730     fHistGoodTrackMatchesdPhidEta = new TH2F("fHistGoodTrackMatchesdPhidEta", "fHistGoodTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
1731
1732     fHistHadronDepositsAll = new TH1F("fHistHadronDepositsAll","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
1733     fHistHadronDepositsReco = new TH1F("fHistHadronDepositsReco","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
1734     //,10,0,100
1735       Int_t nMult = 20;
1736       Float_t nMultCuts[21] = { 0, 5,10,15,20, 25,30,35,40,45,
1737                                50,55,60,65,70, 75,80,85,90,95,
1738                                 100};
1739       Int_t nCent = 20;
1740       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};
1741
1742       fHistHadronDepositsAllCent = new TH2F("fHistHadronDepositsAllCent","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
1743       fHistHadronDepositsAllCent500MeV = new TH2F("fHistHadronDepositsAllCent500MeV","All Hadrons which deposited energy in calorimeter pT>500MeV",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
1744       fHistHadronDepositsRecoCent = new TH2F("fHistHadronDepositsRecoCent","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
1745
1746       fHistHadronsAllCent = new TH2F("fHistHadronsAllCent","All Hadrons vs cluster mult",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
1747
1748       fHistMultChVsSignalVsMult = new TH3F("fHistMultChVsSignalVsMult","Charged particle Multiplicity vs Signal particle multiplicity vs Cluster Mult",nMult,nMultCuts,nMult,nMultCuts,nMult,nMultCuts);
1749       fHistNeutralRemovedSecondaryEtVsCent = new TH2F("fHistNeutralRemovedSecondaryEtVsCent","Neutral Removed Secondaries E_{T} vs centrality",100,0.0,2.0,20,-0.5,19.5);
1750       fHistChargedRemovedSecondaryEtVsCent = new TH2F("fHistChargedRemovedSecondaryEtVsCent","Charged Removed Secondaries E_{T} vs centrality",100,0.0,2.0,20,-0.5,19.5);
1751       fHistNeutralNotRemovedSecondaryEtVsCent = new TH2F("fHistNeutralNotRemovedSecondaryEtVsCent","Neutral NotRemoved Secondaries E_{T} vs centrality",100,0.0,2.0,20,-0.5,19.5);
1752       fHistChargedNotRemovedSecondaryEtVsCent = new TH2F("fHistChargedNotRemovedSecondaryEtVsCent","Charged NotRemoved Secondaries E_{T} vs centrality",100,0.0,2.0,20,-0.5,19.5);
1753       fHistNeutralRemovedSecondaryNumVsNCluster = new TH2F("fHistNeutralRemovedSecondaryNumVsNCluster","Neutral Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,250,0,250);
1754       fHistChargedRemovedSecondaryNumVsNCluster = new TH2F("fHistChargedRemovedSecondaryNumVsNCluster","Charged Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,250,0,250);
1755       fHistNeutralNotRemovedSecondaryNumVsNCluster = new TH2F("fHistNeutralNotRemovedSecondaryNumVsNCluster","Neutral NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,250,0,250);
1756       fHistChargedNotRemovedSecondaryNumVsNCluster = new TH2F("fHistChargedNotRemovedSecondaryNumVsNCluster","Charged NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,250,0,250);
1757       fHistNeutralRemovedSecondaryNumVsCent = new TH2F("fHistNeutralRemovedSecondaryNumVsCent","Neutral Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
1758       fHistChargedRemovedSecondaryNumVsCent = new TH2F("fHistChargedRemovedSecondaryNumVsCent","Charged Removed Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
1759       fHistNeutralNotRemovedSecondaryNumVsCent = new TH2F("fHistNeutralNotRemovedSecondaryNumVsCent","Neutral NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
1760       fHistChargedNotRemovedSecondaryNumVsCent = new TH2F("fHistChargedNotRemovedSecondaryNumVsCent","Charged NotRemoved Secondaries Number vs N_{cluster}",20,-0.5,19.5,20,-0.5,19.5);
1761       fHistNeutronsEtVsCent = new TH2F("fHistNeutronsEtVsCent","Neutrons and anti-neutrons - deposited ET vs Centrality bin",100,0,4.0,20,-0.5,19.5);
1762       fHistNeutronsNumVsCent = new TH2F("fHistNeutronsNumVsCent","Neutrons and anti-neutrons - number vs Centrality bin",20,-0.5,19.5,20,-0.5,19.5);
1763       fHistNotNeutronsNumVsCent = new TH2F("fHistNotNeutronsNumVsCent","Neutral particles not otherwise classified - number vs Centrality bin",20,-0.5,19.5,20,-0.5,19.5);
1764       Int_t nbinsEt = 125;
1765       Float_t maxEtRange = 125;
1766       Float_t maxEtRangeShort = 25;
1767       Float_t minEtRange = 0;
1768       Int_t nbinsMult = 100;
1769       Float_t maxMult = 3000;
1770       Float_t minMult = 0;
1771       Int_t nbinsCl = 175;
1772       Float_t maxCl = 350;
1773       Float_t minCl = 0;
1774       fHistPiKPDepositedVsNch = new TH2F("fHistPiKPDepositedVsNch","#pi,K,p E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
1775       fHistPiKPNotTrackMatchedDepositedVsNch = new TH2F("fHistPiKPNotTrackMatchedDepositedVsNch","#pi,K,p E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
1776     fHistNeutronsDepositedVsNch = new TH2F("fHistNeutronsDepositedVsNch","n deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1777     fHistAntiNeutronsDepositedVsNch = new TH2F("fHistAntiNeutronsDepositedVsNch","#bar{n} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1778     fHistProtonsDepositedVsNch = new TH2F("fHistProtonsDepositedVsNch","p deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1779     fHistAntiProtonsDepositedVsNch = new TH2F("fHistAntiProtonsDepositedVsNch","#bar{p} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1780     fHistProtonsNotTrackMatchedDepositedVsNch = new TH2F("fHistProtonsNotTrackMatchedDepositedVsNch","p not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1781     fHistAntiProtonsNotTrackMatchedDepositedVsNch = new TH2F("fHistAntiProtonsNotTrackMatchedDepositedVsNch","#bar{p} not track matched deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1782     fHistSecondariesVsNch = new TH2F("fHistSecondariesVsNch","secondaries deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1783     fHistSecondariesVsNcl = new TH2F("fHistSecondariesVsNcl","secondaries deposited in calorimeter vs number of clusters",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
1784     fHistSecondariesEffCorrVsNch = new TH2F("fHistSecondariesEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRangeShort,nbinsMult,minMult,maxMult);
1785     fHistSecondariesEffCorrVsNcl = new TH2F("fHistSecondariesEffCorrVsNcl","efficiency corrected secondaries deposited in calorimeter vs number of clusters",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
1786     fHistSecondariesOutOfAccEffCorrVsNch = new TH2F("fHistSecondariesOutOfAccEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs number of clusters for secondary particles out of detector acceptance",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
1787     fHistSecondariesDetectorCoverEffCorrVsNch = new TH2F("fHistSecondariesDetectorCoverEffCorrVsNch","efficiency corrected secondaries deposited in calorimeter vs number of clusters for secondaries from the detector cover",nbinsEt,minEtRange,maxEtRangeShort,nbinsCl,minCl,maxCl);
1788     fHistCentVsNchVsNcl = new TH3F("fHistCentVsNchVsNcl","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
1789     //float maxpos = 500;
1790     // int nbinspos = 200;
1791 //     fHistSecondaryPositionInDetector = new TH3F("fHistSecondaryPositionInDetector","Position of secondaries",nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos);
1792 //     fHistSecondaryPositionInDetector->GetXaxis()->SetTitle("X");
1793 //     fHistSecondaryPositionInDetector->GetYaxis()->SetTitle("Y");
1794 //     fHistSecondaryPositionInDetector->GetZaxis()->SetTitle("Z");
1795 //     fHistSecondaryPositionInDetectorMultiple = new TH3F("fHistSecondaryPositionInDetectorMultiple","Position of secondaries",nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos,nbinspos,-maxpos,maxpos);
1796 //     fHistSecondaryPositionInDetectorMultiple->GetXaxis()->SetTitle("X");
1797 //     fHistSecondaryPositionInDetectorMultiple->GetYaxis()->SetTitle("Y");
1798 //     fHistSecondaryPositionInDetectorMultiple->GetZaxis()->SetTitle("Z");
1799     fClusterPositionWeird =  new TH2F("fClusterPositionWeird", "Position of weird secondary clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
1800
1801    fSecondaryClusterEnergy = new TH1F("fSecondaryClusterEnergy","fSecondaryClusterEnergy", 100, 0, 5);
1802 }
1803
1804 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
1805 {   //fill the output list
1806     AliAnalysisEt::FillOutputList(list);
1807
1808   if(!fIsMC) return;
1809     if(fCuts->GetHistMakeTree())
1810     {
1811         list->Add(fPrimaryTree);
1812     }
1813
1814     list->Add(fHistRemovedOrNot);
1815
1816     list->Add(fHistEtNonRemovedProtons);
1817     list->Add(fHistEtNonRemovedAntiProtons);
1818     list->Add(fHistEtNonRemovedPiPlus);
1819     list->Add(fHistEtNonRemovedPiMinus);
1820     list->Add(fHistEtNonRemovedKaonPlus);
1821     list->Add(fHistEtNonRemovedKaonMinus);
1822     list->Add(fHistEtNonRemovedK0s);
1823     list->Add(fHistEtNonRemovedK0L);
1824     list->Add(fHistEtNonRemovedLambdas);
1825     list->Add(fHistEtNonRemovedElectrons);
1826     list->Add(fHistEtNonRemovedPositrons);
1827     list->Add(fHistEtNonRemovedMuPlus);
1828     list->Add(fHistEtNonRemovedMuMinus);
1829     list->Add(fHistEtNonRemovedNeutrons);
1830     list->Add(fHistEtNonRemovedAntiNeutrons);
1831     list->Add(fHistEtNonRemovedGammas);
1832     list->Add(fHistEtNonRemovedGammasFromPi0);
1833
1834     list->Add(fHistEtRemovedGammas);
1835     list->Add(fHistEtRemovedNeutrons);
1836     list->Add(fHistEtRemovedAntiNeutrons);
1837
1838     list->Add(fHistEtRemovedCharged);
1839     list->Add(fHistEtRemovedNeutrals);
1840
1841     list->Add(fHistEtNonRemovedCharged);
1842     list->Add(fHistEtNonRemovedNeutrals);
1843
1844     list->Add(fHistMultNonRemovedProtons);
1845     list->Add(fHistMultNonRemovedAntiProtons);
1846     list->Add(fHistMultNonRemovedPiPlus);
1847     list->Add(fHistMultNonRemovedPiMinus);
1848     list->Add(fHistMultNonRemovedKaonPlus);
1849     list->Add(fHistMultNonRemovedKaonMinus);
1850     list->Add(fHistMultNonRemovedK0s);
1851     list->Add(fHistMultNonRemovedK0L);
1852     list->Add(fHistMultNonRemovedLambdas);
1853     list->Add(fHistMultNonRemovedElectrons);
1854     list->Add(fHistMultNonRemovedPositrons);
1855     list->Add(fHistMultNonRemovedMuPlus);
1856     list->Add(fHistMultNonRemovedMuMinus);
1857     list->Add(fHistMultNonRemovedNeutrons);
1858     list->Add(fHistMultNonRemovedAntiNeutrons);
1859     list->Add(fHistMultNonRemovedGammas);
1860
1861     list->Add(fHistMultRemovedGammas);
1862     list->Add(fHistMultRemovedNeutrons);
1863     list->Add(fHistMultRemovedAntiNeutrons);
1864
1865     list->Add(fHistMultRemovedCharged);
1866     list->Add(fHistMultRemovedNeutrals);
1867
1868     list->Add(fHistMultNonRemovedCharged);
1869     list->Add(fHistMultNonRemovedNeutrals);
1870
1871     list->Add(fHistTrackMultvsNonRemovedCharged);
1872     list->Add(fHistTrackMultvsNonRemovedNeutral);
1873     list->Add(fHistTrackMultvsRemovedGamma);
1874
1875     list->Add(fHistClusterMultvsNonRemovedCharged);
1876     list->Add(fHistClusterMultvsNonRemovedNeutral);
1877     list->Add(fHistClusterMultvsRemovedGamma);
1878
1879     //list->Add(fHistDecayVertexNonRemovedCharged);
1880     //list->Add(fHistDecayVertexNonRemovedNeutral);
1881     //list->Add(fHistDecayVertexRemovedCharged);
1882     //list->Add(fHistDecayVertexRemovedNeutral);
1883
1884     list->Add(fHistDxDzNonRemovedCharged);
1885     list->Add(fHistDxDzRemovedCharged);
1886     list->Add(fHistDxDzNonRemovedNeutral);
1887     list->Add(fHistDxDzRemovedNeutral);
1888
1889     if(fCalcForKaonCorrection){
1890       list->Add(fHistK0EDepositsVsPtInAcceptance);
1891       list->Add(fHistK0EGammaVsPtInAcceptance);
1892       list->Add(fHistK0EDepositsVsPtOutOfAcceptance);
1893       list->Add(fHistK0EGammaVsPtOutOfAcceptance);
1894       list->Add(fHistSimKaonsInAcceptance);
1895       list->Add(fHistSimK0SInAcceptance);
1896       list->Add(fHistSimK0LInAcceptance);
1897       list->Add(fHistSimKPlusInAcceptance);
1898       list->Add(fHistSimKMinusInAcceptance);
1899       list->Add(fHistSimKaonsOutOfAcceptance);
1900       list->Add(fHistSimKaonsInAcceptanceWithDepositsPrimaries);
1901       list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries);
1902       list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries);
1903     }
1904
1905     list->Add(fHistPiPlusMult);
1906     list->Add(fHistPiMinusMult);
1907     list->Add(fHistPiZeroMult);
1908     list->Add(fHistPiPlusMultAcc);
1909     list->Add(fHistPiMinusMultAcc);
1910     list->Add(fHistPiZeroMultAcc);
1911
1912     list->Add(fHistGammasFound);
1913     list->Add(fHistGammasGenerated);
1914     list->Add(fHistGammasFoundOutOfAccCent);
1915     list->Add(fHistGammasFoundCent);
1916     list->Add(fHistGammasGeneratedCent);
1917     list->Add(fHistChargedTracksCut);
1918     list->Add(fHistChargedTracksAccepted);
1919     list->Add(fHistGammasCut);
1920     list->Add(fHistGammasAccepted);
1921     if(fCalcTrackMatchVsMult){
1922       list->Add(fHistChargedTracksCutMult);
1923       list->Add(fHistChargedTracksAcceptedMult);
1924       list->Add(fHistChargedTracksAcceptedLowPtCent);
1925       list->Add(fHistChargedTracksAcceptedLowPtCent500MeV);
1926       list->Add(fHistChargedTracksAcceptedLowPtCentNoAntiProtons);
1927       list->Add(fHistChargedTracksAcceptedLowPtCentAntiProtons);
1928       list->Add(fHistGammasCutMult);
1929       list->Add(fHistGammasAcceptedMult);
1930     }
1931     list->Add(fHistBadTrackMatches);
1932     list->Add(fHistMatchedTracksEvspTBkgd);
1933     list->Add(fHistMatchedTracksEvspTSignal);
1934     if(fCalcTrackMatchVsMult){
1935       list->Add(fHistMatchedTracksEvspTBkgdPeripheral);
1936       list->Add(fHistMatchedTracksEvspTSignalPeripheral);
1937       list->Add(fHistMatchedTracksEvspTBkgdvsCent);
1938       list->Add(fHistMatchedTracksEvspTSignalvsCent);
1939       list->Add(fHistMatchedTracksEvspTBkgdvsCentEffCorr);
1940       list->Add(fHistMatchedTracksEvspTSignalvsCentEffCorr);
1941       list->Add(fHistChargedTracksCutPeripheral);
1942       list->Add(fHistChargedTracksAcceptedPeripheral);
1943       list->Add(fHistGammasCutPeripheral);
1944       list->Add(fHistGammasAcceptedPeripheral);
1945     }
1946     list->Add(fHistBadTrackMatchesdPhidEta);
1947     list->Add(fHistGoodTrackMatchesdPhidEta);
1948     list->Add(fHistHadronDepositsAll);
1949     list->Add(fHistHadronDepositsReco);
1950     list->Add(fHistHadronDepositsAllCent);
1951     list->Add(fHistHadronDepositsAllCent500MeV);
1952     list->Add(fHistHadronDepositsRecoCent);
1953     list->Add(fHistHadronsAllCent);
1954     list->Add(fHistMultChVsSignalVsMult);
1955     list->Add(fHistNeutralRemovedSecondaryEtVsCent);
1956     list->Add(fHistChargedRemovedSecondaryEtVsCent);
1957     list->Add(fHistNeutralNotRemovedSecondaryEtVsCent);
1958     list->Add(fHistChargedNotRemovedSecondaryEtVsCent);
1959     list->Add(fHistNeutralRemovedSecondaryNumVsNCluster);
1960     list->Add(fHistChargedRemovedSecondaryNumVsNCluster);
1961     list->Add(fHistNeutralNotRemovedSecondaryNumVsNCluster);
1962     list->Add(fHistChargedNotRemovedSecondaryNumVsNCluster);
1963     list->Add(fHistNeutralRemovedSecondaryNumVsCent);
1964     list->Add(fHistChargedRemovedSecondaryNumVsCent);
1965     list->Add(fHistNeutralNotRemovedSecondaryNumVsCent);
1966     list->Add(fHistChargedNotRemovedSecondaryNumVsCent);
1967     list->Add(fHistNeutronsEtVsCent);
1968     list->Add(fHistNeutronsNumVsCent);
1969     list->Add(fHistNotNeutronsNumVsCent);
1970     list->Add(fHistPiKPDepositedVsNch);
1971     list->Add(fHistPiKPNotTrackMatchedDepositedVsNch);
1972     list->Add(fHistNeutronsDepositedVsNch);
1973     list->Add(fHistAntiNeutronsDepositedVsNch);
1974     list->Add(fHistProtonsDepositedVsNch);
1975     list->Add(fHistAntiProtonsDepositedVsNch);
1976     list->Add(fHistProtonsNotTrackMatchedDepositedVsNch);
1977     list->Add(fHistAntiProtonsNotTrackMatchedDepositedVsNch);
1978     list->Add(fHistSecondariesVsNch);
1979     list->Add(fHistSecondariesVsNcl);
1980     list->Add(fHistSecondariesEffCorrVsNch);
1981     list->Add(fHistSecondariesEffCorrVsNch);
1982     list->Add(fHistSecondariesOutOfAccEffCorrVsNch);
1983     list->Add(fHistSecondariesDetectorCoverEffCorrVsNch);
1984     list->Add(fHistCentVsNchVsNcl);
1985     //list->Add(fHistSecondaryPositionInDetector);
1986     list->Add(fClusterPositionWeird);
1987     //list->Add(fHistSecondaryPositionInDetectorMultiple);
1988     list->Add(fSecondaryClusterEnergy);
1989
1990
1991 }
1992
1993
1994 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1995 {
1996     //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
1997     AliESDtrack *esdTrack = new AliESDtrack(part);
1998     // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1999
2000     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
2001
2002     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2003
2004     bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
2005     delete esdTrack;
2006
2007     return status;
2008 }
2009
2010 void AliAnalysisEtMonteCarlo::FillHistograms()
2011 {   // let base class fill its histograms, and us fill the local ones
2012     AliAnalysisEt::FillHistograms();
2013   if(!fIsMC) return;
2014     //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
2015
2016     fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
2017     fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
2018     fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
2019     fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
2020     fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
2021     fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
2022     fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
2023     fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
2024     fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
2025     fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
2026     fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
2027     fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
2028     fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
2029     fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
2030     fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
2031     fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
2032     fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
2033
2034     fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
2035     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
2036     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
2037
2038 //     fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
2039 //                                             +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
2040 //                                             +fEtRemovedProtons.
2041 //                              fCentClass);
2042 //     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
2043 //
2044 //     fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
2045 //                                             +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
2046 //                                             +fEtNonRemovedProtons,
2047 //                              fCentClass);
2048 //     fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
2049
2050     fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
2051     fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
2052     fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
2053     fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
2054
2055     fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
2056     fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
2057     fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
2058     fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
2059
2060
2061     fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
2062     fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
2063     fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
2064     fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
2065     fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
2066     fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
2067     fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
2068     fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
2069     fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
2070     fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
2071     fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
2072     fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
2073     fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
2074     fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
2075     fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
2076     fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
2077
2078     fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
2079     fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
2080     fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
2081
2082     fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
2083                                             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
2084                                             +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
2085                                             +fMultNonRemovedProtons);
2086
2087     fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
2088                                             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
2089
2090     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
2091                                        fMultRemovedGammas);
2092
2093     fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
2094             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
2095             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
2096             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
2097
2098     fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
2099             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
2100
2101     fHistClusterMultvsRemovedGamma->Fill(fNClusters,
2102                                          fMultRemovedGammas);
2103
2104 }
2105
2106
2107
2108
2109 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
2110 { // print family tree
2111     TParticle *part = stack->Particle(partIdx);
2112 //     if(part->GetPdgCode() == fgK0SCode)
2113     {
2114         std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
2115         std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
2116         std::cout << "Energy: " << part->Energy() << std::endl;
2117         Float_t vtx = TMath::Sqrt( TMath::Power(part->Vx(),2) + TMath::Power(part->Vy(),2) + TMath::Power(part->Vz(),2) );
2118         std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() <<"|Vtx| "<<vtx << std::endl;
2119     }
2120     return PrintMothers(partIdx, stack, 1);
2121 }
2122
2123 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
2124 { // print mothers
2125     char *tabs = new char[gen+1];
2126     for(Int_t i = 0; i < gen; ++i)
2127     {
2128         //std::cout << i << std::endl;
2129         tabs[i] = '\t';
2130     }
2131     tabs[gen] = '\0';
2132     Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2133     if(mothIdx < 0)
2134     {
2135       delete [] tabs;
2136       return 0;
2137     }
2138     TParticle *mother = stack->Particle(mothIdx);
2139 //     if(mother->GetPdgCode() == fgK0SCode)
2140     {
2141         //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
2142         std::cout << tabs << "Index: " << mothIdx << std::endl;
2143         std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
2144         std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
2145         std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
2146         if(mother->GetFirstMother() >= 0)
2147         {
2148           std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
2149           if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
2150           std::cout << std::endl;
2151         }
2152         Float_t vtx = TMath::Sqrt( TMath::Power(mother->Vx(),2) + TMath::Power(mother->Vy(),2) + TMath::Power(mother->Vz(),2) );
2153         std::cout<<tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() <<"|Vtx| "<<vtx << std::endl;
2154     }
2155     if(mother->GetPdgCode() == fgK0SCode)
2156     {
2157 //      std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
2158     }
2159 //  std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
2160 //  std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
2161 //  std::cout << "Energy: " << mother->Energy() << std::endl;
2162 //  std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
2163
2164     delete [] tabs;
2165     return PrintMothers(mothIdx, stack, gen+1) + 1;
2166 }
2167
2168 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
2169 { // get primary mother
2170     if(partIdx >= 0)
2171     {
2172         //return stack->GetPrimary(partIdx);
2173
2174         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2175         if(mothIdx < 0) return -1;
2176         TParticle *mother = stack->Particle(mothIdx);
2177         if(mother)
2178         {
2179             if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
2180             else return GetPrimMother(mothIdx, stack);
2181         }
2182         else
2183         {
2184             return -1;
2185         }
2186     }
2187     return -1;
2188 }
2189
2190 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
2191 { // get K0 in family
2192     if(partIdx >= 0)
2193     {
2194         if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
2195         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
2196         if(mothIdx < 0) return -1;
2197         TParticle *mother = stack->Particle(mothIdx);
2198         if(mother)
2199         {
2200           if(mother->GetPdgCode() == fgK0SCode)
2201             {
2202                 return mothIdx;
2203             }
2204             return GetK0InFamily(mothIdx, stack);
2205         }
2206         else
2207         {
2208             return -1;
2209         }
2210     }
2211     return -1;
2212 }
2213