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