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