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