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