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
  39         ,fHistDecayVertexNonRemovedCharged(0)
  40         ,fHistDecayVertexRemovedCharged(0)
  41         ,fHistDecayVertexNonRemovedNeutral(0)
  42         ,fHistDecayVertexRemovedNeutral(0)
  43
  44         ,fHistRemovedOrNot(0)
  45         ,fHistEtNonRemovedProtons(0)
  46         ,fHistEtNonRemovedAntiProtons(0)
  47         ,fHistEtNonRemovedPiPlus(0)
  48         ,fHistEtNonRemovedPiMinus(0)
  49         ,fHistEtNonRemovedKaonPlus(0)
  50         ,fHistEtNonRemovedKaonMinus(0)
  51         ,fHistEtNonRemovedK0s(0)
  52         ,fHistEtNonRemovedLambdas(0)
  53         ,fHistEtNonRemovedElectrons(0)
  54         ,fHistEtNonRemovedPositrons(0)
  55         ,fHistEtNonRemovedMuPlus(0)
  56         ,fHistEtNonRemovedMuMinus(0)
  57         ,fHistEtNonRemovedNeutrons(0)
  58         ,fHistEtNonRemovedAntiNeutrons(0)
  59         ,fHistEtNonRemovedGammas(0)
  60         ,fHistEtNonRemovedGammasFromPi0(0)
  61         ,fHistEtRemovedGammas(0)
  62         ,fHistEtRemovedNeutrons(0)
  63         ,fHistEtRemovedAntiNeutrons(0)
  64         ,fHistMultNonRemovedProtons(0)
  65         ,fHistMultNonRemovedAntiProtons(0)
  66         ,fHistMultNonRemovedPiPlus(0)
  67         ,fHistMultNonRemovedPiMinus(0)
  68         ,fHistMultNonRemovedKaonPlus(0)
  69         ,fHistMultNonRemovedKaonMinus(0)
  70         ,fHistMultNonRemovedK0s(0)
  71         ,fHistMultNonRemovedLambdas(0)
  72         ,fHistMultNonRemovedElectrons(0)
  73         ,fHistMultNonRemovedPositrons(0)
  74         ,fHistMultNonRemovedMuPlus(0)
  75         ,fHistMultNonRemovedMuMinus(0)
  76         ,fHistMultNonRemovedNeutrons(0)
  77         ,fHistMultNonRemovedAntiNeutrons(0)
  78         ,fHistMultNonRemovedGammas(0)
  79         ,fHistMultRemovedGammas(0)
  80         ,fHistMultRemovedNeutrons(0)
  81         ,fHistMultRemovedAntiNeutrons(0)
  82         ,fHistTrackMultvsNonRemovedCharged(0)
  83         ,fHistTrackMultvsNonRemovedNeutral(0)
  84         ,fHistTrackMultvsRemovedGamma(0)
  85         ,fHistClusterMultvsNonRemovedCharged(0)
  86         ,fHistClusterMultvsNonRemovedNeutral(0)
  87         ,fHistClusterMultvsRemovedGamma(0)
  88         ,fHistMultvsNonRemovedChargedE(0)
  89         ,fHistMultvsNonRemovedNeutralE(0)
  90         ,fHistMultvsRemovedGammaE(0)
  91         ,fEtNonRemovedProtons(0)
  92         ,fEtNonRemovedAntiProtons(0)
  93         ,fEtNonRemovedPiPlus(0)
  94         ,fEtNonRemovedPiMinus(0)
  95         ,fEtNonRemovedKaonPlus(0)
  96         ,fEtNonRemovedKaonMinus(0)
  97         ,fEtNonRemovedK0s(0)
  98         ,fEtNonRemovedLambdas(0)
  99         ,fEtNonRemovedElectrons(0)
 100         ,fEtNonRemovedPositrons(0)
 101         ,fEtNonRemovedMuMinus(0)
 102         ,fEtNonRemovedMuPlus(0)
 103         ,fEtNonRemovedGammas(0)
 104         ,fEtNonRemovedGammasFromPi0(0)
 105         ,fEtNonRemovedNeutrons(0)
 106         ,fEtNonRemovedAntiNeutrons(0)
 107         ,fEtRemovedGammas(0)
 108         ,fEtRemovedNeutrons(0)
 109         ,fEtRemovedAntiNeutrons(0)
 110         ,fMultNonRemovedProtons(0)
 111         ,fMultNonRemovedAntiProtons(0)
 112         ,fMultNonRemovedPiPlus(0)
 113         ,fMultNonRemovedPiMinus(0)
 114         ,fMultNonRemovedKaonPlus(0)
 115         ,fMultNonRemovedKaonMinus(0)
 116         ,fMultNonRemovedK0s(0)
 117         ,fMultNonRemovedLambdas(0)
 118         ,fMultNonRemovedElectrons(0)
 119         ,fMultNonRemovedPositrons(0)
 120         ,fMultNonRemovedMuMinus(0)
 121         ,fMultNonRemovedMuPlus(0)
 122         ,fMultNonRemovedGammas(0)
 123         ,fMultNonRemovedNeutrons(0)
 124         ,fMultNonRemovedAntiNeutrons(0)
 125         ,fMultRemovedGammas(0)
 126         ,fMultRemovedNeutrons(0)
 127         ,fMultRemovedAntiNeutrons(0)
 128         ,fTrackMultInAcc(0)
 129         ,fHistDxDzNonRemovedCharged(0)
 130         ,fHistDxDzRemovedCharged(0)
 131         ,fHistDxDzNonRemovedNeutral(0)
 132         ,fHistDxDzRemovedNeutral(0)
 133         ,fHistPiPlusMult(0)
 134         ,fHistPiMinusMult(0)
 135         ,fHistPiZeroMult(0)
 136         ,fHistPiPlusMultAcc(0)
 137         ,fHistPiMinusMultAcc(0)
 138         ,fHistPiZeroMultAcc(0)
 139         ,fPiPlusMult(0)
 140         ,fPiMinusMult(0)
 141         ,fPiZeroMult(0)
 142         ,fPiPlusMultAcc(0)
 143         ,fPiMinusMultAcc(0)
 144         ,fPiZeroMultAcc(0)
 145         ,fNeutralRemoved(0)
 146         ,fChargedRemoved(0)
 147         ,fChargedNotRemoved(0)
 148         ,fNeutralNotRemoved(0)
 149         ,fEnergyNeutralRemoved(0)
 150         ,fEnergyChargedRemoved(0)
 151         ,fEnergyChargedNotRemoved(0)
 152         ,fEnergyNeutralNotRemoved(0)
 153
 154 {
 155     fTrackDistanceCut = 7.0;
 156 }
 157
 158 // dtor
 159 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
 160 {//Destructor
 161         delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
 162         delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
 163         delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
 164         delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
 165
 166         delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
 167
 168         delete fHistEtNonRemovedProtons; // enter comment here
 169         delete fHistEtNonRemovedAntiProtons; // enter comment here
 170         delete fHistEtNonRemovedPiPlus; // enter comment here
 171         delete fHistEtNonRemovedPiMinus; // enter comment here
 172         delete fHistEtNonRemovedKaonPlus; // enter comment here
 173         delete fHistEtNonRemovedKaonMinus; // enter comment here
 174         delete fHistEtNonRemovedK0s; // enter comment here
 175         delete fHistEtNonRemovedLambdas; // enter comment here
 176         delete fHistEtNonRemovedElectrons; // enter comment here
 177         delete fHistEtNonRemovedPositrons; // enter comment here
 178         delete fHistEtNonRemovedMuPlus; // enter comment here
 179         delete fHistEtNonRemovedMuMinus; // enter comment here
 180         delete fHistEtNonRemovedNeutrons; // enter comment here
 181         delete fHistEtNonRemovedAntiNeutrons; // enter comment here
 182         delete fHistEtNonRemovedGammas; // enter comment here
 183         delete fHistEtNonRemovedGammasFromPi0; // enter comment here
 184
 185         delete fHistEtRemovedGammas; // enter comment here
 186         delete fHistEtRemovedNeutrons; // enter comment here
 187         delete fHistEtRemovedAntiNeutrons; // enter comment here
 188
 189
 190         delete fHistMultNonRemovedProtons; // enter comment here
 191         delete fHistMultNonRemovedAntiProtons; // enter comment here
 192         delete fHistMultNonRemovedPiPlus; // enter comment here
 193         delete fHistMultNonRemovedPiMinus; // enter comment here
 194         delete fHistMultNonRemovedKaonPlus; // enter comment here
 195         delete fHistMultNonRemovedKaonMinus; // enter comment here
 196         delete fHistMultNonRemovedK0s; // enter comment here
 197         delete fHistMultNonRemovedLambdas; // enter comment here
 198         delete fHistMultNonRemovedElectrons; // enter comment here
 199         delete fHistMultNonRemovedPositrons; // enter comment here
 200         delete fHistMultNonRemovedMuPlus; // enter comment here
 201         delete fHistMultNonRemovedMuMinus; // enter comment here
 202         delete fHistMultNonRemovedNeutrons; // enter comment here
 203         delete fHistMultNonRemovedAntiNeutrons; // enter comment here
 204         delete fHistMultNonRemovedGammas; // enter comment here
 205
 206         delete fHistMultRemovedGammas; // enter comment here
 207         delete fHistMultRemovedNeutrons; // enter comment here
 208         delete fHistMultRemovedAntiNeutrons; // enter comment here
 209
 210         delete fHistTrackMultvsNonRemovedCharged; // enter comment here
 211         delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
 212         delete fHistTrackMultvsRemovedGamma; // enter comment here
 213
 214         delete fHistClusterMultvsNonRemovedCharged; // enter comment here
 215         delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
 216         delete fHistClusterMultvsRemovedGamma; // enter comment here
 217
 218         delete fHistMultvsNonRemovedChargedE; // enter comment here
 219         delete fHistMultvsNonRemovedNeutralE; // enter comment here
 220         delete fHistMultvsRemovedGammaE; // enter comment here
 221         delete fHistDxDzNonRemovedCharged; // enter comment here
 222         delete fHistDxDzRemovedCharged; // enter comment here
 223         delete fHistDxDzNonRemovedNeutral; // enter comment here
 224         delete fHistDxDzRemovedNeutral; // enter comment here
 225
 226         delete fHistPiPlusMult; // enter comment here
 227         delete fHistPiMinusMult; // enter comment here
 228         delete fHistPiZeroMult; // enter comment here
 229
 230         delete fHistPiPlusMultAcc; // enter comment here
 231         delete fHistPiMinusMultAcc; // enter comment here
 232         delete fHistPiZeroMultAcc; // enter comment here
 233 }
 234
 235 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 236 { // analyse MC event
 237     ResetEventValues();
 238
 239     fPiPlusMult = 0;
 240     fPiMinusMult = 0;
 241     fPiZeroMult = 0;
 242     if (fCentrality)
 243     {
 244         fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
 245
 246     }
 247     fSparseTracks[0] = 0;
 248     fSparseTracks[1] = 0;
 249     fSparseTracks[2] = 0;
 250     fSparseTracks[3] = 0;
 251     fSparseTracks[4] = 0;
 252     fSparseTracks[5] = 0;
 253     fSparseTracks[6] = fCentClass;
 254
 255
 256     // Get us an mc event
 257     if (!ev) {
 258         AliFatal("ERROR: Event does not exist");
 259         return 0;
 260     }
 261     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
 262     if (!event) {
 263         AliFatal("ERROR: MC Event does not exist");
 264         return 0;
 265     }
 266
 267     Double_t protonMass =fgProtonMass;
 268
 269     // Hijing header
 270     AliGenEventHeader* genHeader = event->GenEventHeader();
 271     if (!genHeader) {
 272         Printf("ERROR: Event generation header does not exist");
 273         return 0;
 274     }
 275     AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 276     if (hijingGenHeader) {
 277         fImpactParameter = hijingGenHeader->ImpactParameter();
 278         fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
 279         fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
 280         /*
 281           printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
 282           hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
 283           printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
 284           hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
 285           printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
 286           hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
 287           printf("NN %d NNw %d NwN %d, NwNw %d\n",
 288           hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
 289         */
 290     }
 291
 292     /* // placeholder if we want to get some Pythia info later
 293        AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 294        if (pythiaGenHeader) { // not Hijing; try with Pythia
 295        printf("Pythia: ProcessType %d  GetPtHard %g \n",
 296        pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
 297        }
 298     */
 299
 300     // Let's play with the stack!
 301     AliStack *stack = event->Stack();
 302
 303     Int_t nPrim = stack->GetNtrack();
 304
 305     Int_t partCount = 0;
 306     for (Int_t iPart = 0; iPart < nPrim; iPart++)
 307     {
 308
 309         TParticle *part = stack->Particle(iPart);
 310 //         TParticle *partMom = 0;
 311 //         TParticle *partMomLast = 0;
 312
 313         if (!part)
 314         {
 315             Printf("ERROR: Could not get particle %d", iPart);
 316             continue;
 317         }
 318
 319 //         Int_t iPartMom = part->GetMother(0);
 320 //         Int_t iPartLastMom = part->GetMother(1);
 321
 322         TParticlePDG *pdg = part->GetPDG(0);
 323         //TParticlePDG *pdgMom = 0;
 324 //         TParticlePDG *pdgMomLast = 0;
 325
 326         if (!pdg)
 327         {
 328             //Printf("ERROR: Could not get particle PDG %d", iPart);
 329             continue;
 330         }
 331
 332 //         if (iPartMom>0)
 333 //         {
 334 //             partMom = stack->Particle(iPartMom);
 335 //             //pdgMom = partMom->GetPDG(0);
 336 //         }
 337
 338 //         if (iPartLastMom>0)
 339 //         {
 340 //             partMomLast = stack->Particle(iPartLastMom);
 341 //             pdgMomLast = partMomLast->GetPDG(0);
 342 //         }
 343
 344
 345         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
 346
 347         // Check if it is a primary particle
 348         //if (!stack->IsPhysicalPrimary(iPart)) continue;
 349
 350         //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
 351
 352         // Check for reasonable (for now neutral and singly charged) charge on the particle
 353         //TODO:Maybe not only singly charged?
 354         if ((stack->IsPhysicalPrimary(iPart))&&( TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 || TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3))
 355         {
 356
 357             fMultiplicity++;
 358
 359             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 360             {
 361                 // calculate E_T
 362                 if (
 363                     TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
 364                     TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
 365                     TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
 366                     TMath::Abs(pdg->PdgCode()) == fgXiCode ||
 367                     TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
 368                     TMath::Abs(pdg->PdgCode()) == fgOmegaCode
 369                 )
 370                 {
 371                     if (pdg->PdgCode() > 0) {
 372                         particleMassPart = - protonMass;
 373                     }
 374                     if (pdg->PdgCode() < 0) {
 375                         particleMassPart = protonMass;
 376                     }
 377                 }
 378                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 379
 380                 fSparseTracks[0] = pdg->PdgCode();
 381                 fSparseTracks[1] = pdg->Charge()*3;
 382                 fSparseTracks[2] = pdg->Mass();
 383                 fSparseTracks[3] = et;
 384                 fSparseTracks[4] = part->Pt();
 385                 fSparseTracks[5] = part->Eta();
 386                 fSparseHistTracks->Fill(fSparseTracks);
 387
 388                 // Fill up total E_T counters for each particle species
 389                 if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 390                 {
 391                     fProtonEt += et;
 392                 }
 393                 if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 394                 {
 395                     fPionEt += et;
 396                     if (pdg->PdgCode() == fgPiPlusCode)
 397                     {
 398                         fPiPlusMult++;
 399                     }
 400                     else
 401                     {
 402                         fPiMinusMult++;
 403                     }
 404                 }
 405                 if (pdg->PdgCode() == fgGammaCode)
 406                 {
 407                     fPiZeroMult++;
 408                 }
 409                 if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 410                 {
 411                     fChargedKaonEt += et;
 412                 }
 413                 if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 414                 {
 415                     fMuonEt += et;
 416                 }
 417                 if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 418                 {
 419                     fElectronEt += et;
 420                 }
 421
 422                 // some neutrals also
 423                 if (pdg->PdgCode() == fgNeutronCode)
 424                 {
 425                     fNeutronEt += et;
 426                 }
 427                 if (pdg->PdgCode() == fgAntiNeutronCode)
 428                 {
 429                     fAntiNeutronEt += et;
 430                 }
 431                 if (pdg->PdgCode() == fgGammaCode)
 432                 {
 433                     fGammaEt += et;
 434                 }
 435
 436                 // Neutral particles
 437                 if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
 438                 {
 439
 440                     if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
 441
 442                     // inside EMCal acceptance
 443                     if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 444                     {
 445                         fNeutralMultiplicity++;
 446                         //std::cout << pdg->PdgCode() << ", " << iPart << ", " <<  part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
 447
 448                         //if(part->GetDaughter(0) > 0) std::cout << stack->Particle(part->GetDaughter(0))->GetPdgCode() << " " << stack->Particle(part->GetDaughter(1))->GetPdgCode() << std::endl;
 449                         if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
 450                         if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
 451                         if (part->Energy() > 0.05) partCount++;
 452                     }
 453                 }
 454                 //Charged particles
 455                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
 456                 {
 457                     fChargedMultiplicity++;
 458                     fTotChargedEt += et;
 459
 460                     // inside EMCal acceptance
 461                     if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 462                     {
 463                         fTotChargedEtAcc += et;
 464                         fTotEtAcc += et;
 465
 466                         if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 467                         {
 468                             fProtonEtAcc += et;
 469                         }
 470                         if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 471                         {
 472                             fPionEtAcc += et;
 473                         }
 474                         if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 475                         {
 476                             fChargedKaonEtAcc += et;
 477                         }
 478                         if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 479                         {
 480                             fMuonEtAcc += et;
 481                         }
 482
 483                         if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 484                         {
 485                           fElectronEtAcc += et;
 486                         }
 487                     } // inside EMCal acceptance
 488
 489                     //    if (TrackHitsCalorimeter(part, event->GetMagneticField()))
 490                     if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 491                     {
 492                         if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 493                         else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 494                     }
 495                 }
 496             }
 497
 498         }
 499     }
 500
 501     fTotEt = fTotChargedEt + fTotNeutralEt;
 502     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 503
 504
 505     //std::cout << "Event done! # of particles: " << partCount << std::endl;
 506     return 0;
 507 }
 508 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 509 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 510 { // analyse MC and real event info
 511     //if(!mcEvent || !realEvent){
 512     if (!ev || !ev2) {
 513         AliFatal("ERROR: Event does not exist");
 514         return 0;
 515     }
 516
 517     fSparseClusters[0] = 0;
 518     fSparseClusters[1] = 0;
 519     fSparseClusters[2] = 0;
 520     fSparseClusters[3] = 0;
 521     fSparseClusters[4] = 0;
 522     fSparseClusters[5] = 0;
 523     fSparseClusters[6] = 0;
 524     fSparseClusters[7] = 0;
 525     fSparseClusters[8] = 0;
 526     fSparseClusters[9] = fCentClass;
 527     fSparseClusters[10] = 0;
 528
 529     //AnalyseEvent(mcEvent);
 530     AnalyseEvent(ev);
 531     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 532     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
 533     if (!mcEvent || !realEvent) {
 534         AliFatal("ERROR: mcEvent or realEvent does not exist");
 535         return 0;
 536     }
 537
 538     AliStack *stack = mcEvent->Stack();
 539
 540     // get all detector clusters
 541     TRefArray* caloClusters = new TRefArray();
 542     if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
 543     else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
 544     else {
 545         AliFatal("Detector ID has not been specified");
 546         return -1;
 547     }
 548
 549     // Track loop to fill a pT spectrum
 550     for (Int_t iTracks = 0; iTracks < realEvent->GetNumberOfTracks(); iTracks++)
 551     {
 552         AliESDtrack* track = realEvent->GetTrack(iTracks);
 553         if (!track)
 554         {
 555             Printf("ERROR: Could not receive track %d", iTracks);
 556             continue;
 557         }
 558
 559         if ( track->Phi() < fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180. && track->Phi() > fCuts->GetGeometryPhosPhiAccMinCut()  * TMath::Pi()/180. && TMath::Abs(track->Eta()) < fCuts->GetGeometryPhosEtaAccCut())
 560         {
 561             fTrackMultInAcc++;
 562         }
 563     }
 564
 565     Int_t nCluster = caloClusters->GetEntries();
 566
 567     // loop the clusters
 568     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 569     {
 570         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 571         //Float_t caloE = caloCluster->E();
 572
 573         UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
 574         TParticle *part  = stack->Particle(iPart);
 575         TParticle *partMom = 0;
 576
 577         if (!part)
 578         {
 579             Printf("No MC particle %d", iCluster);
 580             continue;
 581         }
 582
 583         // compare MC and Rec energies for all particles
 584         //fHistAllERecEMC->Fill(part->Energy(),caloE);
 585
 586         TParticlePDG *pdg = part->GetPDG(0);
 587         TParticlePDG *pdgMom = 0;
 588
 589         if (!pdg)
 590         {
 591             Printf("ERROR: Could not get particle PDG %d", iPart);
 592             continue;
 593         }
 594
 595         Int_t iPartMom = part->GetMother(0);
 596
 597         if (iPartMom>0)
 598         {
 599             partMom = stack->Particle(iPartMom);
 600             pdgMom = partMom->GetPDG(0);
 601         }
 602
 603         // Check if it is a primary particle
 604         //if (!stack->IsPhysicalPrimary(iPart)) continue;
 605         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
 606         {
 607             if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
 608                 continue;
 609         } // end of primary particle check
 610
 611         if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
 612
 613         Double_t clEt = CalculateTransverseEnergy(caloCluster);
 614         if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
 615         Float_t pos[3];
 616
 617         caloCluster->GetPosition(pos);
 618         TVector3 cp(pos);
 619
 620         fSparseClusters[0] = pdg->PdgCode();
 621         fSparseClusters[1] = pdg->Charge()/3;
 622         fSparseClusters[2] = pdg->Mass();
 623         fSparseClusters[3] = clEt;
 624         fSparseClusters[4] = caloCluster->E();
 625         fSparseClusters[5] = cp.Eta();
 626         fSparseClusters[6] = part->Energy() * TMath::Sin(part->Theta());
 627         fSparseClusters[7] = part->Pt();
 628         fSparseClusters[8] = part->Eta();
 629         fSparseClusters[9] = fCentClass;
 630         fSparseClusters[10] = 0;
 631         fSparseHistClusters->Fill(fSparseClusters);
 632
 633         //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
 634
 635         if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
 636         {
 637
 638             if (pdg->Charge() != 0)
 639             {
 640                 //std::cout << "Removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
 641                 fChargedRemoved++;
 642                 fEnergyChargedRemoved += caloCluster->E();
 643                 fHistRemovedOrNot->Fill(0.0, fCentClass);
 644                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
 645                 //fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
 646                 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 647
 648
 649             }
 650             else
 651             {
 652                 //std::cout << "Removing neutral: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
 653                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
 654                 fNeutralRemoved++;
 655                 fEnergyNeutralRemoved += caloCluster->E();
 656                 fHistRemovedOrNot->Fill(1.0, fCentClass);
 657                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
 658                 //fHistDecayVertexRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
 659                 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 660
 661                 if (pdg->PdgCode() == fgGammaCode)
 662                 {
 663                     fEtRemovedGammas += clEt;
 664                     fMultRemovedGammas++;
 665                 }
 666                 else if (pdg->PdgCode() == fgNeutronCode)
 667                 {
 668                     fEtRemovedNeutrons += clEt;
 669                     fMultRemovedNeutrons++;
 670                 }
 671                 else if (pdg->PdgCode() == fgAntiNeutronCode)
 672                 {
 673                     fEtRemovedAntiNeutrons += clEt;
 674                     fMultRemovedAntiNeutrons++;
 675                 }
 676                 //else std::cout << "Hmm, what is this (neutral: " << pdg->PdgCode() << std::endl;
 677             }
 678         }
 679         else
 680         {
 681
 682             if (pdg->Charge() != 0)
 683             {
 684                 //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
 685                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
 686                 fChargedNotRemoved++;
 687                 fEnergyChargedNotRemoved += caloCluster->E();
 688                 fHistRemovedOrNot->Fill(2.0, fCentClass);
 689                 //std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
 690                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
 691                 //fHistDecayVertexNonRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
 692                 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 693                 if (pdg->PdgCode() == fgProtonCode)
 694                 {
 695                     //std::cout << clEt << std::endl;
 696                     fEtNonRemovedProtons += clEt;
 697                     fMultNonRemovedProtons++;
 698                 }
 699                 else if (pdg->PdgCode() == fgAntiProtonCode)
 700                 {
 701                     //std::cout << clEt << std::endl;
 702                     fEtNonRemovedAntiProtons += clEt;
 703                     fMultNonRemovedAntiProtons++;
 704                 }
 705                 else if (pdg->PdgCode() == fgPiPlusCode)
 706                 {
 707                     //std::cout << "PI+" <<  clEt << std::endl;
 708                     fEtNonRemovedPiPlus += clEt;
 709                     fMultNonRemovedPiPlus++;
 710                 }
 711                 else if (pdg->PdgCode() == fgPiMinusCode)
 712                 {
 713                     // std::cout << "PI-"  << clEt << std::endl;
 714                     fEtNonRemovedPiMinus += clEt;
 715                     fMultNonRemovedPiMinus++;
 716                 }
 717                 else if (pdg->PdgCode() == fgKPlusCode)
 718                 {
 719                     //std::cout << clEt << std::endl;
 720                     fEtNonRemovedKaonPlus += clEt;
 721                     fMultNonRemovedKaonPlus++;
 722                 }
 723                 else if (pdg->PdgCode() == fgKMinusCode)
 724                 {
 725                     //std::cout << clEt << std::endl;
 726                     fEtNonRemovedKaonMinus += clEt;
 727                     fMultNonRemovedKaonMinus++;
 728                 }
 729                 else if (pdg->PdgCode() == fgEPlusCode)
 730                 {
 731                     //std::cout << clEt << std::endl;
 732                     if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
 733                     {
 734                         fEtNonRemovedPositrons += clEt;
 735                         fMultNonRemovedPositrons++;
 736                     }
 737                 }
 738                 else if (pdg->PdgCode() == fgEMinusCode)
 739                 {
 740                     //std::cout << clEt << std::endl;
 741                     if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
 742                     {
 743                         fEtNonRemovedElectrons += clEt;
 744                         fMultNonRemovedElectrons++;
 745                     }
 746                 }
 747                 else if (pdg->PdgCode() == fgMuPlusCode)
 748                 {
 749                     std::cout << clEt << std::endl;
 750                     fEtNonRemovedMuPlus += clEt;
 751                     fMultNonRemovedMuPlus++;
 752                 }
 753                 else if (pdg->PdgCode() == fgMuMinusCode)
 754                 {
 755                     std::cout << clEt << std::endl;
 756                     fEtNonRemovedMuMinus += clEt;
 757                     fMultNonRemovedMuMinus++;
 758                 }
 759
 760             }
 761             else
 762             {
 763                 //std::cout << "Accepted: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
 764                 //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
 765                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
 766                 fNeutralNotRemoved++;
 767                 fEnergyNeutralNotRemoved += caloCluster->E();
 768                 fHistRemovedOrNot->Fill(3.0, fCentClass);
 769                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
 770                 //fHistDecayVertexNonRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
 771                 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 772                 if (pdg->PdgCode() == fgGammaCode)
 773                 {
 774                     fEtNonRemovedGammas += clEt;
 775                     fMultNonRemovedGammas++;
 776                     if (pdgMom)
 777                     {
 778                         if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
 779                         {
 780 //                      std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() <<  ", E: " << part->Energy() << std::endl;
 781                             fEtNonRemovedGammasFromPi0 += clEt;
 782                         }
 783                     }
 784                 }
 785                 else if (pdg->PdgCode() == fgNeutronCode)
 786                 {
 787                     fEtNonRemovedNeutrons += clEt;
 788                     fMultNonRemovedNeutrons++;
 789                 }
 790                 else if (pdg->PdgCode() == fgAntiNeutronCode)
 791                 {
 792                     fEtNonRemovedAntiNeutrons += clEt;
 793                     fMultNonRemovedAntiNeutrons++;
 794                 }
 795                 else if (pdg->PdgCode() == fgK0LCode || pdg->PdgCode() == fgK0SCode)
 796                 {
 797                     fEtNonRemovedK0s += clEt;
 798                     fMultNonRemovedK0s++;
 799
 800                 }
 801                 else if (TMath::Abs(pdg->PdgCode()) == fgLambdaCode)
 802                 {
 803                     fEtNonRemovedLambdas += clEt;
 804                     fMultNonRemovedLambdas++;
 805                 }
 806                 else std::cout << "Hmm, what is this (neutral not removed): " << pdg->PdgCode() << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
 807             }
 808         }
 809     } // end of loop over clusters
 810
 811     delete caloClusters;
 812     //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
 813     //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
 814     //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
 815     //std::cout << "Number of non-removed charged: " << fChargedNotRemoved << std::endl;
 816     //std::cout << "Number of non-removed neutral: " << fNeutralNotRemoved << std::endl;
 817
 818 //  std::cout << "Energy of removed neutrals: " << fEnergyNeutralRemoved << std::endl;
 819 //  std::cout << "Energy of removed charged: " << fEnergyChargedRemoved << std::endl;
 820 //  std::cout << "Energy of non-removed charged: " << fEnergyChargedNotRemoved << std::endl;
 821 //  std::cout << "Energy of non-removed neutral: " << fEnergyNeutralNotRemoved << std::endl;
 822
 823     FillHistograms();
 824     return 0;
 825 }
 826
 827 void AliAnalysisEtMonteCarlo::Init()
 828 { // init
 829     AliAnalysisEt::Init();
 830 }
 831
 832 void AliAnalysisEtMonteCarlo::ResetEventValues()
 833 { // reset event values
 834     AliAnalysisEt::ResetEventValues();
 835
 836     // collision geometry defaults for p+p:
 837     fImpactParameter = 0;
 838     fNcoll = 1;
 839     fNpart = 2;
 840
 841     fEtNonRemovedProtons = 0;
 842     fEtNonRemovedAntiProtons = 0;
 843     fEtNonRemovedPiPlus = 0;
 844     fEtNonRemovedPiMinus = 0;
 845     fEtNonRemovedKaonPlus = 0;
 846     fEtNonRemovedKaonMinus = 0;
 847     fEtNonRemovedK0s = 0;
 848     fEtNonRemovedLambdas = 0;
 849     fEtNonRemovedElectrons = 0;
 850     fEtNonRemovedPositrons = 0;
 851     fEtNonRemovedMuPlus = 0;
 852     fEtNonRemovedMuMinus = 0;
 853     fEtNonRemovedNeutrons = 0;
 854     fEtNonRemovedAntiNeutrons = 0;
 855     fEtNonRemovedGammas = 0;
 856     fEtNonRemovedGammasFromPi0 = 0;
 857
 858     fEtRemovedGammas = 0;
 859     fEtRemovedNeutrons = 0;
 860     fEtRemovedAntiNeutrons = 0;
 861
 862     fMultNonRemovedProtons = 0;
 863     fMultNonRemovedAntiProtons = 0;
 864     fMultNonRemovedPiPlus = 0;
 865     fMultNonRemovedPiMinus = 0;
 866     fMultNonRemovedKaonPlus = 0;
 867     fMultNonRemovedKaonMinus = 0;
 868     fMultNonRemovedK0s = 0;
 869     fMultNonRemovedLambdas = 0;
 870     fMultNonRemovedElectrons = 0;
 871     fMultNonRemovedPositrons = 0;
 872     fMultNonRemovedMuPlus = 0;
 873     fMultNonRemovedMuMinus = 0;
 874     fMultNonRemovedNeutrons = 0;
 875     fMultNonRemovedAntiNeutrons = 0;
 876     fMultNonRemovedGammas = 0;
 877
 878     fMultRemovedGammas = 0;
 879     fMultRemovedNeutrons = 0;
 880     fMultRemovedAntiNeutrons = 0;
 881
 882     fTrackMultInAcc = 0;
 883
 884 }
 885
 886 void AliAnalysisEtMonteCarlo::CreateHistograms()
 887 { // histogram related additions
 888     AliAnalysisEt::CreateHistograms();
 889     if (fTree) {
 890         fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
 891         fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
 892         fTree->Branch("fNpart",&fNpart,"fNpart/I");
 893     }
 894
 895     //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
 896     //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
 897     //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
 898     //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
 899
 900     fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
 901
 902     fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
 903     fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
 904     fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
 905     fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
 906     fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
 907     fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
 908     fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
 909     fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
 910     fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
 911     fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
 912     fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
 913     fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
 914     fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
 915     fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
 916
 917     fHistEtNonRemovedGammas = new  TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
 918     fHistEtNonRemovedGammasFromPi0 = new  TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
 919
 920     fHistEtRemovedGammas = new  TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
 921     fHistEtRemovedNeutrons = new  TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
 922     fHistEtRemovedAntiNeutrons = new  TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
 923
 924     fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 925     fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 926     fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 927     fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 928     fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 929     fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 930     fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
 931     fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
 932     fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 933     fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 934     fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 935     fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
 936     fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 937     fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 938
 939     fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
 940
 941     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
 942     fHistMultRemovedNeutrons = new  TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 943     fHistMultRemovedAntiNeutrons = new  TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 944
 945     fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 946     fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 947     fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 948
 949     fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 950     fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 951     fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
 952
 953     fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
 954     fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
 955     fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
 956     fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
 957
 958     fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
 959     fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
 960     fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
 961
 962     fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
 963     fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
 964     fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
 965
 966 }
 967
 968 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
 969 {//fill the output list
 970     AliAnalysisEt::FillOutputList(list);
 971
 972     list->Add(fHistRemovedOrNot);
 973
 974     list->Add(fHistEtNonRemovedProtons);
 975     list->Add(fHistEtNonRemovedAntiProtons);
 976     list->Add(fHistEtNonRemovedPiPlus);
 977     list->Add(fHistEtNonRemovedPiMinus);
 978     list->Add(fHistEtNonRemovedKaonPlus);
 979     list->Add(fHistEtNonRemovedKaonMinus);
 980     list->Add(fHistEtNonRemovedK0s);
 981     list->Add(fHistEtNonRemovedLambdas);
 982     list->Add(fHistEtNonRemovedElectrons);
 983     list->Add(fHistEtNonRemovedPositrons);
 984     list->Add(fHistEtNonRemovedMuPlus);
 985     list->Add(fHistEtNonRemovedMuMinus);
 986     list->Add(fHistEtNonRemovedNeutrons);
 987     list->Add(fHistEtNonRemovedAntiNeutrons);
 988     list->Add(fHistEtNonRemovedGammas);
 989     list->Add(fHistEtNonRemovedGammasFromPi0);
 990
 991     list->Add(fHistEtRemovedGammas);
 992     list->Add(fHistEtRemovedNeutrons);
 993     list->Add(fHistEtRemovedAntiNeutrons);
 994
 995
 996     list->Add(fHistMultNonRemovedProtons);
 997     list->Add(fHistMultNonRemovedAntiProtons);
 998     list->Add(fHistMultNonRemovedPiPlus);
 999     list->Add(fHistMultNonRemovedPiMinus);
1000     list->Add(fHistMultNonRemovedKaonPlus);
1001     list->Add(fHistMultNonRemovedKaonMinus);
1002     list->Add(fHistMultNonRemovedK0s);
1003     list->Add(fHistMultNonRemovedLambdas);
1004     list->Add(fHistMultNonRemovedElectrons);
1005     list->Add(fHistMultNonRemovedPositrons);
1006     list->Add(fHistMultNonRemovedMuPlus);
1007     list->Add(fHistMultNonRemovedMuMinus);
1008     list->Add(fHistMultNonRemovedNeutrons);
1009     list->Add(fHistMultNonRemovedAntiNeutrons);
1010     list->Add(fHistMultNonRemovedGammas);
1011
1012     list->Add(fHistMultRemovedGammas);
1013     list->Add(fHistMultRemovedNeutrons);
1014     list->Add(fHistMultRemovedAntiNeutrons);
1015
1016     list->Add(fHistTrackMultvsNonRemovedCharged);
1017     list->Add(fHistTrackMultvsNonRemovedNeutral);
1018     list->Add(fHistTrackMultvsRemovedGamma);
1019
1020     list->Add(fHistClusterMultvsNonRemovedCharged);
1021     list->Add(fHistClusterMultvsNonRemovedNeutral);
1022     list->Add(fHistClusterMultvsRemovedGamma);
1023
1024     //list->Add(fHistDecayVertexNonRemovedCharged);
1025     //list->Add(fHistDecayVertexNonRemovedNeutral);
1026     //list->Add(fHistDecayVertexRemovedCharged);
1027     //list->Add(fHistDecayVertexRemovedNeutral);
1028
1029     list->Add(fHistDxDzNonRemovedCharged);
1030     list->Add(fHistDxDzRemovedCharged);
1031     list->Add(fHistDxDzNonRemovedNeutral);
1032     list->Add(fHistDxDzRemovedNeutral);
1033
1034     list->Add(fHistPiPlusMult);
1035     list->Add(fHistPiMinusMult);
1036     list->Add(fHistPiZeroMult);
1037     list->Add(fHistPiPlusMultAcc);
1038     list->Add(fHistPiMinusMultAcc);
1039     list->Add(fHistPiZeroMultAcc);
1040
1041 }
1042
1043
1044 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1045 {
1046     //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
1047     AliESDtrack *esdTrack = new AliESDtrack(part);
1048     // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1049
1050     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1051
1052     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1053
1054     bool status = prop &&
1055                   TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
1056                   esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
1057                   esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
1058     delete esdTrack;
1059
1060     return status;
1061 }
1062
1063 void AliAnalysisEtMonteCarlo::FillHistograms()
1064 { // let base class fill its histograms, and us fill the local ones
1065     AliAnalysisEt::FillHistograms();
1066     //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1067
1068     fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1069     fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1070     fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1071     fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
1072     fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0s, fCentClass);
1073     fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1074     fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1075     fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1076     fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1077     fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1078     fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1079     fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1080     fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1081     fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1082     fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1083     fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1084
1085     fHistEtRemovedGammas->Fill(fEtRemovedGammas, fCentClass);
1086     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1087     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1088
1089     fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1090     fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1091     fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1092     fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1093     fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
1094     fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1095     fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1096     fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1097     fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1098     fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1099     fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1100     fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1101     fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1102     fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1103     fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1104
1105     fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1106     fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1107     fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1108
1109     fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1110                                             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1111                                             +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1112                                             +fMultNonRemovedProtons);
1113
1114     fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1115                                             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
1116
1117     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1118                                        fMultRemovedGammas);
1119
1120     fHistClusterMultvsNonRemovedCharged->Fill(fNeutralMultiplicity,
1121             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1122             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1123             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1124
1125     fHistClusterMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1126             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
1127
1128     fHistClusterMultvsRemovedGamma->Fill(fTrackMultInAcc,
1129                                          fMultRemovedGammas);
1130
1131 }
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147