]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALJetFinderInputSimPrep.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fgMCEvGen changed to fMCEvGen.
[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"
45#include "AliGenEventHeader.h"
46#include "AliGenPythiaEventHeader.h"
47#include "AliGenerator.h"
48#include "AliHeader.h"
49
50
51ClassImp(AliEMCALJetFinderInputSimPrep)
52
53AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
54{
55if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
56
57 fDebug = 0;
58 fSmearType = kSmearEffic;
59 fEMCALType = kHits;
60 fTrackType = kCharged;
61 fEfficiency = 0.90;
62
63}
64AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
65{
66if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
67
68}
69
70void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
71{
72if (fDebug > 1) Info("Reset","Beginning Reset");
73 switch (resettype){
74
75 case kResetData:
76 fInputObject.Reset(resettype);
77 break;
78 case kResetTracks:
79 fInputObject.Reset(resettype);
80 break;
81 case kResetDigits:
82 fInputObject.Reset(resettype);
83 break;
84 case kResetParameters:
85 fSmearType = kSmearEffic;
86 fEMCALType = kHits;
87 fTrackType = kCharged;
88 break;
89 case kResetAll:
90 fSmearType = kSmearEffic;
91 fEMCALType = kHits;
92 fTrackType = kCharged;
93 fInputObject.Reset(resettype);
94 break;
95 case kResetPartons:
96 Warning("FillFromFile", "kResetPartons not implemented") ;
97 break;
98 case kResetParticles:
99 Warning("FillFromFile", "kResetParticles not implemented") ;
100 break;
101 case kResetJets:
102 Warning("FillFromFile", "kResetJets not implemented") ;
103 break;
104 }// end switch
105
106}
107
108Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
109{
110 // gObjectTable->Print();
111 // Test that file exists - Getter doesn't like bogus filenames
112if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
113 fFileType = filetype;
114 TFile file(filename->Data());
115 if (!file.IsOpen()){
116 Error("FillFromFile","Could not open file in FillFromFile");
117 return -1;
118 }
119 file.Close();
120 /*
121 gAlice = static_cast<AliRun *>(file.Get("gAlice"));
122 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
123 if (gAlice->GetEvent(EventNumber) < 0){
124 Error("FillFromFile","Could not open event in FillFromFile");
125 file.Close();
126 return -1;
127 }*/
128 AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
129 gime->Event(EventNumber) ;
130
131 if ( fEMCALType == kHits ||
132 fEMCALType == kTimeCut ) FillHits();
133 if ( fTrackType != kNoTracks ) FillTracks();
134 if ( fFileType != kData){
135 FillPartons();
136// FillParticles();
137 }
138 //gime->CloseFile();
139// gObjectTable->Print();
140 delete gime;
141 return 0;
142}
143
144void AliEMCALJetFinderInputSimPrep::FillHits() // Fill from the hits to input object from simulation
145{
146if (fDebug > 1) Info("FillHits","Beginning FillHits");
147
148// Access hit information
149 AliEMCALHit *mHit;
150 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
151 Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
152 AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
153// Float_t samplingF = geom->GetSampling();
154 Float_t samplingF = 11.6;
155 Info("AliEMCALJetFinderInputSimPrep","Sampling fraction is %f",samplingF);
156 TTree *treeH = AliEMCALGetter::Instance()->TreeH();
157 Int_t ntracks = (Int_t) treeH->GetEntries();
158//
159// Loop over tracks
160//
161 Int_t nbytes = 0;
162 Double_t etH = 0.0;
163
164 for (Int_t track=0; track<ntracks;track++) {
165 gAlice->ResetHits();
166 nbytes += treeH->GetEvent(track);
167//
168// Loop over hits
169//
170 for(mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
171 mHit;
172 mHit=(AliEMCALHit*) pEMCAL->NextHit())
173 {
174 if (fEMCALType == kTimeCut &&
175 (mHit->GetTime()>fTimeCut) ) continue;
176 Float_t x = mHit->X(); // x-pos of hit
177 Float_t y = mHit->Y(); // y-pos
178 Float_t z = mHit->Z(); // z-pos
179 Float_t eloss = mHit->GetEnergy(); // deposited energy
180
181 Float_t r = TMath::Sqrt(x*x+y*y);
182 Float_t theta = TMath::ATan2(r,z);
183 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
184 Float_t phi = TMath::ATan2(y,x);
185 etH = samplingF*eloss*TMath::Sin(theta);
186 if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
187 if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
188 {
189 if (fDebug >5)
190 {
191 Error("FillHits","Hit not inside EMCAL!");
192 mHit->Dump();
193 }
194 continue;
195 }
196 fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
197 } // Hit Loop
198 } // Track Loop
199
200
201}
202void AliEMCALJetFinderInputSimPrep::FillTracks() // Fill from particles simulating a TPC to input object from simulation
203{
204
205 if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
206
207 TParticlePDG* pdgP = 0;
208 TParticle *MPart;
209 Int_t npart = (gAlice->GetHeader())->GetNprimary();
210 Float_t bfield,rEMCAL;
211
212 if (fDebug > 1) Info("FillTracks","Defining the geometry");
213
214 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
215 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
216 fEtaMax = geom->GetArm1EtaMax();
217 fEtaMin = geom->GetArm1EtaMin();
218 fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
219 fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
220
221 if (fDebug > 1) Info("FillTracks","Starting particle loop");
222
223 for (Int_t part = 0; part < npart; part++) {
224 MPart = gAlice->Particle(part);
225 //if (part%10) gObjectTable->Print();
226 pdgP = MPart->GetPDG();
227
228 if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
229
230 if (fFileType == kPythia) {
231 if (MPart->GetStatusCode() != 1) continue;
232 } else if (fFileType == kHijing) {
233 if (MPart->GetFirstDaughter() >= 0 && MPart->GetFirstDaughter() < npart) continue;
234 }
235
236 if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",MPart->Eta(),MPart->Phi());
237 if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
238
239 if (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin) continue;
240 if (MPart->Phi() > fPhiMax || MPart->Phi() < fPhiMin ) continue;
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();
263 Float_t rB = 3335.6 * MPart->Pt() / bfield; // [cm] (case of |charge|=1)
264 if (2.*rB < rEMCAL) continue; // track curls away
265
266 //if (part%10) gObjectTable->Print();
267 switch(fTrackType)
268 {
269
270 case kAll: // All Stable particles to be included
271 if (fDebug > 5) Info("FillTracks","Storing track");
272 if (fSmearType == kSmear ||
273 fSmearType == kSmearEffic ){
274 Smear(MPart);/*
275 TParticle *tmp = Smear(MPart);
276 fInputObject.AddTrack(Smear(MPart));
277 delete tmp;*/
278 }else{
279 fInputObject.AddTrack(*MPart);
280 }
281 break;
282 case kEM: // All Electromagnetic particles
9411bca0 283 if (MPart->GetPdgCode() == kElectron ||
284 MPart->GetPdgCode() == kMuonPlus ||
285 MPart->GetPdgCode() == kMuonMinus ||
f7d5860b 286 MPart->GetPdgCode() == kPositron ){
287 if (fDebug > 5) Info("FillTracks","Storing electron or positron");
288 if (fSmearType == kSmear ||
289 fSmearType == kSmearEffic ){
290 Smear(MPart);/*
291 TParticle *tmp = Smear(MPart);
292 fInputObject.AddTrack(tmp);
293 delete tmp;*/
294 }else{
295 fInputObject.AddTrack(*MPart);
296 }
297 }
298 if ( MPart->GetPdgCode() == kGamma ){
299 fInputObject.AddTrack(*MPart);
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 ){
309 Smear(MPart);/*
310 TParticle *tmp = Smear(MPart);
311 fInputObject.AddTrack(tmp);
312 delete tmp;*/
313 }else{
314 fInputObject.AddTrack(*MPart);
315 }
316 }
317 break;
318 case kNeutral: // All particles with no charge
319 if (pdgP->Charge() == 0){
320 fInputObject.AddTrack(*MPart);
321 if (fDebug > 5) Info("FillTracks","Storing neutral");
322 }
323 break;
324 case kHadron: //All hadrons
67a57da8 325 if (MPart->GetPdgCode() != kElectron &&
326 MPart->GetPdgCode() != kPositron &&
327 MPart->GetPdgCode() != kMuonPlus &&
328 MPart->GetPdgCode() != kMuonMinus &&
f7d5860b 329 MPart->GetPdgCode() != kGamma )
330 {
331 if (fDebug > 5) Info("FillTracks","Storing hadron");
332 if (pdgP->Charge() == 0){
333 fInputObject.AddTrack(*MPart);
334 }else{
335 if (fSmearType == kSmear ||
336 fSmearType == kSmearEffic ){
337 Smear(MPart);/*
338 TParticle *tmp = Smear(MPart);
339 fInputObject.AddTrack(tmp);
340 delete tmp;*/
341 }else{
342 fInputObject.AddTrack(*MPart);
343 }
344 }
345 }
346 break;
347 case kChargedHadron: // only charged hadrons
67a57da8 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 ){
357 Smear(MPart);/*
358 TParticle *tmp = Smear(MPart);
359 fInputObject.AddTrack(tmp);
360 delete tmp;*/
361 }else{
362 fInputObject.AddTrack(*MPart);
363 }
364 }
365 break;
366 case kNoTracks:
367 break;
368 default:
369 break;
370 delete MPart;
371 } //end of switch
372// Info("FillTracks","After Particle Storage");
373 //if (part%10) gObjectTable->Print();
374
375 } //end of particle loop
376 //gObjectTable->Print();
377}
378
379void AliEMCALJetFinderInputSimPrep::FillPartons() // Fill partons to input object from simulation
380{
381if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
382
383 AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
384 Float_t triggerJetValues[4];
385 AliEMCALParton tempParton;
386
387 if (fFileType == kPythia)
388 {
389 Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
390 if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
391 for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
392 ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
393 TLorentzVector tempLParton;
394 tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
395 tempParton.SetEnergy(tempLParton.Energy());
396 tempParton.SetEta(tempLParton.Eta());
397 tempParton.SetPhi(tempLParton.Phi());
398 FillPartonTracks(&tempParton);
399 fInputObject.AddParton(&tempParton);
400 }
401
402 }
403}
404
405void AliEMCALJetFinderInputSimPrep::FillParticles() // Fill particles to input object from simulation
406{
407if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
408
409 Int_t npart = (gAlice->GetHeader())->GetNprimary();
410 TParticlePDG* pdgP = 0;
411
412 for (Int_t part = 0; part < npart; part++) {
413 TParticle *MPart = gAlice->Particle(part);
414 pdgP = MPart->GetPDG();
415
416 if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
417
418 if (fFileType == kPythia) {
419 if (MPart->GetStatusCode() != 1) continue;
420 } else if (fFileType == kHijing) {
421 if (MPart->GetFirstDaughter() >= 0 && MPart->GetFirstDaughter() < npart) continue;
422 }
423
424 if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
425
426 if (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin) continue;
427 if (MPart->Phi() > fPhiMax || MPart->Phi() < fPhiMin ) continue;
428
429
430/*
431 {kProton, kProtonBar, kElectron, kPositron,
432 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
433 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
434 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
435 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
436 0,kNuMu,kNuMuBar
437*/
438 switch(fTrackType)
439 {
440
441 case kAll: // All Stable particles to be included
442 if (fDebug > 5) Info("FillParticles","Storing particle");
443 fInputObject.AddParticle(MPart);
444 break;
445 case kEM: // All Electromagnetic particles
446 if (MPart->GetPdgCode() == kElectron ||
447 MPart->GetPdgCode() == kPositron ||
448 MPart->GetPdgCode() == kGamma){
449 if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
450 fInputObject.AddParticle(MPart);
451 }
452 break;
453 case kCharged: // All Particles with non-zero charge
454 if (pdgP->Charge() != 0) {
455 if (fDebug > 5) Info("FillParticles","Storing charged particle");
456 fInputObject.AddParticle(MPart);
457 }
458 break;
459 case kNeutral: // All particles with no charge
460 if (pdgP->Charge() == 0){
461 if (fDebug > 5) Info("FillParticles","Storing neutral particle");
462 fInputObject.AddParticle(MPart);
463 }
464 break;
465 case kHadron: //All hadrons
466 if (MPart->GetPdgCode() != kElectron &&
467 MPart->GetPdgCode() != kPositron &&
468 MPart->GetPdgCode() != kGamma )
469 {
470
471 if (fDebug > 5) Info("FillParticles","Storing hadron");
472 fInputObject.AddParticle(MPart);
473 }
474 break;
475 case kChargedHadron: // only charged hadrons
476 if (MPart->GetPdgCode() != kElectron &&
477 MPart->GetPdgCode() != kPositron &&
478 MPart->GetPdgCode() != kGamma &&
479 pdgP->Charge() != 0 ){
480 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
481 fInputObject.AddParticle(MPart);
482 }
483 break;
484 case kNoTracks:
485 break;
486 default:
487 break;
488 } //end of switch
489 }// end of loop over particles
490
491}
492void AliEMCALJetFinderInputSimPrep::FillDigits() // Fill digits to input object
493{
494
495}
496
497void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
498{
499if (fDebug > 5) Info("Smear","Beginning Smear");
500
501 Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
502 Float_t tmpt = (tmp/particle->P()) * particle->Pt();
503 if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
504 TLorentzVector tmpvector;
505 tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
506// create a new particle with smearing momentum - and no daughter or parent information and the same
507// vertex
508
509 TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
510 tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
511 particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
512 fInputObject.AddTrack(tmparticle);
513 return;
514}
515
516void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
517{
518
519 if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
520 Int_t npart = (gAlice->GetHeader())->GetNprimary();
521 Int_t ntracks =0;
522 TParticlePDG *getpdg;
523 TParticle *tempPart;
524 for (Int_t part = 0; part < npart; part++)
525 {
526 tempPart = gAlice->Particle(part);
527 if (tempPart->GetStatusCode() != 1) continue;
528 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
529 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
530 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
531 continue;
532 }
533 getpdg = tempPart->GetPDG();
534 if (getpdg->Charge() == 0.0 ) {
535 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
536 continue;
537 }
538 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
539 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
540 }
541 Float_t *Energy = new Float_t[ntracks];
542 Float_t *Eta = new Float_t[ntracks];
543 Float_t *Phi = new Float_t[ntracks];
544 Int_t *Pdg = new Int_t[ntracks];
545 ntracks=0;
546 for (Int_t part = 0; part < npart; part++)
547 {
548 tempPart = gAlice->Particle(part);
549 if (tempPart->GetStatusCode() != 1) continue;
550 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
551 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
552 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
553 continue;
554 }
555 if (tempPart->GetStatusCode() != 1) continue;
556 getpdg = tempPart->GetPDG();
557 if (getpdg->Charge() == 0.0 ) {
558 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
559 continue;
560 }
561 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
562 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
563 {
564 Energy[ntracks] = tempPart->Pt();
565 Eta[ntracks] = tempPart->Eta();
566 Phi[ntracks] = tempPart->Phi();
567 Pdg[ntracks] = tempPart->GetPdgCode();
568 ntracks++;
569 }
570 }
571 parton->SetTrackList(ntracks,Energy,Eta,Phi,Pdg);
572}
573