Line data Source code
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 : /* $Id$ */
17 :
18 : //_________________________________________________________________________
19 : //*-- Implementation version v1 of EMCAL Manager class
20 : //*-- An object of this class does not produce digits
21 : //*-- It is the one to use if you do want to produce outputs in TREEH
22 : //*--
23 : //*-- Author: Sahal Yacoob (LBL /UCT)
24 : //*-- : Jennifer Klay (LBL)
25 : // This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis.
26 : // This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL.
27 :
28 : // 15/02/2002 .... Yves Schutz
29 : // 1. fSamplingFraction and fLayerToPreshowerRatio have been removed
30 : // 2. Timing signal is collected and added to hit
31 :
32 : // --- ROOT system ---
33 : #include <TClonesArray.h>
34 : #include <TParticle.h>
35 : #include <TVirtualMC.h>
36 :
37 : // --- Standard library ---
38 :
39 : // --- AliRoot header files ---
40 : #include "AliEMCALv1.h"
41 : #include "AliEMCALHit.h"
42 : #include "AliEMCALGeometry.h"
43 : #include "AliRun.h"
44 : #include "AliMC.h"
45 : #include "AliStack.h"
46 : #include <iostream>
47 :
48 : using std::cout;
49 : using std::endl;
50 42 : ClassImp(AliEMCALv1)
51 :
52 :
53 : //______________________________________________________________________
54 : AliEMCALv1::AliEMCALv1()
55 12 : : AliEMCALv0(),
56 12 : fCurPrimary(-1),
57 12 : fCurParent(-1),
58 12 : fCurTrack(-1),
59 12 : fTimeCut(30e-09)
60 36 : {
61 : // default ctor
62 12 : }
63 :
64 : //______________________________________________________________________
65 : AliEMCALv1::AliEMCALv1(const char *name, const char *title,
66 : const Bool_t checkGeoAndRun)
67 1 : : AliEMCALv0(name,title,checkGeoAndRun),
68 1 : fCurPrimary(-1),
69 1 : fCurParent(-1),
70 1 : fCurTrack(-1),
71 1 : fTimeCut(30e-09)
72 2 : {
73 : // Standard Creator.
74 :
75 3 : fHits= new TClonesArray("AliEMCALHit",1000);
76 : // gAlice->GetMCApp()->AddHitList(fHits); // 20-dec-04 - advice of Andreas
77 :
78 1 : fNhits = 0;
79 1 : fIshunt = 2; // All hits are associated with particles entering the calorimeter
80 1 : }
81 :
82 : //______________________________________________________________________
83 26 : AliEMCALv1::~AliEMCALv1(){
84 : // dtor
85 :
86 13 : if ( fHits ) {
87 1 : fHits->Clear();
88 2 : delete fHits;
89 1 : fHits = 0;
90 1 : }
91 13 : }
92 :
93 : //______________________________________________________________________
94 : void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy,
95 : Int_t id, Float_t * hits,Float_t * p){
96 : // Add a hit to the hit list.
97 : // An EMCAL hit is the sum of all hits in a tower section
98 : // originating from the same entering particle
99 : Int_t hitCounter=0;
100 :
101 : AliEMCALHit *newHit=0;
102 : AliEMCALHit *curHit=0;
103 : Bool_t deja = kFALSE;
104 :
105 0 : newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
106 0 : for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) {
107 0 : curHit = (AliEMCALHit*) (*fHits)[hitCounter];
108 : // We add hits with the same tracknumber, while GEANT treats
109 : // primaries succesively
110 0 : if(curHit->GetPrimary() != primary)
111 : break;
112 0 : if( *curHit == *newHit ) {
113 0 : *curHit = *curHit + *newHit;
114 : deja = kTRUE;
115 0 : } // end if
116 : } // end for hitCounter
117 :
118 0 : if ( !deja ) {
119 0 : new((*fHits)[fNhits]) AliEMCALHit(*newHit);
120 0 : fNhits++;
121 0 : } // end if
122 :
123 0 : delete newHit;
124 0 : }
125 : //______________________________________________________________________
126 : void AliEMCALv1::StepManager(void){
127 : // Accumulates hits as long as the track stays in a tower
128 :
129 0 : Int_t id[2]={0,0}; // (phi, Eta) indices
130 : // position wrt MRS and energy deposited
131 0 : Float_t xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
132 0 : Float_t pmom[4]={0.,0.,0.,0.};
133 0 : TLorentzVector pos; // Lorentz vector of the track current position.
134 0 : TLorentzVector mom; // Lorentz vector of the track current momentum.
135 0 : Int_t tracknumber = gAlice->GetMCApp()->GetCurrentTrackNumber();
136 : static Float_t ienergy = 0;
137 0 : Int_t copy = 0;
138 :
139 0 : AliEMCALGeometry * geom = GetGeometry() ;
140 :
141 : TParticle *part=0;
142 : Int_t parent=-1;
143 :
144 0 : static Int_t idXPHI = TVirtualMC::GetMC()->VolId("XPHI");
145 0 : if(TVirtualMC::GetMC()->CurrentVolID(copy) == idXPHI ) { // We are in a Scintillator Layer
146 : Float_t depositedEnergy ;
147 :
148 0 : if( ((depositedEnergy = TVirtualMC::GetMC()->Edep()) > 0.) && (TVirtualMC::GetMC()->TrackTime() < fTimeCut)){// Track is inside a scintillator and deposits some energy
149 0 : if (fCurPrimary==-1)
150 0 : fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber);
151 :
152 0 : if (fCurParent==-1 || tracknumber != fCurTrack) {
153 : // Check parentage
154 : //Int_t parent=tracknumber;
155 : parent=tracknumber;
156 0 : if (fCurParent != -1) {
157 0 : while (parent != fCurParent && parent != -1) {
158 : //TParticle *part=gAlice->GetMCApp()->Particle(parent);
159 0 : part=gAlice->GetMCApp()->Particle(parent);
160 0 : parent=part->GetFirstMother();
161 : }
162 : }
163 0 : if (fCurParent==-1 || parent==-1) {
164 : //Int_t parent=tracknumber;
165 : //TParticle *part=gAlice->GetMCApp()->Particle(parent);
166 : parent=tracknumber;
167 0 : part=gAlice->GetMCApp()->Particle(parent);
168 0 : while (parent != -1 && geom->IsInEMCAL(part->Vx(),part->Vy(),part->Vz())) {
169 0 : parent=part->GetFirstMother();
170 0 : if (parent!=-1)
171 0 : part=gAlice->GetMCApp()->Particle(parent);
172 : }
173 0 : fCurParent=parent;
174 0 : if (fCurParent==-1)
175 0 : Error("StepManager","Cannot find parent");
176 : else {
177 : //TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
178 0 : part=gAlice->GetMCApp()->Particle(fCurParent);
179 0 : ienergy = part->Energy();
180 : }
181 0 : while (parent != -1) {
182 0 : part=gAlice->GetMCApp()->Particle(parent);
183 0 : part->SetBit(kKeepBit);
184 0 : parent=part->GetFirstMother();
185 : }
186 : }
187 0 : fCurTrack=tracknumber;
188 0 : }
189 0 : TVirtualMC::GetMC()->TrackPosition(pos);
190 0 : xyzte[0] = pos[0];
191 0 : xyzte[1] = pos[1];
192 0 : xyzte[2] = pos[2];
193 0 : xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
194 :
195 0 : TVirtualMC::GetMC()->TrackMomentum(mom);
196 0 : pmom[0] = mom[0];
197 0 : pmom[1] = mom[1];
198 0 : pmom[2] = mom[2];
199 0 : pmom[3] = mom[3];
200 :
201 0 : TVirtualMC::GetMC()->CurrentVolOffID(1, id[0]); // get the POLY copy number;
202 0 : TVirtualMC::GetMC()->CurrentVolID(id[1]); // get the phi number inside the layer
203 :
204 0 : Int_t tower = (id[0]-1) % geom->GetNZ() + 1 + (id[1] - 1) * geom->GetNZ() ;
205 0 : Int_t layer = static_cast<Int_t>((id[0]-1)/(geom->GetNZ())) + 1 ;
206 : Int_t absid = tower;
207 0 : Int_t nlayers = geom->GetNECLayers();
208 0 : if ((layer > nlayers)||(layer<1))
209 0 : Fatal("StepManager", "Wrong calculation of layer number: layer = %d > %d\n", layer, nlayers) ;
210 :
211 : Float_t lightYield = depositedEnergy ;
212 : // Apply Birk's law (copied from G3BIRK)
213 :
214 0 : if (TVirtualMC::GetMC()->TrackCharge()!=0) { // Check
215 : Float_t birkC1Mod = 0;
216 0 : if (fBirkC0==1){ // Apply correction for higher charge states
217 0 : if (TMath::Abs(TVirtualMC::GetMC()->TrackCharge())>=2)
218 0 : birkC1Mod=fBirkC1*7.2/12.6;
219 : else
220 0 : birkC1Mod=fBirkC1;
221 : }
222 : Float_t dedxcm=0.;
223 0 : if (TVirtualMC::GetMC()->TrackStep()>0)
224 0 : dedxcm=1000.*TVirtualMC::GetMC()->Edep()/TVirtualMC::GetMC()->TrackStep();
225 : else
226 : dedxcm=0;
227 0 : lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm);
228 0 : }
229 :
230 : // use sampling fraction to get original energy --HG
231 0 : xyzte[4] = lightYield * geom->GetSampling();
232 :
233 0 : if (gDebug == 2)
234 0 : printf("StepManager: id0 = %d, id1 = %d, absid = %d tower = %d layer = %d energy = %f\n", id[0], id[1], absid, tower, layer, xyzte[4]) ;
235 :
236 0 : AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid, xyzte, pmom);
237 0 : } // there is deposited energy
238 0 : }
239 0 : }
240 :
241 : //___________________________________________________________
242 : void AliEMCALv1::RemapTrackHitIDs(Int_t *map) {
243 : // remap track index numbers for primary and parent indices
244 : // (Called by AliStack::PurifyKine)
245 224 : if (Hits()==0)
246 : return;
247 112 : TIter hitIter(Hits());
248 : Int_t iHit=0;
249 1728 : while (AliEMCALHit *hit=dynamic_cast<AliEMCALHit*>(hitIter()) ) {
250 292 : if (map[hit->GetIparent()]==-99)
251 0 : cout << "Remapping, found -99 for parent id " << hit->GetIparent() << ", " << map[hit->GetIparent()] << ", iHit " << iHit << endl;
252 292 : hit->SetIparent(map[hit->GetIparent()]);
253 292 : if (map[hit->GetPrimary()]==-99)
254 0 : cout << "Remapping, found -99 for primary id " << hit->GetPrimary() << ", " << map[hit->GetPrimary()] << ", iHit " << iHit << endl;
255 292 : hit->SetPrimary(map[hit->GetPrimary()]);
256 292 : iHit++;
257 292 : }
258 224 : }
259 :
260 : //___________________________________________________________
261 : void AliEMCALv1::FinishPrimary() {
262 : // finish primary
263 224 : fCurPrimary=-1;
264 112 : fCurParent=-1;
265 112 : fCurTrack=-1;
266 112 : }
|
