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