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 "TVector3.h"
  33 using namespace std;
  34
  35 ClassImp(AliAnalysisEtMonteCarlo);
  36
  37
  38 // ctor
  39 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
  40     ,fImpactParameter(0)
  41     ,fNcoll(0)
  42     ,fNpart(0)
  43     ,fPrimaryTree(0)
  44     ,fTotEtWithSecondaryRemoved(0)
  45     ,fTotEtSecondaryFromEmEtPrimary(0)
  46     ,fTotEtSecondary(0)
  47     ,fPrimaryCode(0)
  48     ,fPrimaryCharge(0)
  49     ,fPrimaryE(0)
  50     ,fPrimaryEt(0)
  51     ,fPrimaryPx(0)
  52     ,fPrimaryPy(0)
  53     ,fPrimaryPz(0)
  54     ,fPrimaryVx(0)
  55     ,fPrimaryVy(0)
  56     ,fPrimaryVz(0)
  57     ,fPrimaryAccepted(0)
  58     ,fPrimaryMatched(0)
  59     ,fDepositedCode(0)
  60     ,fDepositedE(0)
  61     ,fDepositedEt(0)
  62     ,fDepositedCharge(0)
  63     ,fDepositedVx(0)
  64     ,fDepositedVy(0)
  65     ,fDepositedVz(0)
  66     ,fSecondary(kFALSE)
  67     ,fReconstructedE(0)
  68     ,fReconstructedEt(0)
  69     ,fTotPx(0)
  70     ,fTotPy(0)
  71     ,fTotPz(0)
  72     ,fClusterMult(0)
  73     ,fHistDecayVertexNonRemovedCharged(0)
  74     ,fHistDecayVertexRemovedCharged(0)
  75     ,fHistDecayVertexNonRemovedNeutral(0)
  76     ,fHistDecayVertexRemovedNeutral(0)
  77
  78     ,fHistRemovedOrNot(0)
  79     ,fHistEtNonRemovedProtons(0)
  80     ,fHistEtNonRemovedAntiProtons(0)
  81     ,fHistEtNonRemovedPiPlus(0)
  82     ,fHistEtNonRemovedPiMinus(0)
  83     ,fHistEtNonRemovedKaonPlus(0)
  84     ,fHistEtNonRemovedKaonMinus(0)
  85     ,fHistEtNonRemovedK0s(0)
  86     ,fHistEtNonRemovedK0L(0)
  87     ,fHistEtNonRemovedLambdas(0)
  88     ,fHistEtNonRemovedElectrons(0)
  89     ,fHistEtNonRemovedPositrons(0)
  90     ,fHistEtNonRemovedMuPlus(0)
  91     ,fHistEtNonRemovedMuMinus(0)
  92     ,fHistEtNonRemovedNeutrons(0)
  93     ,fHistEtNonRemovedAntiNeutrons(0)
  94     ,fHistEtNonRemovedGammas(0)
  95     ,fHistEtNonRemovedGammasFromPi0(0)
  96     ,fHistEtRemovedGammas(0)
  97     ,fHistEtRemovedNeutrons(0)
  98     ,fHistEtRemovedAntiNeutrons(0)
  99     ,fHistEtRemovedCharged(0)
 100     ,fHistEtRemovedNeutrals(0)
 101     ,fHistEtNonRemovedCharged(0)
 102     ,fHistEtNonRemovedNeutrals(0)
 103     ,fHistMultNonRemovedProtons(0)
 104     ,fHistMultNonRemovedAntiProtons(0)
 105     ,fHistMultNonRemovedPiPlus(0)
 106     ,fHistMultNonRemovedPiMinus(0)
 107     ,fHistMultNonRemovedKaonPlus(0)
 108     ,fHistMultNonRemovedKaonMinus(0)
 109     ,fHistMultNonRemovedK0s(0)
 110     ,fHistMultNonRemovedK0L(0)
 111     ,fHistMultNonRemovedLambdas(0)
 112     ,fHistMultNonRemovedElectrons(0)
 113     ,fHistMultNonRemovedPositrons(0)
 114     ,fHistMultNonRemovedMuPlus(0)
 115     ,fHistMultNonRemovedMuMinus(0)
 116     ,fHistMultNonRemovedNeutrons(0)
 117     ,fHistMultNonRemovedAntiNeutrons(0)
 118     ,fHistMultNonRemovedGammas(0)
 119     ,fHistMultRemovedGammas(0)
 120     ,fHistMultRemovedNeutrons(0)
 121     ,fHistMultRemovedAntiNeutrons(0)
 122     ,fHistMultRemovedCharged(0)
 123     ,fHistMultRemovedNeutrals(0)
 124     ,fHistMultNonRemovedCharged(0)
 125     ,fHistMultNonRemovedNeutrals(0)
 126     ,fHistTrackMultvsNonRemovedCharged(0)
 127     ,fHistTrackMultvsNonRemovedNeutral(0)
 128     ,fHistTrackMultvsRemovedGamma(0)
 129     ,fHistClusterMultvsNonRemovedCharged(0)
 130     ,fHistClusterMultvsNonRemovedNeutral(0)
 131     ,fHistClusterMultvsRemovedGamma(0)
 132     ,fHistMultvsNonRemovedChargedE(0)
 133     ,fHistMultvsNonRemovedNeutralE(0)
 134     ,fHistMultvsRemovedGammaE(0)
 135                                                   ,fHistK0EDepositsVsPtForEOver000MeV(0)
 136                                                   ,fHistK0EDepositsVsPtForEOver100MeV(0)
 137                                                   ,fHistK0EDepositsVsPtForEOver150MeV(0)
 138                                                   ,fHistK0EDepositsVsPtForEOver200MeV(0)
 139                                                   ,fHistK0EDepositsVsPtForEOver250MeV(0)
 140                                                   ,fHistK0EDepositsVsPtForEOver300MeV(0)
 141                                                   ,fHistK0EDepositsVsPtForEOver350MeV(0)
 142                                                   ,fHistK0EDepositsVsPtForEOver400MeV(0)
 143                                                   ,fHistK0EDepositsVsPtForEOver450MeV(0)
 144                                                   ,fHistK0EDepositsVsPtForEOver500MeV(0)
 145                                                   ,fHistK0EGammaVsPtForEOver000MeV(0)
 146                                                   ,fHistK0EGammaVsPtForEOver100MeV(0)
 147                                                   ,fHistK0EGammaVsPtForEOver150MeV(0)
 148                                                   ,fHistK0EGammaVsPtForEOver200MeV(0)
 149                                                   ,fHistK0EGammaVsPtForEOver250MeV(0)
 150                                                   ,fHistK0EGammaVsPtForEOver300MeV(0)
 151                                                   ,fHistK0EGammaVsPtForEOver350MeV(0)
 152                                                   ,fHistK0EGammaVsPtForEOver400MeV(0)
 153                                                   ,fHistK0EGammaVsPtForEOver450MeV(0)
 154                                                   ,fHistK0EGammaVsPtForEOver500MeV(0)
 155     ,fEtNonRemovedProtons(0)
 156     ,fEtNonRemovedAntiProtons(0)
 157     ,fEtNonRemovedPiPlus(0)
 158     ,fEtNonRemovedPiMinus(0)
 159     ,fEtNonRemovedKaonPlus(0)
 160     ,fEtNonRemovedKaonMinus(0)
 161     ,fEtNonRemovedK0S(0)
 162     ,fEtNonRemovedK0L(0)
 163     ,fEtNonRemovedLambdas(0)
 164     ,fEtNonRemovedElectrons(0)
 165     ,fEtNonRemovedPositrons(0)
 166     ,fEtNonRemovedMuMinus(0)
 167     ,fEtNonRemovedMuPlus(0)
 168     ,fEtNonRemovedGammas(0)
 169     ,fEtNonRemovedGammasFromPi0(0)
 170     ,fEtNonRemovedNeutrons(0)
 171     ,fEtNonRemovedAntiNeutrons(0)
 172     ,fEtRemovedProtons(0)
 173     ,fEtRemovedAntiProtons(0)
 174     ,fEtRemovedPiPlus(0)
 175     ,fEtRemovedPiMinus(0)
 176     ,fEtRemovedKaonPlus(0)
 177     ,fEtRemovedKaonMinus(0)
 178     ,fEtRemovedK0s(0)
 179     ,fEtRemovedK0L(0)
 180     ,fEtRemovedLambdas(0)
 181     ,fEtRemovedElectrons(0)
 182     ,fEtRemovedPositrons(0)
 183     ,fEtRemovedMuMinus(0)
 184     ,fEtRemovedMuPlus(0)
 185     ,fEtRemovedGammasFromPi0(0)
 186     ,fEtRemovedGammas(0)
 187     ,fEtRemovedNeutrons(0)
 188     ,fEtRemovedAntiNeutrons(0)
 189     ,fMultNonRemovedProtons(0)
 190     ,fMultNonRemovedAntiProtons(0)
 191     ,fMultNonRemovedPiPlus(0)
 192     ,fMultNonRemovedPiMinus(0)
 193     ,fMultNonRemovedKaonPlus(0)
 194     ,fMultNonRemovedKaonMinus(0)
 195     ,fMultNonRemovedK0s(0)
 196     ,fMultNonRemovedK0L(0)
 197     ,fMultNonRemovedLambdas(0)
 198     ,fMultNonRemovedElectrons(0)
 199     ,fMultNonRemovedPositrons(0)
 200     ,fMultNonRemovedMuMinus(0)
 201     ,fMultNonRemovedMuPlus(0)
 202     ,fMultNonRemovedGammas(0)
 203     ,fMultNonRemovedNeutrons(0)
 204     ,fMultNonRemovedAntiNeutrons(0)
 205     ,fMultRemovedProtons(0)
 206     ,fMultRemovedAntiProtons(0)
 207     ,fMultRemovedPiPlus(0)
 208     ,fMultRemovedPiMinus(0)
 209     ,fMultRemovedKaonPlus(0)
 210     ,fMultRemovedKaonMinus(0)
 211     ,fMultRemovedK0s(0)
 212     ,fMultRemovedK0L(0)
 213     ,fMultRemovedLambdas(0)
 214     ,fMultRemovedElectrons(0)
 215     ,fMultRemovedPositrons(0)
 216     ,fMultRemovedMuMinus(0)
 217     ,fMultRemovedMuPlus(0)
 218     ,fMultRemovedGammas(0)
 219     ,fMultRemovedNeutrons(0)
 220     ,fMultRemovedAntiNeutrons(0)
 221     ,fTrackMultInAcc(0)
 222     ,fHistDxDzNonRemovedCharged(0)
 223     ,fHistDxDzRemovedCharged(0)
 224     ,fHistDxDzNonRemovedNeutral(0)
 225     ,fHistDxDzRemovedNeutral(0)
 226     ,fHistPiPlusMult(0)
 227     ,fHistPiMinusMult(0)
 228     ,fHistPiZeroMult(0)
 229     ,fHistPiPlusMultAcc(0)
 230     ,fHistPiMinusMultAcc(0)
 231     ,fHistPiZeroMultAcc(0)
 232 //     ,fPiPlusMult(0)
 233 //     ,fPiMinusMult(0)
 234     ,fPiZeroMult(0)
 235     ,fPiPlusMultAcc(0)
 236     ,fPiMinusMultAcc(0)
 237     ,fPiZeroMultAcc(0)
 238     ,fNeutralRemoved(0)
 239     ,fChargedRemoved(0)
 240     ,fChargedNotRemoved(0)
 241     ,fNeutralNotRemoved(0)
 242     ,fGammaRemoved(0)
 243     ,fSecondaryNotRemoved(0)
 244     ,fEnergyNeutralRemoved(0)
 245     ,fEnergyChargedRemoved(0)
 246     ,fEnergyChargedNotRemoved(0)
 247     ,fEnergyNeutralNotRemoved(0)
 248     ,fEnergyGammaRemoved(0)
 249     ,fNClusters(0)
 250     ,fTotNeutralEtAfterMinEnergyCut(0)
 251     ,fHistGammasFound(0)
 252                                                   ,fHistGammasGenerated(0)
 253                                                   ,fHistChargedTracksCut(0)
 254                                                   ,fHistChargedTracksAccepted(0)
 255                                                   ,fHistGammasCut(0)
 256                                                   ,fHistGammasAccepted(0)
 257                                                   ,fHistBadTrackMatches(0)
 258                                                   ,fHistMatchedTracksEvspTBkgd(0)
 259                                                   ,fHistMatchedTracksEvspTSignal(0)
 260                                                   ,fHistMatchedTracksEvspTBkgdPeripheral(0)
 261                                                   ,fHistMatchedTracksEvspTSignalPeripheral(0)
 262                                                   ,fHistChargedTracksCutPeripheral(0)
 263                                                   ,fHistChargedTracksAcceptedPeripheral(0)
 264                                                   ,fHistGammasCutPeripheral(0)
 265                                                   ,fHistGammasAcceptedPeripheral(0)
 266                                                   ,fHistBadTrackMatchesdPhidEta(0)
 267                                                   ,fHistGoodTrackMatchesdPhidEta(0)
 268 {
 269 }
 270
 271 // dtor
 272 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
 273 {   //Destructor
 274
 275   if(fPrimaryTree){
 276     fPrimaryTree->Clear();
 277     delete fPrimaryTree;
 278   }
 279     delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
 280     delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
 281     delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
 282     delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
 283
 284     delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
 285
 286     delete fHistEtNonRemovedProtons; // enter comment here
 287     delete fHistEtNonRemovedAntiProtons; // enter comment here
 288     delete fHistEtNonRemovedPiPlus; // enter comment here
 289     delete fHistEtNonRemovedPiMinus; // enter comment here
 290     delete fHistEtNonRemovedKaonPlus; // enter comment here
 291     delete fHistEtNonRemovedKaonMinus; // enter comment here
 292     delete fHistEtNonRemovedK0s; // enter comment here
 293     delete fHistEtNonRemovedK0L; // enter comment here
 294     delete fHistEtNonRemovedLambdas; // enter comment here
 295     delete fHistEtNonRemovedElectrons; // enter comment here
 296     delete fHistEtNonRemovedPositrons; // enter comment here
 297     delete fHistEtNonRemovedMuPlus; // enter comment here
 298     delete fHistEtNonRemovedMuMinus; // enter comment here
 299     delete fHistEtNonRemovedNeutrons; // enter comment here
 300     delete fHistEtNonRemovedAntiNeutrons; // enter comment here
 301     delete fHistEtNonRemovedGammas; // enter comment here
 302     delete fHistEtNonRemovedGammasFromPi0; // enter comment here
 303
 304     delete fHistEtRemovedGammas; // enter comment here
 305     delete fHistEtRemovedNeutrons; // enter comment here
 306     delete fHistEtRemovedAntiNeutrons; // enter comment here
 307
 308
 309     delete fHistMultNonRemovedProtons; // enter comment here
 310     delete fHistMultNonRemovedAntiProtons; // enter comment here
 311     delete fHistMultNonRemovedPiPlus; // enter comment here
 312     delete fHistMultNonRemovedPiMinus; // enter comment here
 313     delete fHistMultNonRemovedKaonPlus; // enter comment here
 314     delete fHistMultNonRemovedKaonMinus; // enter comment here
 315     delete fHistMultNonRemovedK0s; // enter comment here
 316     delete fHistMultNonRemovedK0L; // enter comment here
 317     delete fHistMultNonRemovedLambdas; // enter comment here
 318     delete fHistMultNonRemovedElectrons; // enter comment here
 319     delete fHistMultNonRemovedPositrons; // enter comment here
 320     delete fHistMultNonRemovedMuPlus; // enter comment here
 321     delete fHistMultNonRemovedMuMinus; // enter comment here
 322     delete fHistMultNonRemovedNeutrons; // enter comment here
 323     delete fHistMultNonRemovedAntiNeutrons; // enter comment here
 324     delete fHistMultNonRemovedGammas; // enter comment here
 325
 326     delete fHistMultRemovedGammas; // enter comment here
 327     delete fHistMultRemovedNeutrons; // enter comment here
 328     delete fHistMultRemovedAntiNeutrons; // enter comment here
 329
 330     delete fHistTrackMultvsNonRemovedCharged; // enter comment here
 331     delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
 332     delete fHistTrackMultvsRemovedGamma; // enter comment here
 333
 334     delete fHistClusterMultvsNonRemovedCharged; // enter comment here
 335     delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
 336     delete fHistClusterMultvsRemovedGamma; // enter comment here
 337
 338     delete fHistMultvsNonRemovedChargedE; // enter comment here
 339     delete fHistMultvsNonRemovedNeutralE; // enter comment here
 340     delete fHistMultvsRemovedGammaE; // enter comment here
 341
 342     delete fHistK0EDepositsVsPtForEOver000MeV; // enter comment here
 343     delete fHistK0EDepositsVsPtForEOver100MeV; // enter comment here
 344     delete fHistK0EDepositsVsPtForEOver150MeV; // enter comment here
 345     delete fHistK0EDepositsVsPtForEOver200MeV; // enter comment here
 346     delete fHistK0EDepositsVsPtForEOver250MeV; // enter comment here
 347     delete fHistK0EDepositsVsPtForEOver300MeV; // enter comment here
 348     delete fHistK0EDepositsVsPtForEOver350MeV; // enter comment here
 349     delete fHistK0EDepositsVsPtForEOver400MeV; // enter comment here
 350     delete fHistK0EDepositsVsPtForEOver450MeV; // enter comment here
 351     delete fHistK0EDepositsVsPtForEOver500MeV; // enter comment here
 352     delete fHistK0EGammaVsPtForEOver000MeV; // enter comment here
 353     delete fHistK0EGammaVsPtForEOver100MeV; // enter comment here
 354     delete fHistK0EGammaVsPtForEOver150MeV; // enter comment here
 355     delete fHistK0EGammaVsPtForEOver200MeV; // enter comment here
 356     delete fHistK0EGammaVsPtForEOver250MeV; // enter comment here
 357     delete fHistK0EGammaVsPtForEOver300MeV; // enter comment here
 358     delete fHistK0EGammaVsPtForEOver350MeV; // enter comment here
 359     delete fHistK0EGammaVsPtForEOver400MeV; // enter comment here
 360     delete fHistK0EGammaVsPtForEOver450MeV; // enter comment here
 361     delete fHistK0EGammaVsPtForEOver500MeV; // enter comment here
 362
 363     delete fHistDxDzNonRemovedCharged; // enter comment here
 364     delete fHistDxDzRemovedCharged; // enter comment here
 365     delete fHistDxDzNonRemovedNeutral; // enter comment here
 366     delete fHistDxDzRemovedNeutral; // enter comment here
 367
 368     delete fHistPiPlusMult; // enter comment here
 369     delete fHistPiMinusMult; // enter comment here
 370     delete fHistPiZeroMult; // enter comment here
 371
 372     delete fHistPiPlusMultAcc; // enter comment here
 373     delete fHistPiMinusMultAcc; // enter comment here
 374     delete fHistPiZeroMultAcc; // enter comment here
 375     delete fHistGammasFound; // enter comment here
 376     delete fHistGammasGenerated; // enter comment here
 377     delete fHistChargedTracksCut;
 378     delete fHistChargedTracksAccepted;
 379     delete fHistGammasCut;
 380     delete fHistGammasAccepted;
 381     delete fHistBadTrackMatches;
 382     delete fHistMatchedTracksEvspTBkgd;
 383     delete fHistMatchedTracksEvspTSignal;
 384     delete fHistMatchedTracksEvspTBkgdPeripheral;
 385     delete fHistMatchedTracksEvspTSignalPeripheral;
 386     delete fHistChargedTracksCutPeripheral;
 387     delete fHistChargedTracksAcceptedPeripheral;
 388     delete fHistGammasCutPeripheral;
 389     delete fHistGammasAcceptedPeripheral;
 390     delete fHistBadTrackMatchesdPhidEta;
 391     delete fHistGoodTrackMatchesdPhidEta;
 392 }
 393
 394 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 395 {   // analyse MC event
 396     ResetEventValues();
 397
 398
 399     fPiPlusMult = 0;
 400     fPiMinusMult = 0;
 401     fPiZeroMult = 0;
 402     if (fCentrality)
 403     {
 404         fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
 405
 406     }
 407
 408     // Get us an mc event
 409     if (!ev) {
 410         AliFatal("ERROR: Event does not exist");
 411         return 0;
 412     }
 413     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
 414     if (!event) {
 415         AliFatal("ERROR: MC Event does not exist");
 416         return 0;
 417     }
 418
 419     Double_t protonMass =fgProtonMass;
 420
 421     // Hijing header
 422     AliGenEventHeader* genHeader = event->GenEventHeader();
 423     if (!genHeader) {
 424         Printf("ERROR: Event generation header does not exist");
 425         return 0;
 426     }
 427     AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 428     if (hijingGenHeader) {
 429         fImpactParameter = hijingGenHeader->ImpactParameter();
 430         fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
 431         fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
 432     }
 433
 434     // Let's play with the stack!
 435     AliStack *stack = event->Stack();
 436
 437     Int_t nPrim = stack->GetNtrack();
 438
 439     Int_t partCount = 0;
 440     for (Int_t iPart = 0; iPart < nPrim; iPart++)
 441     {
 442
 443         TParticle *part = stack->Particle(iPart);
 444
 445
 446
 447         if (!part)
 448         {
 449             Printf("ERROR: Could not get particle %d", iPart);
 450             continue;
 451         }
 452         TParticlePDG *pdg = part->GetPDG(0);
 453
 454         if (!pdg)
 455         {
 456             Printf("ERROR: Could not get particle PDG %d", iPart);
 457             continue;
 458         }
 459
 460         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
 461         Int_t code = pdg->PdgCode();
 462
 463         if(stack->IsPhysicalPrimary(iPart))
 464         {
 465           fTotPx += part->Px();
 466           fTotPy += part->Py();
 467           fTotPz += part->Pz();
 468         }
 469
 470         // Check for reasonable (for now neutral and singly charged) charge on the particle
 471         if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
 472         {
 473
 474             fMultiplicity++;
 475 //             PrintFamilyTree(iPart, stack);
 476 //
 477             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 478             {
 479                 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
 480                 // calculate E_T
 481                 if (
 482                     TMath::Abs(code) == fgProtonCode ||
 483                     TMath::Abs(code) == fgNeutronCode ||
 484                     TMath::Abs(code) == fgLambdaCode ||
 485                     TMath::Abs(code) == fgXiCode ||
 486                     TMath::Abs(code) == fgXi0Code ||
 487                     TMath::Abs(code) == fgOmegaCode
 488                 )
 489                 {
 490                     if (code > 0) {
 491                         particleMassPart = - protonMass;
 492                     }
 493                     if (code < 0) {
 494                         particleMassPart = protonMass;
 495                     }
 496                 }
 497                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 498
 499
 500                 // Fill up total E_T counters for each particle species
 501                 if (code == fgProtonCode || code == fgAntiProtonCode)
 502                 {
 503                 }
 504                 if (code == fgPiPlusCode || code == fgPiMinusCode)
 505                 {
 506                     if (code == fgPiPlusCode)
 507                     {
 508                     }
 509                     else
 510                     {
 511                     }
 512                 }
 513                 if (code == fgGammaCode)
 514                 {
 515                 }
 516                 if (code == fgKPlusCode)
 517                 {
 518                 }
 519                 if(code == fgKMinusCode)
 520                 {
 521                 }
 522                 if (code == fgMuPlusCode || code == fgMuMinusCode)
 523                 {
 524                 }
 525                 if (code == fgEPlusCode || code == fgEMinusCode)
 526                 {
 527                 }
 528                 // some neutrals also
 529                 if (code == fgNeutronCode)
 530                 {
 531                 }
 532                 if (code == fgAntiNeutronCode)
 533                 {
 534                 }
 535                 if (code == fgGammaCode)
 536                 {
 537                 }
 538
 539                 // Neutral particles
 540                 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
 541
 542                 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
 543                 {
 544                 //  PrintFamilyTree(iPart,stack);
 545                     //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
 546                     //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
 547
 548                     // inside EMCal acceptance
 549
 550                     //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
 551
 552                     if(fSelector->CutGeometricalAcceptance(*part) )
 553                     {
 554                         fNeutralMultiplicity++;
 555                         fTotNeutralEt += et;
 556                         if(fSelector->PassMinEnergyCut(*part))
 557                         {
 558                             fTotNeutralEtAfterMinEnergyCut += et;
 559                         }
 560                         if (part->Energy() > 0.05) partCount++;
 561                     }
 562                 }
 563                 //Charged particles
 564                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
 565                 {
 566
 567                     // inside EMCal acceptance
 568                     if (fSelector->CutGeometricalAcceptance(*part))
 569                     {
 570
 571                         fChargedMultiplicity++;
 572
 573                         fTotChargedEt += et;
 574
 575                         if (code == fgProtonCode || code == fgAntiProtonCode)
 576                         {
 577                         }
 578                         if (code == fgPiPlusCode || code == fgPiMinusCode)
 579                         {
 580                         }
 581                         if (code == fgKPlusCode || code == fgKMinusCode)
 582                         {
 583                         }
 584                         if (code == fgMuPlusCode || code == fgMuMinusCode)
 585                         {
 586                         }
 587
 588                         if (code == fgEPlusCode || code == fgEMinusCode)
 589                         {
 590                             fTotNeutralEt += et; // calling electrons neutral
 591                             fTotChargedEt -= et;
 592                         }
 593                     } // inside EMCal acceptance
 594
 595                     if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 596                     {
 597                         if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 598                         else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 599                     }
 600                 }
 601             }
 602
 603         }
 604     }
 605    // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
 606     fTotEt = fTotChargedEt + fTotNeutralEt;
 607     //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
 608     //std::cout << "Event done! # of particles: " << partCount << std::endl;
 609     return 0;
 610 }
 611 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 612 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 613 {   // analyse MC and real event info
 614     //if(!mcEvent || !realEvent){
 615     if (!ev || !ev2) {
 616         AliFatal("ERROR: Event does not exist");
 617         return 0;
 618     }
 619     AliAnalysisEt::AnalyseEvent(ev);
 620     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 621     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
 622     if (!mcEvent || !realEvent) {
 623         AliFatal("ERROR: mcEvent or realEvent does not exist");
 624
 625     }
 626
 627     std::vector<Int_t> foundGammas;
 628
 629     fSelector->SetEvent(realEvent);
 630
 631     AnalyseEvent(ev);
 632
 633     AliStack *stack = mcEvent->Stack();
 634
 635     // get all detector clusters
 636     //  TRefArray* caloClusters = new TRefArray();
 637
 638 //    if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
 639     //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
 640     //else {
 641     //AliFatal("Detector ID has not been specified");
 642     //return -1;
 643 //    }
 644
 645 //Note that this only returns clusters for the selected detector.  fSelector actually calls the right GetClusters... for the detector
 646 //It does not apply any cuts on these clusters
 647     TRefArray *caloClusters = fSelector->GetClusters();
 648
 649     Int_t nCluster = caloClusters->GetEntries();
 650     fClusterMult = nCluster;
 651     fNClusters = 0;
 652     // loop the clusters
 653     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 654     {
 655         Int_t cf = 0;
 656         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 657         //Float_t caloE = caloCluster->E()
 658         fNClusters++;
 659         const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
 660         TParticle *part  =  stack->Particle(iPart);
 661
 662         if (!part)
 663         {
 664             Printf("No MC particle %d", iCluster);
 665             continue;
 666         }
 667
 668         int primIdx = iPart;
 669         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
 670         {
 671             primIdx = GetPrimMother(iPart, stack);
 672             if(primIdx != stack->GetPrimary(iPart))
 673             {
 674               //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
 675               //PrintFamilyTree(iPart, stack);
 676             }
 677             //if it is from a K0S
 678             if(primIdx < 0)
 679             {
 680               std::cout << "What!? No primary?" << std::endl;
 681               PrintFamilyTree(iPart, stack);
 682               //continue;
 683               //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.
 684               primIdx = iPart;
 685             }
 686
 687         } // end of primary particle check
 688
 689         //const int primCode = stack->Particle(primIdx)->GetPdgCode();
 690         TParticlePDG *pdg = part->GetPDG();
 691         if (!pdg)
 692         {
 693             Printf("ERROR: Could not get particle PDG %d", iPart);
 694             continue;
 695         }
 696
 697         //Int_t code = pdg->PdgCode();
 698 //      if(primCode == fgGammaCode)
 699 //      {
 700
 701         Float_t pos[3];
 702         caloCluster->GetPosition(pos);
 703         TVector3 cp(pos);
 704         for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
 705         {
 706           Int_t pIdx = caloCluster->GetLabelAt(i);
 707
 708
 709           TParticle *p = stack->Particle(pIdx);
 710           if(!stack->IsPhysicalPrimary(pIdx))
 711           {
 712
 713             TVector3 decayVtx(p->Vx(), p->Vy(), p->Vz());
 714              if((decayVtx-cp).Mag() < 30.0)
 715               {
 716                 continue;
 717               }
 718             pIdx = GetPrimMother(pIdx, stack);
 719           }
 720           if((stack->Particle(pIdx)->GetPdgCode() == fgGammaCode))
 721           {
 722                if(fSelector->PassDistanceToBadChannelCut(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
 723                 {
 724                  // std::cout << "Gamma primary: " << primIdx << std::endl;
 725                   foundGammas.push_back(pIdx);
 726                 }
 727           }
 728         }
 729         fCutFlow->Fill(cf++);
 730         if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
 731         Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster);
 732 //      if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
 733         if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
 734
 735
 736         fCutFlow->Fill(cf++);
 737
 738         TParticle *primPart = stack->Particle(primIdx);
 739         fPrimaryCode = primPart->GetPdgCode();
 740         fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
 741
 742         fPrimaryE = primPart->Energy();
 743         fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
 744         fPrimaryPx = primPart->Px();
 745         fPrimaryPy = primPart->Py();
 746         fPrimaryPz = primPart->Pz();
 747         //cout<<"I have a cluster and it's good energy "<<caloCluster->E()<<" simulated "<<fPrimaryE<<endl;
 748         fPrimaryVx = primPart->Vx();
 749         fPrimaryVy = primPart->Vy();
 750         fPrimaryVz = primPart->Vz();
 751
 752         fPrimaryAccepted = false;
 753         fPrimaryMatched = false;
 754
 755         fDepositedCode = part->GetPdgCode();
 756         fDepositedE = part->Energy();
 757         fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
 758         fDepositedCharge = (Int_t) part->GetPDG()->Charge();
 759
 760         fDepositedVx = part->Vx();
 761         fDepositedVy = part->Vy();
 762         fDepositedVz = part->Vz();
 763
 764         //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
 765         fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
 766 //      if(fSecondary)
 767 //      {
 768 //        //std::cout << "Have secondary!" << std::endl;
 769 //        //PrintFamilyTree(iPart, stack);
 770 //      }
 771         fReconstructedE = caloCluster->E();
 772         fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
 773
 774         pdg = primPart->GetPDG(0);
 775         //Int_t code = primPart->GetPdgCode();
 776
 777         Bool_t written = kFALSE;
 778
 779         Bool_t nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
 780
 781         if(fSecondary){//all particles from secondary interactions
 782           written = kTRUE;
 783           if(nottrackmatched){//secondaries not removed
 784             if (fDepositedCharge != 0){//charged track not removed
 785               fChargedNotRemoved++;
 786               fEnergyChargedNotRemoved += clEt;
 787               fHistRemovedOrNot->Fill(2.0, fCentClass);
 788               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 789             }
 790             else{
 791               fSecondaryNotRemoved++;
 792             }
 793           }
 794           else{//secondaries removed
 795             if (fDepositedCharge != 0){
 796               fChargedRemoved++;
 797               fEnergyChargedRemoved += clEt;
 798               fHistRemovedOrNot->Fill(0.0, fCentClass);
 799               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 800               fHistChargedTracksCut->Fill(fDepositedEt);
 801               if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
 802               fHistMatchedTracksEvspTBkgd->Fill(part->Pt(),fReconstructedE);
 803               if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(part->P(),fReconstructedE);}
 804               Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
 805               if(caloCluster->GetLabel()!=trackindex){
 806                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
 807                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
 808                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
 809               }
 810               else{
 811                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
 812               }
 813             }
 814             else{//neutral energy removed
 815               fNeutralRemoved++;
 816               fEnergyNeutralRemoved += clEt;
 817               fHistRemovedOrNot->Fill(1.0, fCentClass);
 818               fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 819             }
 820           }
 821         }
 822         else{
 823
 824           if (fDepositedCharge != 0 && fDepositedCode!=fgEMinusCode && fDepositedCode!=fgEPlusCode){//if the particle hitting the calorimeter is pi/k/p/mu
 825             written = kTRUE;
 826             if(nottrackmatched){//not removed but should be
 827               fChargedNotRemoved++;
 828               fEnergyChargedNotRemoved += clEt;
 829               fHistRemovedOrNot->Fill(2.0, fCentClass);
 830               fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 831               fHistChargedTracksAccepted->Fill(fDepositedEt);
 832               if(fClusterMult<25){fHistChargedTracksAcceptedPeripheral->Fill(fDepositedEt);}
 833             }
 834             else{//removed and should have been
 835               Int_t trackindex = (caloCluster->GetLabelsArray())->At(0);
 836               if(caloCluster->GetLabel()!=trackindex){
 837                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
 838                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
 839                 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
 840               }
 841               else{
 842                 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
 843               }
 844               fChargedRemoved++;
 845               fEnergyChargedRemoved += clEt;
 846               fHistRemovedOrNot->Fill(0.0, fCentClass);
 847               fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 848               fHistChargedTracksCut->Fill(fDepositedEt);
 849               if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
 850               fHistMatchedTracksEvspTBkgd->Fill(part->P(),fReconstructedE);
 851               if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(part->P(),fReconstructedE);}
 852             }
 853           }
 854           //K0L and any neutral particles from the decay of K+/- or K0S
 855           if(!written && (fPrimaryCode==fgKPlusCode || fPrimaryCode==fgKMinusCode || fPrimaryCode==fgK0SCode ||fPrimaryCode==fgK0LCode)){
 856             written = kTRUE;//At this point we are not tracking them but we don't count them as neutrals accidentally removed
 857           }
 858
 859           if(!written && (fDepositedCode==fgGammaCode || fDepositedCode==fgEMinusCode || fDepositedCode ==fgEPlusCode)){//if the particle hitting the calorimeter is gamma, electron and not from a kaon
 860             written = kTRUE;
 861             if(nottrackmatched){//Not removed and not supposed to be removed - signal
 862               fEtNonRemovedGammas += clEt;
 863               fMultNonRemovedGammas++;
 864               fNeutralNotRemoved--;
 865               fEnergyNeutralNotRemoved -= clEt;
 866               fHistGammasAccepted->Fill(fDepositedEt);
 867               if(fClusterMult<25){fHistGammasAcceptedPeripheral->Fill(fDepositedEt);}
 868             }
 869             else{//removed but shouldn't have been
 870               Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
 871               if(caloCluster->GetLabel()!=trackindex){
 872                 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
 873                 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
 874 //              cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
 875 //              PrintFamilyTree(trackindex, stack);
 876 //              cout<<"Cluster"<<endl;
 877               }
 878               fGammaRemoved++;
 879               fGammaRemovedEt+=clEt;
 880               fHistGammasCut->Fill(fDepositedEt);
 881               if(fClusterMult<25){fHistGammasCutPeripheral->Fill(fDepositedEt);}
 882               fHistMatchedTracksEvspTSignal->Fill(part->P(),fReconstructedE);
 883               if(fClusterMult<25){fHistMatchedTracksEvspTSignalPeripheral->Fill(part->P(),fReconstructedE);}
 884             }
 885           }
 886           //all other cases - neutron, anti-neutron, not aware of other cases
 887           if(!written){
 888             fNeutralNotRemoved++;
 889             fEnergyNeutralNotRemoved += clEt;
 890             fHistRemovedOrNot->Fill(3.0, fCentClass);
 891             fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 892           }
 893         }
 894         fPrimaryTree->Fill();
 895     } // end of loop over clusters
 896
 897
 898     std::sort(foundGammas.begin(), foundGammas.end());
 899     for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
 900     {
 901
 902         if(!stack->IsPhysicalPrimary(iPart)) continue;
 903
 904         TParticle *part = stack->Particle(iPart);
 905
 906         if (!part)
 907         {
 908             Printf("ERROR: Could not get particle %d", iPart);
 909             continue;
 910         }
 911         TParticlePDG *pdg = part->GetPDG(0);
 912
 913         if (!pdg)
 914         {
 915             Printf("ERROR: Could not get particle PDG %d", iPart);
 916             continue;
 917         }
 918
 919         if(pdg->PdgCode()==fgGammaCode && fSelector->CutGeometricalAcceptance(*part))// TMath::Abs(part->Eta()) < 0.12)
 920         {
 921           fHistGammasGenerated->Fill(part->Energy());
 922           if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
 923           {
 924             fHistGammasFound->Fill(part->Energy());
 925           }
 926         }
 927
 928     }
 929     Float_t etCuts[10] = {0.0,0.1,0.15,0.2,0.25, 0.3,0.35,0.4,0.45,0.5};
 930     Int_t nEtCuts = 10;
 931     //loop over simulated particles in order to find K0S
 932     for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++){
 933       TParticle *part = stack->Particle(iPart);
 934       if (!part){
 935         //Printf("ERROR: Could not get particle %d", iPart);
 936         continue;
 937       }
 938       TParticlePDG *pdg = part->GetPDG(0);
 939       if (!pdg){
 940         //Printf("ERROR: Could not get particle PDG %d", iPart);
 941         continue;
 942       }
 943
 944       if(stack->IsPhysicalPrimary(iPart)){//if it is a K0 it might have decayed into four pions
 945         //fgK0SCode, fgGammaCode, fgPi0Code
 946         Int_t code = pdg->PdgCode();
 947         if(code == fgK0SCode){//this is a K0S
 948           if(stack->Particle(iPart)->GetNDaughters()==2){//it has two daughters
 949             Int_t daughterID = stack->Particle(iPart)->GetFirstDaughter();
 950             if(stack->Particle(daughterID)->GetPDG(0)->PdgCode()==fgPi0Code){//and it decayed into a pi0
 951               //cout<<"I am a gamma from a pi0 from a K0S"<<endl;
 952               Float_t pTk0 = stack->Particle(iPart)->Pt();
 953               Int_t gammaDaughterIDs[6] = {0,0,0,0,0,0};
 954               Float_t gammaEts[4] = {0.0,0.0,0.0,0.0};
 955               Float_t totalGammaEts[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 956               Float_t totalClusterEts[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 957               gammaDaughterIDs[4] = stack->Particle(iPart)->GetDaughter(0);
 958               gammaDaughterIDs[5] = stack->Particle(iPart)->GetDaughter(1);
 959               if(gammaDaughterIDs[4]>0){
 960                 gammaDaughterIDs[0] = stack->Particle(gammaDaughterIDs[4])->GetDaughter(0);
 961                 gammaDaughterIDs[1] = stack->Particle(gammaDaughterIDs[4])->GetDaughter(1);
 962               }
 963               if(gammaDaughterIDs[5]>0){
 964                 gammaDaughterIDs[2] = stack->Particle(gammaDaughterIDs[5])->GetDaughter(0);
 965                 gammaDaughterIDs[3] = stack->Particle(gammaDaughterIDs[5])->GetDaughter(1);
 966               }
 967               for(int k=0;k<4;k++){
 968                 //if( TMath::Abs(stack->Particle(gammaDaughterIDs[k])->Eta()) <= fCuts->GetGeometryEmcalEtaAccCut() ){//only add the energy if it's within the detector acceptance
 969                 if(gammaDaughterIDs[k]==-1) continue;
 970                 if( fSelector->CutGeometricalAcceptance( * stack->Particle(gammaDaughterIDs[k]))){//only add the energy if it's within the detector acceptance
 971                   //cout<<"Found a gamma "<<" K0S eta "<<stack->Particle(iPart)->Eta()<<" gamma daughter  eta "<< stack->Particle(gammaDaughterIDs[k])->Eta()<<endl;
 972                   gammaEts[k] = stack->Particle(gammaDaughterIDs[4])->Energy() * TMath::Sin(stack->Particle(gammaDaughterIDs[4])->Theta() );
 973                 }
 974                 else{
 975                   //cout<<"Eta rejected "<<" K0S eta "<<stack->Particle(iPart)->Eta()<<" gamma daughter  eta "<< stack->Particle(gammaDaughterIDs[k])->Eta()<<" eta cut "<<fCuts->GetGeometryEmcalEtaAccCut()<<endl;
 976                 }
 977               }
 978               //does not always have all of the daughters
 979               if(gammaDaughterIDs[0]==gammaDaughterIDs[1]){
 980                 gammaDaughterIDs[1] = -1;
 981                 gammaEts[1] = 0.0;
 982                 //cout<<"Duplicates.  This has "<<stack->Particle(gammaDaughterIDs[4])->GetNDaughters()<<" daughters"<<endl;
 983               }
 984               if(gammaDaughterIDs[2]==gammaDaughterIDs[3]){
 985                 gammaDaughterIDs[3] = -1;
 986                 gammaEts[3] = 0.0;
 987                 //cout<<"Duplicates.  This has "<<stack->Particle(gammaDaughterIDs[5])->GetNDaughters()<<" daughters"<<endl;
 988               }
 989               for(int l=0;l<nEtCuts;l++){//loop over cut values
 990                 for(int k=0;k<4;k++){//loop over gamma daughter energies
 991                   if(gammaEts[k]>=etCuts[l]){
 992                     totalGammaEts[l] += gammaEts[k];//if gamma et is above the cut off energy add it
 993                   }
 994                 }
 995               }
 996
 997               for (int iCluster = 0; iCluster < nCluster; iCluster++ ){//if this cluster is from any of the decay daughters of our K0S
 998                 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 999                 const UInt_t myPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
1000                 for(int j=0;j<6;j++){
1001                   if( myPart==((UInt_t) gammaDaughterIDs[j]) ){
1002                     //Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster);
1003                     Double_t clEt = caloCluster->E();
1004                     //cout<<"Found a matching cluster!!!  Energy: "<<clEt<<endl;
1005                     //cout<<" it has energy "<<clEt;
1006                     for(int l=0;l<nEtCuts;l++){//loop over cut values
1007                       //cout<<", cut = "<<etCuts[l];
1008                       if(clEt>=etCuts[l]){
1009                         //cout<<", "<<clEt<<">="<<etCuts[l];
1010                         totalClusterEts[l] += clEt;//if cluster et is above the cut off energy add it
1011                       }
1012                     }
1013                     //cout<<endl;
1014
1015                   }
1016                 }
1017               }
1018               //now we have the total energy measured for each case
1019               fHistK0EDepositsVsPtForEOver000MeV->Fill(pTk0,totalClusterEts[0]);
1020               fHistK0EDepositsVsPtForEOver100MeV->Fill(pTk0,totalClusterEts[1]);
1021               fHistK0EDepositsVsPtForEOver150MeV->Fill(pTk0,totalClusterEts[2]);
1022               fHistK0EDepositsVsPtForEOver200MeV->Fill(pTk0,totalClusterEts[3]);
1023               fHistK0EDepositsVsPtForEOver250MeV->Fill(pTk0,totalClusterEts[4]);
1024               fHistK0EDepositsVsPtForEOver300MeV->Fill(pTk0,totalClusterEts[5]);
1025               fHistK0EDepositsVsPtForEOver350MeV->Fill(pTk0,totalClusterEts[6]);
1026               fHistK0EDepositsVsPtForEOver400MeV->Fill(pTk0,totalClusterEts[7]);
1027               fHistK0EDepositsVsPtForEOver450MeV->Fill(pTk0,totalClusterEts[8]);
1028               fHistK0EDepositsVsPtForEOver500MeV->Fill(pTk0,totalClusterEts[9]);
1029               fHistK0EGammaVsPtForEOver000MeV->Fill(pTk0,totalGammaEts[0]);
1030               fHistK0EGammaVsPtForEOver100MeV->Fill(pTk0,totalGammaEts[1]);
1031               fHistK0EGammaVsPtForEOver150MeV->Fill(pTk0,totalGammaEts[2]);
1032               fHistK0EGammaVsPtForEOver200MeV->Fill(pTk0,totalGammaEts[3]);
1033               fHistK0EGammaVsPtForEOver250MeV->Fill(pTk0,totalGammaEts[4]);
1034               fHistK0EGammaVsPtForEOver300MeV->Fill(pTk0,totalGammaEts[5]);
1035               fHistK0EGammaVsPtForEOver350MeV->Fill(pTk0,totalGammaEts[6]);
1036               fHistK0EGammaVsPtForEOver400MeV->Fill(pTk0,totalGammaEts[7]);
1037               fHistK0EGammaVsPtForEOver450MeV->Fill(pTk0,totalGammaEts[8]);
1038               fHistK0EGammaVsPtForEOver500MeV->Fill(pTk0,totalGammaEts[9]);
1039             }
1040           }
1041         }
1042       }
1043     }
1044
1045     FillHistograms();
1046     return 0;
1047 }
1048
1049 void AliAnalysisEtMonteCarlo::Init()
1050 {   // init
1051     AliAnalysisEt::Init();
1052 }
1053
1054 void AliAnalysisEtMonteCarlo::ResetEventValues()
1055 {   // reset event values
1056     AliAnalysisEt::ResetEventValues();
1057
1058     fTotEtSecondary = 0;
1059     fTotEtSecondaryFromEmEtPrimary = 0;
1060     fTotEtWithSecondaryRemoved = 0;
1061
1062     // collision geometry defaults for p+p:
1063     fImpactParameter = 0;
1064     fNcoll = 1;
1065     fNpart = 2;
1066
1067     fEtNonRemovedProtons = 0;
1068     fEtNonRemovedAntiProtons = 0;
1069     fEtNonRemovedPiPlus = 0;
1070     fEtNonRemovedPiMinus = 0;
1071     fEtNonRemovedKaonPlus = 0;
1072     fEtNonRemovedKaonMinus = 0;
1073     fEtNonRemovedK0S = 0;
1074     fEtNonRemovedK0L = 0;
1075     fEtNonRemovedLambdas = 0;
1076     fEtNonRemovedElectrons = 0;
1077     fEtNonRemovedPositrons = 0;
1078     fEtNonRemovedMuPlus = 0;
1079     fEtNonRemovedMuMinus = 0;
1080     fEtNonRemovedNeutrons = 0;
1081     fEtNonRemovedAntiNeutrons = 0;
1082     fEtNonRemovedGammas = 0;
1083     fEtNonRemovedGammasFromPi0 = 0;
1084
1085     fEtRemovedProtons = 0;
1086     fEtRemovedAntiProtons = 0;
1087     fEtRemovedPiPlus = 0;
1088     fEtRemovedPiMinus = 0;
1089     fEtRemovedKaonPlus = 0;
1090     fEtRemovedKaonMinus = 0;
1091     fEtRemovedK0s = 0;
1092     fEtRemovedK0L = 0;
1093     fEtRemovedLambdas = 0;
1094     fEtRemovedElectrons = 0;
1095     fEtRemovedPositrons = 0;
1096     fEtRemovedMuPlus = 0;
1097     fEtRemovedMuMinus = 0;
1098     fEtRemovedNeutrons = 0;
1099
1100     fEtRemovedGammasFromPi0 = 0;
1101     fEtRemovedGammas = 0;
1102     fEtRemovedNeutrons = 0;
1103     fEtRemovedAntiNeutrons = 0;
1104
1105     fMultNonRemovedProtons = 0;
1106     fMultNonRemovedAntiProtons = 0;
1107     fMultNonRemovedPiPlus = 0;
1108     fMultNonRemovedPiMinus = 0;
1109     fMultNonRemovedKaonPlus = 0;
1110     fMultNonRemovedKaonMinus = 0;
1111     fMultNonRemovedK0s = 0;
1112     fMultNonRemovedK0L = 0;
1113     fMultNonRemovedLambdas = 0;
1114     fMultNonRemovedElectrons = 0;
1115     fMultNonRemovedPositrons = 0;
1116     fMultNonRemovedMuPlus = 0;
1117     fMultNonRemovedMuMinus = 0;
1118     fMultNonRemovedNeutrons = 0;
1119     fMultNonRemovedAntiNeutrons = 0;
1120     fMultNonRemovedGammas = 0;
1121
1122     fMultRemovedProtons = 0;
1123     fMultRemovedAntiProtons = 0;
1124     fMultRemovedPiPlus = 0;
1125     fMultRemovedPiMinus = 0;
1126     fMultRemovedKaonPlus = 0;
1127     fMultRemovedKaonMinus = 0;
1128     fMultRemovedK0s = 0;
1129     fMultRemovedK0L = 0;
1130     fMultRemovedLambdas = 0;
1131     fMultRemovedElectrons = 0;
1132     fMultRemovedPositrons = 0;
1133     fMultRemovedMuPlus = 0;
1134     fMultRemovedMuMinus = 0;
1135
1136     fMultRemovedGammas = 0;
1137     fMultRemovedNeutrons = 0;
1138     fMultRemovedAntiNeutrons = 0;
1139
1140     fEnergyChargedNotRemoved = 0;
1141     fEnergyChargedRemoved = 0;
1142     fEnergyNeutralNotRemoved = 0;
1143     fEnergyNeutralRemoved = 0;
1144
1145     fChargedNotRemoved = 0;
1146     fChargedRemoved = 0;
1147     fNeutralNotRemoved = 0;
1148     fNeutralRemoved = 0;
1149     fGammaRemoved = 0;
1150     fSecondaryNotRemoved = 0;
1151
1152     fTrackMultInAcc = 0;
1153
1154     fTotNeutralEtAfterMinEnergyCut = 0;
1155
1156     fSecondaryNotRemoved = 0;
1157
1158     fTotPx = 0;
1159     fTotPy = 0;
1160     fTotPz = 0;
1161
1162
1163 }
1164
1165 void AliAnalysisEtMonteCarlo::CreateHistograms()
1166 {   // histogram related additions
1167     AliAnalysisEt::CreateHistograms();
1168     if (fEventSummaryTree) {
1169         fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1170         fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1171         fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1172         fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1173         fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1174         fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1175         fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
1176         fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1177         fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1178         fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1179         fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1180         fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1181         fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1182         fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1183         fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1184         fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
1185 //      fEventSummaryTree->Branch("f
1186     }
1187
1188     //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1189     //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1190     //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1191     //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1192
1193     fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1194
1195     fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1196     fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1197     fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1198     fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1199     fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1200     fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1201     fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
1202     fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
1203     fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1204     fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1205     fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1206     fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1207     fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1208     fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1209     fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1210
1211     fHistEtNonRemovedGammas = new  TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1212     fHistEtNonRemovedGammasFromPi0 = new  TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1213
1214     fHistEtRemovedGammas = new  TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1215     fHistEtRemovedNeutrons = new  TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1216     fHistEtRemovedAntiNeutrons = new  TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1217
1218     fHistEtRemovedCharged = new  TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1219     fHistEtRemovedNeutrals = new  TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1220
1221     fHistEtNonRemovedCharged = new  TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1222     fHistEtNonRemovedNeutrals = new  TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1223
1224     fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1225     fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1226     fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1227     fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1228     fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1229     fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1230     fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
1231     fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
1232     fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1233     fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1234     fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1235     fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1236     fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1237     fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1238     fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1239
1240     fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1241
1242     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1243     fHistMultRemovedNeutrons = new  TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1244     fHistMultRemovedAntiNeutrons = new  TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1245     /*
1246       fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1247       fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1248
1249       fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1250       fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
1251
1252
1253     fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1254     fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1255
1256     fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1257     fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1258
1259     fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1260     fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1261     fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1262
1263     fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1264     fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1265     fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1266
1267     fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1268     fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1269     fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1270     fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1271
1272     fHistK0EDepositsVsPtForEOver000MeV = new TH2F("fHistK0EDepositsVsPtForEOver000MeV","K^{0}_{S} deposits over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1273     fHistK0EDepositsVsPtForEOver100MeV = new TH2F("fHistK0EDepositsVsPtForEOver100MeV","K^{0}_{S} deposits over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1274     fHistK0EDepositsVsPtForEOver150MeV = new TH2F("fHistK0EDepositsVsPtForEOver150MeV","K^{0}_{S} deposits over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1275     fHistK0EDepositsVsPtForEOver200MeV = new TH2F("fHistK0EDepositsVsPtForEOver200MeV","K^{0}_{S} deposits over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1276     fHistK0EDepositsVsPtForEOver250MeV = new TH2F("fHistK0EDepositsVsPtForEOver250MeV","K^{0}_{S} deposits over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1277     fHistK0EDepositsVsPtForEOver300MeV = new TH2F("fHistK0EDepositsVsPtForEOver300MeV","K^{0}_{S} deposits over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1278     fHistK0EDepositsVsPtForEOver350MeV = new TH2F("fHistK0EDepositsVsPtForEOver350MeV","K^{0}_{S} deposits over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1279     fHistK0EDepositsVsPtForEOver400MeV = new TH2F("fHistK0EDepositsVsPtForEOver400MeV","K^{0}_{S} deposits over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1280     fHistK0EDepositsVsPtForEOver450MeV = new TH2F("fHistK0EDepositsVsPtForEOver450MeV","K^{0}_{S} deposits over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1281     fHistK0EDepositsVsPtForEOver500MeV = new TH2F("fHistK0EDepositsVsPtForEOver500MeV","K^{0}_{S} deposits over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1282     fHistK0EGammaVsPtForEOver000MeV = new TH2F("fHistK0EGammaVsPtForEOver000MeV","K^{0}_{S} gamma over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1283     fHistK0EGammaVsPtForEOver100MeV = new TH2F("fHistK0EGammaVsPtForEOver100MeV","K^{0}_{S} gamma over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1284     fHistK0EGammaVsPtForEOver150MeV = new TH2F("fHistK0EGammaVsPtForEOver150MeV","K^{0}_{S} gamma over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1285     fHistK0EGammaVsPtForEOver200MeV = new TH2F("fHistK0EGammaVsPtForEOver200MeV","K^{0}_{S} gamma over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1286     fHistK0EGammaVsPtForEOver250MeV = new TH2F("fHistK0EGammaVsPtForEOver250MeV","K^{0}_{S} gamma over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1287     fHistK0EGammaVsPtForEOver300MeV = new TH2F("fHistK0EGammaVsPtForEOver300MeV","K^{0}_{S} gamma over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1288     fHistK0EGammaVsPtForEOver350MeV = new TH2F("fHistK0EGammaVsPtForEOver350MeV","K^{0}_{S} gamma over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1289     fHistK0EGammaVsPtForEOver400MeV = new TH2F("fHistK0EGammaVsPtForEOver400MeV","K^{0}_{S} gamma over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1290     fHistK0EGammaVsPtForEOver450MeV = new TH2F("fHistK0EGammaVsPtForEOver450MeV","K^{0}_{S} gamma over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1291     fHistK0EGammaVsPtForEOver500MeV = new TH2F("fHistK0EGammaVsPtForEOver500MeV","K^{0}_{S} gamma over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1292
1293     fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1294     fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1295     fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1296
1297     fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1298     fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1299     fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1300
1301     if(fCuts->GetHistMakeTree())
1302     {
1303         TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1304         fPrimaryTree = new TTree(treename, treename);
1305
1306         fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1307         fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1308         fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1309
1310         fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1311         fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1312
1313         fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1314         fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1315
1316         fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1317         fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1318         fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1319
1320         fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1321         fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1322         fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1323
1324         fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
1325         fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
1326
1327
1328         fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1329         fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
1330         fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
1331         fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1332
1333         fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1334         fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1335         fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1336
1337         fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
1338
1339
1340         fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
1341         fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
1342
1343         fPrimaryTree->Branch("fClusterMult", &fClusterMult,  "fClusterMult/I");
1344
1345
1346     }
1347
1348     fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
1349     fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
1350     fHistChargedTracksCut = new TH1F("fHistChargedTracksCut", "fHistChargedTracksCut",100, 0, 5);
1351     fHistChargedTracksAccepted = new TH1F("fHistChargedTracksAccepted", "fHistChargedTracksAccepted",100, 0, 5);
1352     fHistGammasCut = new TH1F("fHistGammasTracksCut", "fHistGammasTracksCut",100, 0, 5);
1353     fHistGammasAccepted = new TH1F("fHistGammasTracksAccepted", "fHistGammasTracksAccepted",100, 0, 5);
1354     fHistBadTrackMatches = new TH1F("fHistBadTrackMatches", "fHistBadTrackMatches",100, 0, 5);
1355     fHistMatchedTracksEvspTBkgd = new TH2F("fHistMatchedTracksEvspTBkgd", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1356     fHistMatchedTracksEvspTSignal = new TH2F("fHistMatchedTracksEvspTSignal", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
1357     fHistMatchedTracksEvspTBkgdPeripheral = new TH2F("fHistMatchedTracksEvspTBkgdPeripheral", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1358     fHistMatchedTracksEvspTSignalPeripheral = new TH2F("fHistMatchedTracksEvspTSignalPeripheral", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
1359     fHistChargedTracksCutPeripheral = new TH1F("fHistChargedTracksCutPeripheral", "fHistChargedTracksCut",100, 0, 5);
1360     fHistChargedTracksAcceptedPeripheral = new TH1F("fHistChargedTracksAcceptedPeripheral", "fHistChargedTracksAccepted",100, 0, 5);
1361     fHistGammasCutPeripheral = new TH1F("fHistGammasTracksCutPeripheral", "fHistGammasTracksCut",100, 0, 5);
1362     fHistGammasAcceptedPeripheral = new TH1F("fHistGammasTracksAcceptedPeripheral", "fHistGammasTracksAccepted",100, 0, 5);
1363     fHistBadTrackMatchesdPhidEta = new TH2F("fHistBadTrackMatchesdPhidEta", "fHistBadTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
1364     fHistGoodTrackMatchesdPhidEta = new TH2F("fHistGoodTrackMatchesdPhidEta", "fHistGoodTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
1365 }
1366
1367 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
1368 {   //fill the output list
1369     AliAnalysisEt::FillOutputList(list);
1370
1371     if(fCuts->GetHistMakeTree())
1372     {
1373         list->Add(fPrimaryTree);
1374     }
1375
1376     list->Add(fHistRemovedOrNot);
1377
1378     list->Add(fHistEtNonRemovedProtons);
1379     list->Add(fHistEtNonRemovedAntiProtons);
1380     list->Add(fHistEtNonRemovedPiPlus);
1381     list->Add(fHistEtNonRemovedPiMinus);
1382     list->Add(fHistEtNonRemovedKaonPlus);
1383     list->Add(fHistEtNonRemovedKaonMinus);
1384     list->Add(fHistEtNonRemovedK0s);
1385     list->Add(fHistEtNonRemovedK0L);
1386     list->Add(fHistEtNonRemovedLambdas);
1387     list->Add(fHistEtNonRemovedElectrons);
1388     list->Add(fHistEtNonRemovedPositrons);
1389     list->Add(fHistEtNonRemovedMuPlus);
1390     list->Add(fHistEtNonRemovedMuMinus);
1391     list->Add(fHistEtNonRemovedNeutrons);
1392     list->Add(fHistEtNonRemovedAntiNeutrons);
1393     list->Add(fHistEtNonRemovedGammas);
1394     list->Add(fHistEtNonRemovedGammasFromPi0);
1395
1396     list->Add(fHistEtRemovedGammas);
1397     list->Add(fHistEtRemovedNeutrons);
1398     list->Add(fHistEtRemovedAntiNeutrons);
1399
1400     list->Add(fHistEtRemovedCharged);
1401     list->Add(fHistEtRemovedNeutrals);
1402
1403     list->Add(fHistEtNonRemovedCharged);
1404     list->Add(fHistEtNonRemovedNeutrals);
1405
1406     list->Add(fHistMultNonRemovedProtons);
1407     list->Add(fHistMultNonRemovedAntiProtons);
1408     list->Add(fHistMultNonRemovedPiPlus);
1409     list->Add(fHistMultNonRemovedPiMinus);
1410     list->Add(fHistMultNonRemovedKaonPlus);
1411     list->Add(fHistMultNonRemovedKaonMinus);
1412     list->Add(fHistMultNonRemovedK0s);
1413     list->Add(fHistMultNonRemovedK0L);
1414     list->Add(fHistMultNonRemovedLambdas);
1415     list->Add(fHistMultNonRemovedElectrons);
1416     list->Add(fHistMultNonRemovedPositrons);
1417     list->Add(fHistMultNonRemovedMuPlus);
1418     list->Add(fHistMultNonRemovedMuMinus);
1419     list->Add(fHistMultNonRemovedNeutrons);
1420     list->Add(fHistMultNonRemovedAntiNeutrons);
1421     list->Add(fHistMultNonRemovedGammas);
1422
1423     list->Add(fHistMultRemovedGammas);
1424     list->Add(fHistMultRemovedNeutrons);
1425     list->Add(fHistMultRemovedAntiNeutrons);
1426
1427     list->Add(fHistMultRemovedCharged);
1428     list->Add(fHistMultRemovedNeutrals);
1429
1430     list->Add(fHistMultNonRemovedCharged);
1431     list->Add(fHistMultNonRemovedNeutrals);
1432
1433     list->Add(fHistTrackMultvsNonRemovedCharged);
1434     list->Add(fHistTrackMultvsNonRemovedNeutral);
1435     list->Add(fHistTrackMultvsRemovedGamma);
1436
1437     list->Add(fHistClusterMultvsNonRemovedCharged);
1438     list->Add(fHistClusterMultvsNonRemovedNeutral);
1439     list->Add(fHistClusterMultvsRemovedGamma);
1440
1441     //list->Add(fHistDecayVertexNonRemovedCharged);
1442     //list->Add(fHistDecayVertexNonRemovedNeutral);
1443     //list->Add(fHistDecayVertexRemovedCharged);
1444     //list->Add(fHistDecayVertexRemovedNeutral);
1445
1446     list->Add(fHistDxDzNonRemovedCharged);
1447     list->Add(fHistDxDzRemovedCharged);
1448     list->Add(fHistDxDzNonRemovedNeutral);
1449     list->Add(fHistDxDzRemovedNeutral);
1450
1451     list->Add(fHistK0EDepositsVsPtForEOver000MeV);
1452     list->Add(fHistK0EDepositsVsPtForEOver100MeV);
1453     list->Add(fHistK0EDepositsVsPtForEOver150MeV);
1454     list->Add(fHistK0EDepositsVsPtForEOver200MeV);
1455     list->Add(fHistK0EDepositsVsPtForEOver250MeV);
1456     list->Add(fHistK0EDepositsVsPtForEOver300MeV);
1457     list->Add(fHistK0EDepositsVsPtForEOver350MeV);
1458     list->Add(fHistK0EDepositsVsPtForEOver400MeV);
1459     list->Add(fHistK0EDepositsVsPtForEOver450MeV);
1460     list->Add(fHistK0EDepositsVsPtForEOver500MeV);
1461     list->Add(fHistK0EGammaVsPtForEOver000MeV);
1462     list->Add(fHistK0EGammaVsPtForEOver100MeV);
1463     list->Add(fHistK0EGammaVsPtForEOver150MeV);
1464     list->Add(fHistK0EGammaVsPtForEOver200MeV);
1465     list->Add(fHistK0EGammaVsPtForEOver250MeV);
1466     list->Add(fHistK0EGammaVsPtForEOver300MeV);
1467     list->Add(fHistK0EGammaVsPtForEOver350MeV);
1468     list->Add(fHistK0EGammaVsPtForEOver400MeV);
1469     list->Add(fHistK0EGammaVsPtForEOver450MeV);
1470     list->Add(fHistK0EGammaVsPtForEOver500MeV);
1471
1472     list->Add(fHistPiPlusMult);
1473     list->Add(fHistPiMinusMult);
1474     list->Add(fHistPiZeroMult);
1475     list->Add(fHistPiPlusMultAcc);
1476     list->Add(fHistPiMinusMultAcc);
1477     list->Add(fHistPiZeroMultAcc);
1478
1479     list->Add(fHistGammasFound);
1480     list->Add(fHistGammasGenerated);
1481     list->Add(fHistChargedTracksCut);
1482     list->Add(fHistChargedTracksAccepted);
1483     list->Add(fHistGammasCut);
1484     list->Add(fHistGammasAccepted);
1485     list->Add(fHistBadTrackMatches);
1486     list->Add(fHistMatchedTracksEvspTBkgd);
1487     list->Add(fHistMatchedTracksEvspTSignal);
1488     list->Add(fHistMatchedTracksEvspTBkgdPeripheral);
1489     list->Add(fHistMatchedTracksEvspTSignalPeripheral);
1490     list->Add(fHistChargedTracksCutPeripheral);
1491     list->Add(fHistChargedTracksAcceptedPeripheral);
1492     list->Add(fHistGammasCutPeripheral);
1493     list->Add(fHistGammasAcceptedPeripheral);
1494     list->Add(fHistBadTrackMatchesdPhidEta);
1495     list->Add(fHistGoodTrackMatchesdPhidEta);
1496
1497 }
1498
1499
1500 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1501 {
1502     //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
1503     AliESDtrack *esdTrack = new AliESDtrack(part);
1504     // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1505
1506     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1507
1508     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1509
1510     bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
1511     delete esdTrack;
1512
1513     return status;
1514 }
1515
1516 void AliAnalysisEtMonteCarlo::FillHistograms()
1517 {   // let base class fill its histograms, and us fill the local ones
1518     AliAnalysisEt::FillHistograms();
1519     //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1520
1521     fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1522     fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1523     fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1524     fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
1525     fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1526     fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
1527     fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1528     fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1529     fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1530     fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1531     fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1532     fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1533     fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1534     fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1535     fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1536     fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1537     fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1538
1539     fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
1540     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1541     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1542
1543 //     fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1544 //                                             +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1545 //                                             +fEtRemovedProtons.
1546 //                              fCentClass);
1547 //     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
1548 //
1549 //     fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
1550 //                                             +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
1551 //                                             +fEtNonRemovedProtons,
1552 //                              fCentClass);
1553 //     fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
1554
1555     fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
1556     fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
1557     fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
1558     fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
1559
1560     fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1561     fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1562     fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1563     fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
1564
1565
1566     fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1567     fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1568     fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1569     fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1570     fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
1571     fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
1572     fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1573     fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1574     fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1575     fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1576     fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1577     fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1578     fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1579     fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1580     fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1581     fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1582
1583     fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1584     fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1585     fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1586
1587     fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1588                                             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1589                                             +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1590                                             +fMultNonRemovedProtons);
1591
1592     fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1593                                             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
1594
1595     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1596                                        fMultRemovedGammas);
1597
1598     fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
1599             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1600             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1601             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1602
1603     fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
1604             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
1605
1606     fHistClusterMultvsRemovedGamma->Fill(fNClusters,
1607                                          fMultRemovedGammas);
1608
1609 }
1610
1611
1612
1613
1614 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
1615 { // print family tree
1616     TParticle *part = stack->Particle(partIdx);
1617 //     if(part->GetPdgCode() == fgK0SCode)
1618     {
1619         std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
1620         std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
1621         std::cout << "Energy: " << part->Energy() << std::endl;
1622         std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
1623     }
1624     return PrintMothers(partIdx, stack, 1);
1625 }
1626
1627 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
1628 { // print mothers
1629     char *tabs = new char[gen+1];
1630     for(Int_t i = 0; i < gen; ++i)
1631     {
1632         //std::cout << i << std::endl;
1633         tabs[i] = '\t';
1634     }
1635     tabs[gen] = '\0';
1636     Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1637     if(mothIdx < 0)
1638     {
1639       delete [] tabs;
1640       return 0;
1641     }
1642     TParticle *mother = stack->Particle(mothIdx);
1643 //     if(mother->GetPdgCode() == fgK0SCode)
1644     {
1645         //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
1646         std::cout << tabs << "Index: " << mothIdx << std::endl;
1647         std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
1648         std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1649         std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
1650         if(mother->GetFirstMother() >= 0)
1651         {
1652           std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
1653           if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
1654           std::cout << std::endl;
1655         }
1656         std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1657     }
1658     if(mother->GetPdgCode() == fgK0SCode)
1659     {
1660 //      std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
1661     }
1662 //  std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
1663 //  std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1664 //  std::cout << "Energy: " << mother->Energy() << std::endl;
1665 //  std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1666
1667     delete [] tabs;
1668     return PrintMothers(mothIdx, stack, gen+1) + 1;
1669 }
1670
1671 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
1672 { // get primary mother
1673     if(partIdx >= 0)
1674     {
1675         //return stack->GetPrimary(partIdx);
1676
1677         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1678         if(mothIdx < 0) return -1;
1679         TParticle *mother = stack->Particle(mothIdx);
1680         if(mother)
1681         {
1682             if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1683             else return GetPrimMother(mothIdx, stack);
1684         }
1685         else
1686         {
1687             return -1;
1688         }
1689     }
1690     return -1;
1691 }
1692
1693 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
1694 { // get K0 in family
1695     if(partIdx >= 0)
1696     {
1697         if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
1698         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1699         if(mothIdx < 0) return -1;
1700         TParticle *mother = stack->Particle(mothIdx);
1701         if(mother)
1702         {
1703           if(mother->GetPdgCode() == fgK0SCode)
1704             {
1705                 return mothIdx;
1706             }
1707             return GetK0InFamily(mothIdx, stack);
1708         }
1709         else
1710         {
1711             return -1;
1712         }
1713     }
1714     return -1;
1715 }
1716