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