bag in AliPHOS::StepMamager fixed, AliPHOSCPVModule.* removed
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv3.cxx
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 v0 of PHOS Manager class 
20 // Layout EMC + PPSD has name GPS2  
21 // The main goal of this version of AliPHOS is to calculte the 
22 //  induced charged in the PIN diode, taking into account light
23 //  tracking in the PbWO4 crystal, induced signal in the 
24 //  PIN due to MIPS particle and electronic noise.
25 // This is done in the StepManager 
26 //                  
27 //*-- Author:  Odd Harald Oddland & Gines Martinez (SUBATECH)
28
29
30 // --- ROOT system ---
31 #include "TRandom.h"
32
33 // --- Standard library ---
34
35 #include <stdio.h>
36 #include <string.h>
37 #include <stdlib.h>
38 #include <strstream.h>
39
40 // --- AliRoot header files ---
41
42 #include "AliPHOSv3.h"
43 #include "AliPHOSHit.h"
44 #include "AliPHOSDigit.h"
45 #include "AliRun.h"
46 #include "AliConst.h"
47 #include "AliMC.h"
48
49 ClassImp(AliPHOSv3)
50
51 //____________________________________________________________________________
52   AliPHOSv3::AliPHOSv3(const char *name, const char *title):
53 AliPHOSv1(name,title)
54 {
55   // ctor 
56
57   // Number of electrons created in the PIN due to light collected in the PbWo4 crystal is calculated using 
58   // following formula
59   // NumberOfElectrons = EnergyLost * LightYield * PINEfficiency * 
60   //                     exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit) *
61   //                     RecalibrationFactor ;
62   // LightYield is obtained as a Poissonian distribution with a mean at 700000 photons per GeV fromValery Antonenko
63   // PINEfficiency is 0.1875 from Odd Harald Odland work
64   // k_0 is 0.0045 from Valery Antonenko 
65
66
67   fLightYieldMean = 700000. ;
68   fIntrinsicPINEfficiency = 0.1875 ;
69   fLightYieldAttenuation = 0.0045 ;
70   fRecalibrationFactor = 6.2 / fLightYieldMean ;
71   fElectronsPerGeV = 2.77e+8 ; 
72 }
73
74 //____________________________________________________________________________
75 AliPHOSv3::AliPHOSv3(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
76   AliPHOSv1(Reconstructioner,name,title)
77 {
78   // ctor 
79
80   // Number of electrons created in the PIN due to light collected in the PbWo4 crystal is calculated using 
81   // following formula
82   // NumberOfElectrons = EnergyLost * LightYield * PINEfficiency * 
83   //                     exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit) *
84   //                     RecalibrationFactor ;
85   // LightYield is obtained as a Poissonian distribution with a mean at 700000 photons per GeV fromValery Antonenko
86   // PINEfficiency is 0.1875 from Odd Harald Odland work
87   // k_0 is 0.0045 from Valery Antonenko 
88
89   fLightYieldMean = 700000.;
90   fIntrinsicPINEfficiency = 0.1875 ;
91   fLightYieldAttenuation = 0.0045 ;
92   fRecalibrationFactor = 6.2 / fLightYieldMean ;
93   fElectronsPerGeV = 2.77e+8 ;
94 }
95 //____________________________________________________________________________
96
97 void AliPHOSv3::StepManager(void)
98 {
99   // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
100
101 //    if (gMC->IsTrackEntering())
102 //      cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
103 //    if (gMC->IsTrackExiting())
104 //      cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
105
106   Int_t          relid[4] ;        // (box, layer, row, column) indices
107   Int_t          absid    ;        // absolute cell ID number
108   Float_t        xyze[4]={0,0,0,0}  ; // position wrt MRS and energy deposited
109   TLorentzVector pos      ;        // Lorentz vector of the track current position
110   TLorentzVector pmom     ;        //momentum of the particle initiated hit
111   Float_t        xyd[2]={0,0}   ;  //local posiiton of the entering
112   Bool_t         entered = kFALSE    ;  
113   Int_t          copy     ;
114
115   Int_t tracknumber =  gAlice->CurrentTrack() ; 
116   Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
117   TString name      =  fGeom->GetName() ; 
118   Int_t trackpid    =  0  ; 
119
120   if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, 
121                                 // but may be without energy deposition
122
123     // Current position of the hit in the local ref. system
124       gMC -> TrackPosition(pos);
125       Float_t xyzm[3], xyzd[3] ;
126       Int_t i;
127       for (i=0; i<3; i++) xyzm[i] = pos[i];
128       gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
129       xyd[0]  = xyzd[0];
130       xyd[1]  =-xyzd[2];
131       
132       // Current momentum of the hit's track in the local ref. system
133       gMC -> TrackMomentum(pmom);
134       Float_t pm[3], pd[3];
135       for (i=0; i<3; i++) pm[i]   = pmom[i];
136       gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
137       pmom[0] = pd[0];
138       pmom[1] =-pd[1];
139       pmom[2] =-pd[2];
140
141       trackpid = gMC->TrackPid();
142       entered = kTRUE ;      // Mark to create hit even withou energy deposition
143
144   }
145
146
147   if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
148
149     if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell 
150     {
151       gMC->TrackPosition(pos) ;
152       xyze[0] = pos[0] ;
153       xyze[1] = pos[1] ;
154       xyze[2] = pos[2] ;
155       xyze[3] = gMC->Edep() ; 
156
157       if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering  PPSD
158         gMC->CurrentVolOffID(5, relid[0]) ;  // get the PHOS Module number
159         if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
160           relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
161         }
162         gMC->CurrentVolOffID(3, relid[1]) ;  // get the Micromegas Module number 
163       // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
164       //   > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
165         gMC->CurrentVolOffID(1, relid[2]) ;  // get the row number of the cell
166         gMC->CurrentVolID(relid[3]) ;        // get the column number 
167
168         // get the absolute Id number
169
170         fGeom->RelToAbsNumbering(relid, absid) ; 
171
172         // add current hit to the hit list      
173           AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
174
175
176       } // there is deposited energy 
177     } // We are inside the gas of the CPV  
178   } // GPS2 configuration
179
180   if ( name == "IHEP" || name == "MIXT" ) {       // ======> CPV is a IHEP's one
181
182     // Yuri Kharlov, 28 September 2000
183
184     if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
185         entered &&
186         gMC->TrackCharge() != 0) {      
187       
188       // Digitize the current CPV hit:
189
190       // 1. find pad response and
191       
192       Int_t moduleNumber;
193       gMC->CurrentVolOffID(3,moduleNumber);
194       moduleNumber--;
195
196
197       TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
198       CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
199       
200       Float_t xmean = 0;
201       Float_t zmean = 0;
202       Float_t qsum  = 0;
203       Int_t   idigit,ndigits;
204
205       // 2. go through the current digit list and sum digits in pads
206
207       ndigits = cpvDigits->GetEntriesFast();
208       for (idigit=0; idigit<ndigits-1; idigit++) {
209         AliPHOSCPVDigit  *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
210         Float_t x1 = cpvDigit1->GetXpad() ;
211         Float_t z1 = cpvDigit1->GetYpad() ;
212         for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
213           AliPHOSCPVDigit  *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
214           Float_t x2 = cpvDigit2->GetXpad() ;
215           Float_t z2 = cpvDigit2->GetYpad() ;
216           if (x1==x2 && z1==z2) {
217             Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
218             cpvDigit2->SetQpad(qsum) ;
219             cpvDigits->RemoveAt(idigit) ;
220           }
221         }
222       }
223       cpvDigits->Compress() ;
224
225       // 3. add digits to temporary hit list fTmpHits
226
227       ndigits = cpvDigits->GetEntriesFast();
228       for (idigit=0; idigit<ndigits; idigit++) {
229         AliPHOSCPVDigit  *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
230         relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
231         relid[1] =-1 ;                                            // means CPV
232         relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
233         relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
234         
235         // get the absolute Id number
236         fGeom->RelToAbsNumbering(relid, absid) ; 
237
238         // add current digit to the temporary hit list
239         xyze[0] = 0. ;
240         xyze[1] = 0. ;
241         xyze[2] = 0. ;
242         xyze[3] = cpvDigit->GetQpad() ;                           // amplitude in a pad
243         primary = -1;                                             // No need in primary for CPV
244         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
245
246         if (cpvDigit->GetQpad() > 0.02) {
247           xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
248           zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
249           qsum  += cpvDigit->GetQpad();
250         }
251       }
252       delete cpvDigits;
253     }
254   } // end of IHEP configuration
255   
256
257   if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
258     gMC->TrackPosition(pos) ;
259     xyze[0] = pos[0] ;
260     xyze[1] = pos[1] ;
261     xyze[2] = pos[2] ;
262     xyze[3] = gMC->Edep() ;
263
264   
265     if ( (xyze[3] != 0) || entered ) {  // Track is inside the crystal and deposits some energy or just entered 
266
267       gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
268
269       if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
270         relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();      
271
272       relid[1] = 0   ;                    // means PBW04
273       gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
274       gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
275       
276       // get the absolute Id number
277       fGeom->RelToAbsNumbering(relid, absid) ; 
278
279       // add current hit to the hit list
280         AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
281
282
283     } // there is deposited energy
284   } // we are inside a PHOS Xtal
285
286   if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode 
287     {
288       gMC->TrackPosition(pos) ;
289       xyze[0] = pos[0] ;
290       xyze[1] = pos[1] ;
291       xyze[2] = pos[2] ;
292       Float_t lostenergy = gMC->Edep() ;
293       xyze[3] = gMC->Edep() ;
294       
295       if ( xyze[3] != 0 ) {
296         gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ;
297         relid[1] = 0   ;                    // means PW04 and PIN
298         gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module
299         gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module
300         
301         // get the absolute Id number
302         
303         Int_t absid ; 
304         fGeom->RelToAbsNumbering(relid,absid) ;
305         
306         // calculating number of electrons in the PIN diode asociated to this hit
307           Float_t nElectrons = lostenergy * fElectronsPerGeV ;
308           xyze[3] = nElectrons * fRecalibrationFactor ;
309           
310           // add current hit to the hit list
311           AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid,pmom,xyd);
312           //printf("PIN volume is  %d, %d, %d, %d \n",relid[0],relid[1],relid[2],relid[3]);
313           //printf("Lost energy in the PIN is %f \n",lostenergy) ;
314       } // there is deposited energy
315     } // we are inside a PHOS XtalPHOS PIN diode
316 }