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