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