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