Changing kAll to kAllP in order to avoid clash with AliDecayer.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderInputSimPrep.cxx
CommitLineData
f7d5860b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
ee6b678f 17/* $Id$ */
f7d5860b 18
19//_________________________________________________________________________
20// Class for JetFinder Input preparation from simulated data
21//
22//*-- Author: Mark Horner (LBL/UCT)
23//
24
25
26
27#include "AliEMCALJetFinderInputSimPrep.h"
28
29#include <TParticle.h>
30#include <TParticlePDG.h>
31#include <TPDGCode.h>
32#include <TFile.h>
33#include <TTree.h>
34#include <TObjectTable.h>
35#include <TMath.h>
36
37
38#include "AliRun.h"
39#include "AliMagF.h"
40#include "AliEMCALFast.h"
41#include "AliEMCAL.h"
42#include "AliEMCALHit.h"
43#include "AliEMCALGeometry.h"
44#include "AliEMCALGetter.h"
b3c7d6bc 45#include "AliEMCALLoader.h"
f7d5860b 46#include "AliGenEventHeader.h"
47#include "AliGenPythiaEventHeader.h"
48#include "AliGenerator.h"
49#include "AliHeader.h"
5d12ce38 50#include "AliMC.h"
f7d5860b 51
52
b3c7d6bc 53
f7d5860b 54ClassImp(AliEMCALJetFinderInputSimPrep)
55
56AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
57{
44f59d68 58 // Default constructor
f7d5860b 59if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
60
61 fDebug = 0;
62 fSmearType = kSmearEffic;
63 fEMCALType = kHits;
64 fTrackType = kCharged;
65 fEfficiency = 0.90;
66
67}
68AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
69{
70if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
71
72}
73
74void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
75{
44f59d68 76 // Method to reset data
f7d5860b 77if (fDebug > 1) Info("Reset","Beginning Reset");
78 switch (resettype){
79
80 case kResetData:
81 fInputObject.Reset(resettype);
82 break;
83 case kResetTracks:
84 fInputObject.Reset(resettype);
85 break;
86 case kResetDigits:
87 fInputObject.Reset(resettype);
88 break;
89 case kResetParameters:
90 fSmearType = kSmearEffic;
91 fEMCALType = kHits;
92 fTrackType = kCharged;
93 break;
94 case kResetAll:
95 fSmearType = kSmearEffic;
96 fEMCALType = kHits;
97 fTrackType = kCharged;
98 fInputObject.Reset(resettype);
99 break;
100 case kResetPartons:
101 Warning("FillFromFile", "kResetPartons not implemented") ;
102 break;
103 case kResetParticles:
104 Warning("FillFromFile", "kResetParticles not implemented") ;
105 break;
106 case kResetJets:
107 Warning("FillFromFile", "kResetJets not implemented") ;
108 break;
109 }// end switch
110
111}
112
113Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
114{
f7d5860b 115if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
116 fFileType = filetype;
f7d5860b 117 AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
b3c7d6bc 118 if (fDebug > 1) Info("FillFromFile","Instantiated Getter with %s as the file",filename->Data());
119 gime->Event(EventNumber,"XH") ;
120 if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
121
f7d5860b 122 if ( fEMCALType == kHits ||
123 fEMCALType == kTimeCut ) FillHits();
124 if ( fTrackType != kNoTracks ) FillTracks();
125 if ( fFileType != kData){
126 FillPartons();
f7d5860b 127 }
b3c7d6bc 128/* gime->EmcalLoader()->UnloadRecParticles();
129 gime->EmcalLoader()->UnloadTracks();
130 gime->Reset();*/
f7d5860b 131 return 0;
132}
133
44f59d68 134void AliEMCALJetFinderInputSimPrep::FillHits()
f7d5860b 135{
44f59d68 136 // Fill from the hits to input object from simulation
f7d5860b 137if (fDebug > 1) Info("FillHits","Beginning FillHits");
138
139// Access hit information
140 AliEMCALHit *mHit;
b3c7d6bc 141// AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
142// Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
f7d5860b 143 AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
b3c7d6bc 144 AliEMCALGetter *gime = AliEMCALGetter::Instance();
145
146 // TTree *treeH = AliEMCALGetter::Instance()->TreeH();
147// Int_t ntracks = (Int_t) treeH->GetEntries();
f7d5860b 148//
149// Loop over tracks
150//
f7d5860b 151 Double_t etH = 0.0;
152
b3c7d6bc 153/* for (Int_t track=0; track<ntracks;track++) {
f7d5860b 154 gAlice->ResetHits();
b3c7d6bc 155 nbytes += treeH->GetEvent(track);*/
f7d5860b 156//
157// Loop over hits
158//
b3c7d6bc 159 for(Int_t i =0; i < gime->Hits()->GetEntries() ;i++)
f7d5860b 160 {
b3c7d6bc 161 mHit = gime->Hit(i);
f7d5860b 162 if (fEMCALType == kTimeCut &&
163 (mHit->GetTime()>fTimeCut) ) continue;
164 Float_t x = mHit->X(); // x-pos of hit
165 Float_t y = mHit->Y(); // y-pos
166 Float_t z = mHit->Z(); // z-pos
167 Float_t eloss = mHit->GetEnergy(); // deposited energy
168
169 Float_t r = TMath::Sqrt(x*x+y*y);
170 Float_t theta = TMath::ATan2(r,z);
171 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
172 Float_t phi = TMath::ATan2(y,x);
5255a1ce 173 etH = eloss*TMath::Sin(theta);
f7d5860b 174 if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
175 if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
176 {
177 if (fDebug >5)
178 {
179 Error("FillHits","Hit not inside EMCAL!");
180 mHit->Dump();
181 }
182 continue;
183 }
184 fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
185 } // Hit Loop
b3c7d6bc 186// } // Track Loop
f7d5860b 187
188
189}
44f59d68 190void AliEMCALJetFinderInputSimPrep::FillTracks()
f7d5860b 191{
44f59d68 192 // Fill from particles simulating a TPC to input object from simulation
f7d5860b 193
194 if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
195
196 TParticlePDG* pdgP = 0;
44f59d68 197 TParticle *mPart;
f7d5860b 198 Int_t npart = (gAlice->GetHeader())->GetNprimary();
b3c7d6bc 199 if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
a7aa0166 200 //AliEMCALGetter *gime = AliEMCALGetter::Instance();
f7d5860b 201 Float_t bfield,rEMCAL;
b3c7d6bc 202/* for (Int_t i =0; i<gime->NPrimaries(); i++)
203 {
204 cout<<gime->Primary(i)->GetFirstMother()<<endl;
205 if (gime->Primary(i)->GetFirstMother() == -1) {
206 counterforus++;
207 gime->Primary(i)->Dump();
208 }
209 }*/
f7d5860b 210
211 if (fDebug > 1) Info("FillTracks","Defining the geometry");
212
213 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
214 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
215 fEtaMax = geom->GetArm1EtaMax();
216 fEtaMin = geom->GetArm1EtaMin();
217 fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
218 fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
219
220 if (fDebug > 1) Info("FillTracks","Starting particle loop");
221
b3c7d6bc 222 if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
f7d5860b 223 for (Int_t part = 0; part < npart; part++) {
44f59d68 224 mPart = gAlice->GetMCApp()->Particle(part);
f7d5860b 225 //if (part%10) gObjectTable->Print();
44f59d68 226 pdgP = mPart->GetPDG();
f7d5860b 227
228 if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
229
230 if (fFileType == kPythia) {
44f59d68 231 if (mPart->GetStatusCode() != 1) continue;
f7d5860b 232 } else if (fFileType == kHijing) {
44f59d68 233 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
f7d5860b 234 }
235
44f59d68 236 if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
f7d5860b 237 if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
238
b3c7d6bc 239 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
240 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
f7d5860b 241
242/*
243 {kProton, kProtonBar, kElectron, kPositron,
244 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
245 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
246 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
247 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
248 0,kNuMu,kNuMuBar
249*/
250
251 if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
252
253 if ((fSmearType == kEfficiency ||
254 fSmearType == kSmearEffic) &&
255 pdgP->Charge()!=0) {
256 if (AliEMCALFast::RandomReject(fEfficiency)) {
257 continue;
258 }
259 }
260
261 bfield = gAlice->Field()->SolenoidField();
262 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
44f59d68 263 Float_t rB = 3335.6 * mPart->Pt() / bfield; // [cm] (case of |charge|=1)
f7d5860b 264 if (2.*rB < rEMCAL) continue; // track curls away
265
266 //if (part%10) gObjectTable->Print();
267 switch(fTrackType)
268 {
269
58c93485 270 case kAllP: // All Stable particles to be included
f7d5860b 271 if (fDebug > 5) Info("FillTracks","Storing track");
272 if (fSmearType == kSmear ||
273 fSmearType == kSmearEffic ){
44f59d68 274 Smear(mPart);/*
f7d5860b 275 TParticle *tmp = Smear(MPart);
276 fInputObject.AddTrack(Smear(MPart));
277 delete tmp;*/
278 }else{
44f59d68 279 fInputObject.AddTrack(*mPart);
f7d5860b 280 }
281 break;
282 case kEM: // All Electromagnetic particles
44f59d68 283 if (mPart->GetPdgCode() == kElectron ||
284 mPart->GetPdgCode() == kMuonPlus ||
285 mPart->GetPdgCode() == kMuonMinus ||
286 mPart->GetPdgCode() == kPositron ){
f7d5860b 287 if (fDebug > 5) Info("FillTracks","Storing electron or positron");
288 if (fSmearType == kSmear ||
289 fSmearType == kSmearEffic ){
44f59d68 290 Smear(mPart);/*
f7d5860b 291 TParticle *tmp = Smear(MPart);
292 fInputObject.AddTrack(tmp);
293 delete tmp;*/
294 }else{
44f59d68 295 fInputObject.AddTrack(*mPart);
f7d5860b 296 }
297 }
44f59d68 298 if ( mPart->GetPdgCode() == kGamma ){
299 fInputObject.AddTrack(*mPart);
f7d5860b 300 if (fDebug > 5) Info("FillTracks","Storing gamma");
301 }
302
303 break;
304 case kCharged: // All Particles with non-zero charge
305 if (pdgP->Charge() != 0) {
306 if (fDebug > 5) Info("FillTracks","Storing charged track");
307 if (fSmearType == kSmear ||
308 fSmearType == kSmearEffic ){
44f59d68 309 Smear(mPart);/*
f7d5860b 310 TParticle *tmp = Smear(MPart);
311 fInputObject.AddTrack(tmp);
312 delete tmp;*/
313 }else{
44f59d68 314 fInputObject.AddTrack(*mPart);
f7d5860b 315 }
316 }
317 break;
318 case kNeutral: // All particles with no charge
319 if (pdgP->Charge() == 0){
44f59d68 320 fInputObject.AddTrack(*mPart);
f7d5860b 321 if (fDebug > 5) Info("FillTracks","Storing neutral");
322 }
323 break;
324 case kHadron: //All hadrons
44f59d68 325 if (mPart->GetPdgCode() != kElectron &&
326 mPart->GetPdgCode() != kPositron &&
327 mPart->GetPdgCode() != kMuonPlus &&
328 mPart->GetPdgCode() != kMuonMinus &&
329 mPart->GetPdgCode() != kGamma )
f7d5860b 330 {
331 if (fDebug > 5) Info("FillTracks","Storing hadron");
332 if (pdgP->Charge() == 0){
44f59d68 333 fInputObject.AddTrack(*mPart);
f7d5860b 334 }else{
335 if (fSmearType == kSmear ||
336 fSmearType == kSmearEffic ){
44f59d68 337 Smear(mPart);/*
f7d5860b 338 TParticle *tmp = Smear(MPart);
339 fInputObject.AddTrack(tmp);
340 delete tmp;*/
341 }else{
44f59d68 342 fInputObject.AddTrack(*mPart);
f7d5860b 343 }
344 }
345 }
346 break;
347 case kChargedHadron: // only charged hadrons
44f59d68 348 if (mPart->GetPdgCode() != kElectron &&
349 mPart->GetPdgCode() != kPositron &&
350 mPart->GetPdgCode() != kGamma &&
351 mPart->GetPdgCode() != kMuonPlus &&
352 mPart->GetPdgCode() != kMuonMinus &&
f7d5860b 353 pdgP->Charge() != 0 ){
354 if (fDebug > 5) Info("FillTracks","Storing charged hadron");
355 if (fSmearType == kSmear ||
356 fSmearType == kSmearEffic ){
44f59d68 357 Smear(mPart);/*
f7d5860b 358 TParticle *tmp = Smear(MPart);
359 fInputObject.AddTrack(tmp);
360 delete tmp;*/
361 }else{
44f59d68 362 fInputObject.AddTrack(*mPart);
f7d5860b 363 }
364 }
365 break;
366 case kNoTracks:
367 break;
ef09619f 368 case kEMChargedPi0:
369 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
370 mPart->GetPdgCode() == kGamma )
371 {
372 if (fDebug > 5) Info("FillTracks","Storing charged track");
373 if (fSmearType == kSmear ||
374 fSmearType == kSmearEffic ){
375 Smear(mPart);/*
376 TParticle *tmp = Smear(MPart);
377 fInputObject.AddTrack(tmp);
378 delete tmp;*/
379 }else{
380 fInputObject.AddTrack(*mPart);
381 }
382 }
383 break;
a4c7a211 384 case kNoNeutronNeutrinoKlong:
385 if ( mPart->GetPdgCode() != kNeutron &&
386 mPart->GetPdgCode() != kNeutronBar &&
387 mPart->GetPdgCode() != kK0Long &&
388 mPart->GetPdgCode() != kNuE &&
389 mPart->GetPdgCode() != kNuEBar &&
390 mPart->GetPdgCode() != kNuMu &&
391 mPart->GetPdgCode() != kNuMuBar &&
392 mPart->GetPdgCode() != kNuTau &&
393 mPart->GetPdgCode() != kNuTauBar )
394 {
395 if (fDebug > 5) Info("FillTracks","Storing charged track");
396 if (fSmearType == kSmear ||
397 fSmearType == kSmearEffic ){
398 Smear(mPart);/*
399 TParticle *tmp = Smear(MPart);
400 fInputObject.AddTrack(tmp);
401 delete tmp;*/
402 }else{
403 fInputObject.AddTrack(*mPart);
404 }
405 }
406 break;
f7d5860b 407 default:
408 break;
44f59d68 409 delete mPart;
f7d5860b 410 } //end of switch
411// Info("FillTracks","After Particle Storage");
412 //if (part%10) gObjectTable->Print();
413
414 } //end of particle loop
415 //gObjectTable->Print();
416}
417
44f59d68 418void AliEMCALJetFinderInputSimPrep::FillPartons()
f7d5860b 419{
44f59d68 420 // Fill partons to
421 // input object from simulation
f7d5860b 422if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
423
424 AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
425 Float_t triggerJetValues[4];
426 AliEMCALParton tempParton;
427
428 if (fFileType == kPythia)
429 {
430 Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
431 if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
432 for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
433 ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
434 TLorentzVector tempLParton;
435 tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
436 tempParton.SetEnergy(tempLParton.Energy());
437 tempParton.SetEta(tempLParton.Eta());
438 tempParton.SetPhi(tempLParton.Phi());
439 FillPartonTracks(&tempParton);
440 fInputObject.AddParton(&tempParton);
441 }
442
443 }
444}
445
44f59d68 446void AliEMCALJetFinderInputSimPrep::FillParticles()
f7d5860b 447{
44f59d68 448 // Fill particles to input object from simulation
f7d5860b 449if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
450
451 Int_t npart = (gAlice->GetHeader())->GetNprimary();
452 TParticlePDG* pdgP = 0;
453
454 for (Int_t part = 0; part < npart; part++) {
44f59d68 455 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
456 pdgP = mPart->GetPDG();
f7d5860b 457
458 if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
459
460 if (fFileType == kPythia) {
44f59d68 461 if (mPart->GetStatusCode() != 1) continue;
f7d5860b 462 } else if (fFileType == kHijing) {
44f59d68 463 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
f7d5860b 464 }
465
466 if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
467
b3c7d6bc 468 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
469 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
f7d5860b 470
471
472/*
473 {kProton, kProtonBar, kElectron, kPositron,
474 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
475 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
476 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
477 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
478 0,kNuMu,kNuMuBar
479*/
480 switch(fTrackType)
481 {
482
58c93485 483 case kAllP: // All Stable particles to be included
f7d5860b 484 if (fDebug > 5) Info("FillParticles","Storing particle");
44f59d68 485 fInputObject.AddParticle(mPart);
f7d5860b 486 break;
487 case kEM: // All Electromagnetic particles
44f59d68 488 if (mPart->GetPdgCode() == kElectron ||
489 mPart->GetPdgCode() == kPositron ||
490 mPart->GetPdgCode() == kGamma){
f7d5860b 491 if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
44f59d68 492 fInputObject.AddParticle(mPart);
f7d5860b 493 }
494 break;
495 case kCharged: // All Particles with non-zero charge
496 if (pdgP->Charge() != 0) {
497 if (fDebug > 5) Info("FillParticles","Storing charged particle");
44f59d68 498 fInputObject.AddParticle(mPart);
f7d5860b 499 }
500 break;
501 case kNeutral: // All particles with no charge
502 if (pdgP->Charge() == 0){
503 if (fDebug > 5) Info("FillParticles","Storing neutral particle");
44f59d68 504 fInputObject.AddParticle(mPart);
f7d5860b 505 }
506 break;
507 case kHadron: //All hadrons
44f59d68 508 if (mPart->GetPdgCode() != kElectron &&
509 mPart->GetPdgCode() != kPositron &&
510 mPart->GetPdgCode() != kGamma )
f7d5860b 511 {
512
513 if (fDebug > 5) Info("FillParticles","Storing hadron");
44f59d68 514 fInputObject.AddParticle(mPart);
f7d5860b 515 }
516 break;
517 case kChargedHadron: // only charged hadrons
44f59d68 518 if (mPart->GetPdgCode() != kElectron &&
519 mPart->GetPdgCode() != kPositron &&
520 mPart->GetPdgCode() != kGamma &&
f7d5860b 521 pdgP->Charge() != 0 ){
522 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
44f59d68 523 fInputObject.AddParticle(mPart);
f7d5860b 524 }
525 break;
526 case kNoTracks:
527 break;
ef09619f 528 case kEMChargedPi0:
529 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
530 mPart->GetPdgCode() == kGamma )
531 {
a4c7a211 532 if (fDebug > 5) Info("FillTracks","Storing kEMChargedPi0 track");
ef09619f 533 if (fSmearType == kSmear ||
534 fSmearType == kSmearEffic ){
535 Smear(mPart);/*
536 TParticle *tmp = Smear(MPart);
537 fInputObject.AddTrack(tmp);
538 delete tmp;*/
539 }else{
540 fInputObject.AddTrack(*mPart);
541 }
542 }
543 break;
a4c7a211 544 case kNoNeutronNeutrinoKlong:
545 if ( mPart->GetPdgCode() != kNeutron &&
546 mPart->GetPdgCode() != kNeutronBar &&
547 mPart->GetPdgCode() != kK0Long &&
548 mPart->GetPdgCode() != kNuE &&
549 mPart->GetPdgCode() != kNuEBar &&
550 mPart->GetPdgCode() != kNuMu &&
551 mPart->GetPdgCode() != kNuMuBar &&
552 mPart->GetPdgCode() != kNuTau &&
553 mPart->GetPdgCode() != kNuTauBar )
554 {
555 if (fDebug > 5) Info("FillTracks","Storing kNoNeutronNeutrinoKlong track");
556 if (fSmearType == kSmear ||
557 fSmearType == kSmearEffic ){
558 Smear(mPart);/*
559 TParticle *tmp = Smear(MPart);
560 fInputObject.AddTrack(tmp);
561 delete tmp;*/
562 }else{
563 fInputObject.AddTrack(*mPart);
564 }
565 }
566 break;
f7d5860b 567 default:
568 break;
569 } //end of switch
570 }// end of loop over particles
571
572}
44f59d68 573void AliEMCALJetFinderInputSimPrep::FillDigits()
f7d5860b 574{
44f59d68 575 // Fill digits to input object
f7d5860b 576
577}
578
579void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
580{
44f59d68 581 // Smear particle momentum
f7d5860b 582if (fDebug > 5) Info("Smear","Beginning Smear");
583
584 Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
585 Float_t tmpt = (tmp/particle->P()) * particle->Pt();
586 if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
587 TLorentzVector tmpvector;
588 tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
589// create a new particle with smearing momentum - and no daughter or parent information and the same
590// vertex
591
592 TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
593 tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
594 particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
595 fInputObject.AddTrack(tmparticle);
596 return;
597}
598
599void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
600{
44f59d68 601 // Populate parton tracks for input distributions
f7d5860b 602 if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
603 Int_t npart = (gAlice->GetHeader())->GetNprimary();
604 Int_t ntracks =0;
605 TParticlePDG *getpdg;
606 TParticle *tempPart;
607 for (Int_t part = 0; part < npart; part++)
608 {
5d12ce38 609 tempPart = gAlice->GetMCApp()->Particle(part);
f7d5860b 610 if (tempPart->GetStatusCode() != 1) continue;
b3c7d6bc 611 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
612 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
f7d5860b 613 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
614 continue;
615 }
616 getpdg = tempPart->GetPDG();
617 if (getpdg->Charge() == 0.0 ) {
618 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
619 continue;
620 }
621 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
622 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
623 }
44f59d68 624 Float_t *energy = new Float_t[ntracks];
625 Float_t *eta = new Float_t[ntracks];
626 Float_t *phi = new Float_t[ntracks];
627 Int_t *pdg = new Int_t[ntracks];
f7d5860b 628 ntracks=0;
629 for (Int_t part = 0; part < npart; part++)
630 {
5d12ce38 631 tempPart = gAlice->GetMCApp()->Particle(part);
f7d5860b 632 if (tempPart->GetStatusCode() != 1) continue;
b3c7d6bc 633 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
634 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
f7d5860b 635 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
636 continue;
637 }
638 if (tempPart->GetStatusCode() != 1) continue;
639 getpdg = tempPart->GetPDG();
640 if (getpdg->Charge() == 0.0 ) {
641 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
642 continue;
643 }
644 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
645 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
646 {
44f59d68 647 energy[ntracks] = tempPart->Pt();
648 eta[ntracks] = tempPart->Eta();
649 phi[ntracks] = tempPart->Phi();
650 pdg[ntracks] = tempPart->GetPdgCode();
f7d5860b 651 ntracks++;
652 }
653 }
44f59d68 654 parton->SetTrackList(ntracks,energy,eta,phi,pdg);
b3c7d6bc 655 delete[] energy;
656 delete[] eta;
657 delete[] phi;
658 delete[] pdg;
f7d5860b 659}
660