Removing ^M at the end of lines
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALv1.cxx
1 /**************************************************************************
2
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4
5  *                                                                        *
6
7  * Author: The ALICE Off-line Project.                                    *
8
9  * Contributors are mentioned in the code where appropriate.              *
10
11  *                                                                        *
12
13  * Permission to use, copy, modify and distribute this software and its   *
14
15  * documentation strictly for non-commercial purposes is hereby granted   *
16
17  * without fee, provided that the above copyright notice appears in all   *
18
19  * copies and that both the copyright notice and this permission notice   *
20
21  * appear in the supporting documentation. The authors make no claims     *
22
23  * about the suitability of this software for any purpose. It is          *
24
25  * provided "as is" without express or implied warranty.                  *
26
27  **************************************************************************/
28
29
30
31 /* $Id$ */
32
33
34
35 //_________________________________________________________________________
36
37 // Implementation version v1 of EMCAL Manager class 
38
39 // An object of this class does not produce digits
40
41 // It is the one to use if you do want to produce outputs in TREEH 
42
43 //                  
44
45 //*-- Author: Sahal Yacoob (LBL /UCT)
46
47 //*--       : Jennifer Klay (LBL)
48
49
50
51 // This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis.
52
53 // This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL.
54
55
56
57 // 15/02/2002 .... Yves Schutz
58
59 //  1. fSamplingFraction and fLayerToPreshowerRatio have been removed
60
61 //  2. Timing signal is collected and added to hit
62
63
64
65 // --- ROOT system ---
66
67 #include "TPGON.h"
68
69 #include "TTUBS.h"
70
71 #include "TNode.h"
72
73 #include "TRandom.h"
74
75 #include "TTree.h"
76
77 #include "TGeometry.h"
78
79 #include "TParticle.h"
80
81
82
83 // --- Standard library ---
84
85
86
87 #include <stdio.h>
88
89 #include <string.h>
90
91 #include <stdlib.h>
92
93 #include <strstream.h>
94
95 #include <iostream.h>
96
97 #include <math.h>
98
99 // --- AliRoot header files ---
100
101
102
103 #include "AliEMCALv1.h"
104
105 #include "AliEMCALHit.h"
106
107 #include "AliEMCALGeometry.h"
108
109 #include "AliConst.h"
110
111 #include "AliRun.h"
112
113 #include "AliMC.h"
114
115
116
117 ClassImp(AliEMCALv1)
118
119
120
121
122
123 //______________________________________________________________________
124
125 AliEMCALv1::AliEMCALv1():AliEMCALv0(){
126
127   // ctor
128
129     fLightYieldMean = 0 ; 
130
131     fIntrinsicAPDEfficiency = fLightFactor = fLightYieldAttenuation =   fAPDFactor = fAPDGain = fRecalibrationFactor = fAPDFactor = 0. ; 
132
133 }
134
135
136
137 //______________________________________________________________________
138
139 AliEMCALv1::AliEMCALv1(const char *name, const char *title):
140
141     AliEMCALv0(name,title){
142
143     // Standard Creator.
144
145
146
147     fHits= new TClonesArray("AliEMCALHit",1000);
148
149     gAlice->AddHitList(fHits);
150
151
152
153     fNhits = 0;
154
155     fIshunt     =  2; // All hits are associated with particles entering the calorimeter
156
157
158
159     //Photoelectron statistics:
160
161     // The light yield is a poissonian distribution of the number of
162
163     // photons created in a plastic layer, calculated using following formula
164
165     // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
166
167     //              exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit)
168
169     // LightYieldMean is parameter calculated to be over 100000 photons per GeV (a guess)
170
171     // APDEfficiency is 0.02655
172
173     // fLightYieldAttenuation is 0.0045 a guess
174
175     // TO BE FIXED
176
177     //***** Need a method in geometry to retrieve the fiber length corresponding to each layer
178
179     //***** See the step manager for the light attenuation calculation 
180
181     // The number of electrons created in the APD is
182
183     // NumberOfElectrons = APDGain * LightYield
184
185     // The APD Gain is 300
186
187     
188
189     fLightYieldMean         = 10000000.;  // This is a guess
190
191     fIntrinsicAPDEfficiency = 0.02655 ;
192
193     fLightFactor            = fLightYieldMean * fIntrinsicAPDEfficiency ; 
194
195     fLightYieldAttenuation  = 0.0045 ; // an other guess 
196
197     fAPDGain                = 300. ;
198
199     fRecalibrationFactor    = 13.418/ fLightYieldMean ;
200
201     fAPDFactor              = (fRecalibrationFactor/100.) * fAPDGain ; 
202
203
204
205 }
206
207
208
209 //______________________________________________________________________
210
211 AliEMCALv1::~AliEMCALv1(){
212
213     // dtor
214
215
216
217     if ( fHits) {
218
219         fHits->Delete();
220
221         delete fHits;
222
223         fHits = 0;
224
225     }
226
227 }
228
229 //______________________________________________________________________
230
231 void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, 
232
233                         Int_t id, Float_t * hits,Float_t * p){
234
235     // Add a hit to the hit list.
236
237     // An EMCAL hit is the sum of all hits in a single segment 
238
239     //   originating from the same enterring particle 
240
241     Int_t hitCounter;
242
243     
244
245     AliEMCALHit *newHit;
246
247     AliEMCALHit *curHit;
248
249     Bool_t deja = kFALSE;
250
251
252
253     newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
254
255     for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) {
256
257         curHit = (AliEMCALHit*) (*fHits)[hitCounter];
258
259         // We add hits with the same tracknumber, while GEANT treats
260
261         // primaries succesively
262
263         if(curHit->GetPrimary() != primary) break;
264
265         if( *curHit == *newHit ) {
266
267             *curHit = *curHit + *newHit;
268
269             deja = kTRUE;
270
271         } // end if
272
273     } // end for hitCounter
274
275
276
277     if ( !deja ) {
278
279         new((*fHits)[fNhits]) AliEMCALHit(*newHit);
280
281         fNhits++;
282
283     } // end if
284
285
286
287     delete newHit;
288
289 }
290
291 //______________________________________________________________________
292
293 void AliEMCALv1::StepManager(void){
294
295   // Accumulates hits as long as the track stays in a single
296
297   // crystal or PPSD gas Cell
298
299
300
301   Int_t          id[2];           // (layer, phi, Eta) indices
302
303   Int_t          absid;
304
305   // position wrt MRS and energy deposited
306
307   Float_t        xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
308
309   Float_t        pmom[4]={0.,0.,0.,0.};
310
311   TLorentzVector pos; // Lorentz vector of the track current position.
312
313   TLorentzVector mom; // Lorentz vector of the track current momentum.
314
315   Int_t tracknumber =  gAlice->CurrentTrack();
316
317   Int_t primary = 0;
318
319   static Int_t iparent = 0;
320
321   static Float_t ienergy = 0;
322
323   Int_t copy = 0;
324
325   
326
327   AliEMCALGeometry * geom = GetGeometry() ; 
328
329
330
331   if(gMC->IsTrackEntering() && (strcmp(gMC->CurrentVolName(),"XALU") == 0)){ // This Particle in enterring the Calorimeter
332
333     gMC->TrackPosition(pos) ;
334
335     xyzte[0] = pos[0] ;
336
337     xyzte[1] = pos[1] ;
338
339     xyzte[2] = pos[2] ;
340
341     if ( (xyzte[0]*xyzte[0] + xyzte[1]*xyzte[1])  
342
343          <  (geom->GetEnvelop(0)+geom->GetGap2Active()+1.5 )*(geom->GetEnvelop(0)+geom->GetGap2Active()+1.5 ) ) {
344
345       iparent = tracknumber;
346
347       gMC->TrackMomentum(mom);
348
349       ienergy = mom[3]; 
350
351       TParticle * part = 0 ;
352
353       Int_t parent = iparent ;
354
355       while ( parent != -1 ) { // <------------- flags this particle to be kept and
356
357         //all the ancestors of this particle
358
359         part = gAlice->Particle(parent) ;
360
361         part->SetBit(kKeepBit);
362
363         parent = part->GetFirstMother() ;
364
365       }
366
367     }
368
369   }
370
371   if(gMC->CurrentVolID(copy) == gMC->VolId("XPHI") ) { // We are in a Scintillator Layer 
372
373     
374
375     Float_t depositedEnergy ; 
376
377     
378
379     if( (depositedEnergy = gMC->Edep()) > 0.){// Track is inside a scintillator and deposits some energy
380
381       
382
383       gMC->TrackPosition(pos);
384
385       xyzte[0] = pos[0];
386
387       xyzte[1] = pos[1];
388
389       xyzte[2] = pos[2];
390
391       xyzte[3] = gMC->TrackTime() ;       
392
393       
394
395       gMC->TrackMomentum(mom);
396
397       pmom[0] = mom[0];
398
399       pmom[1] = mom[1];
400
401       pmom[2] = mom[2];
402
403       pmom[3] = mom[3];
404
405       
406
407       gMC->CurrentVolOffID(1, id[0]); // get the POLY copy number;
408
409       gMC->CurrentVolID(id[1]); // get the phi number inside the layer
410
411       absid = (id[0]-1)*(geom->GetNPhi()) + id[1];
412
413       
414
415       //Calculates the light yield, the number of photons produced in the
416
417       //plastic layer 
418
419       // Here we need to know the fiber lebgth to calculate the attenuation
420
421      
422
423       Float_t lengthOfFiber = 0. ;// should be retrieved from the geometry
424
425
426
427       Float_t lightYield = gRandom->Poisson(fLightFactor * depositedEnergy *
428
429                                             exp(-fLightYieldAttenuation * lengthOfFiber)) ;
430
431       xyzte[4] = fAPDFactor * lightYield  ;
432
433    
434
435       primary = gAlice->GetPrimary(tracknumber);
436
437       AddHit(fIshunt, primary,tracknumber, iparent, ienergy, absid, xyzte, pmom);
438
439     } // there is deposited energy
440
441   }
442
443 }
444