6155bacd7490d72f4353fc9c58aba01a452d9883
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.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 v1 of PHOS Manager class 
20 //---
21 // Layout EMC + PPSD has name GPS2:
22 // Produces cumulated hits
23 //---
24 // Layout EMC + CPV  has name IHEP:
25 // Produces hits for CPV, cumulated hits
26 //---
27 // Layout EMC + CPV + PPSD has name GPS:
28 // Produces hits for CPV, cumulated hits
29 //---
30 //*-- Author: Yves Schutz (SUBATECH)
31
32
33 // --- ROOT system ---
34
35 #include "TBRIK.h"
36 #include "TNode.h"
37 #include "TRandom.h"
38 #include "TTree.h"
39 #include "TParticle.h"
40
41 // --- Standard library ---
42
43 #include <string.h>
44 #include <stdlib.h>
45 #include <strstream.h>
46
47 // --- AliRoot header files ---
48
49 #include "AliPHOSv1.h"
50 #include "AliPHOSHit.h"
51 #include "AliPHOSCPVDigit.h"
52 #include "AliRun.h"
53 #include "AliConst.h"
54 #include "AliMC.h"
55 #include "AliPHOSGeometry.h"
56 #include "AliPHOSQAIntCheckable.h"
57 #include "AliPHOSQAFloatCheckable.h"
58 #include "AliPHOSQAMeanChecker.h"
59
60 ClassImp(AliPHOSv1)
61
62 //____________________________________________________________________________
63 AliPHOSv1::AliPHOSv1():
64 AliPHOSv0()
65 {
66   // default ctor: initialze data memebers
67   fQAHitsMul  = 0 ;
68   fQAHitsMulB = 0 ; 
69   fQATotEner  = 0 ; 
70   fQATotEnerB = 0 ; 
71
72   fLightYieldMean         = 0. ;         
73   fIntrinsicPINEfficiency = 0. ; 
74   fLightYieldAttenuation  = 0. ;  
75   fRecalibrationFactor    = 0. ;    
76   fElectronsPerGeV        = 0. ;
77   fAPDGain                = 0. ;  
78   fLightFactor            = 0. ; 
79   fAPDFactor              = 0. ; 
80
81 }
82
83 //____________________________________________________________________________
84 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
85  AliPHOSv0(name,title) 
86 {
87   //
88   // We store hits :
89   //   - fHits (the "normal" one), which retains the hits associated with
90   //     the current primary particle being tracked
91   //     (this array is reset after each primary has been tracked).
92   //
93
94
95
96   // We do not want to save in TreeH the raw hits
97   // But save the cumulated hits instead (need to create the branch myself)
98   // It is put in the Digit Tree because the TreeH is filled after each primary
99   // and the TreeD at the end of the event (branch is set in FinishEvent() ). 
100   
101   fHits= new TClonesArray("AliPHOSHit",1000) ;
102   gAlice->AddHitList(fHits) ; 
103
104   fNhits = 0 ;
105
106   fIshunt     =  2 ; // All hits are associated with primary particles
107
108   //Photoelectron statistics:
109   // The light yield is a poissonian distribution of the number of
110   // photons created in the PbWo4 crystal, calculated using following formula
111   // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
112   //              exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
113   // LightYieldMean is parameter calculated to be over 47000 photons per GeV
114   // APDEfficiency is 0.02655
115   // k_0 is 0.0045 from Valery Antonenko
116   // The number of electrons created in the APD is
117   // NumberOfElectrons = APDGain * LightYield
118   // The APD Gain is 300
119   fLightYieldMean = 47000;
120   fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
121   fLightYieldAttenuation  = 0.0045 ; 
122   fRecalibrationFactor    = 13.418/ fLightYieldMean ;
123   fElectronsPerGeV        = 2.77e+8 ;
124   fAPDGain                = 300. ;
125   fLightFactor            = fLightYieldMean * fIntrinsicPINEfficiency ; 
126   fAPDFactor              = (fRecalibrationFactor/100.) * fAPDGain ; 
127
128
129   Int_t nb   = GetGeometry()->GetNModules() ; 
130   
131   // create checkables 
132   fQAHitsMul   = new AliPHOSQAIntCheckable("HitsM") ; 
133   fQATotEner   = new AliPHOSQAFloatCheckable("TotEn") ; 
134   fQAHitsMulB  = new TClonesArray("AliPHOSQAIntCheckable",nb) ;
135   fQAHitsMulB->SetOwner() ; 
136   fQATotEnerB  = new TClonesArray("AliPHOSQAFloatCheckable", nb); 
137   fQATotEnerB->SetOwner() ; 
138   char tempo[20]  ; 
139   Int_t i ; 
140   for ( i = 0 ; i < nb ; i++ ) {
141     sprintf(tempo, "HitsMB%d", i+1) ; 
142     new( (*fQAHitsMulB)[i]) AliPHOSQAIntCheckable(tempo) ; 
143     sprintf(tempo, "TotEnB%d", i+1) ; 
144     new( (*fQATotEnerB)[i] ) AliPHOSQAFloatCheckable(tempo) ;
145   }
146
147   AliPHOSQAMeanChecker * hmc  = new AliPHOSQAMeanChecker("HitsMul", 100. ,25.) ; 
148   AliPHOSQAMeanChecker * emc  = new AliPHOSQAMeanChecker("TotEner", 10. ,5.) ; 
149   AliPHOSQAMeanChecker * bhmc = new AliPHOSQAMeanChecker("HitsMulB", 100. ,5.) ; 
150   AliPHOSQAMeanChecker * bemc = new AliPHOSQAMeanChecker("TotEnerB", 2. ,.5) ; 
151
152   // associate checkables and checkers 
153   fQAHitsMul->AddChecker(hmc) ; 
154   fQATotEner->AddChecker(emc) ; 
155   for ( i = 0 ; i < nb ; i++ ) {
156     (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]))->AddChecker(bhmc) ;
157     (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]))->AddChecker(bemc) ; 
158   }
159
160 }
161
162 //____________________________________________________________________________
163 AliPHOSv1::~AliPHOSv1()
164 {
165   // dtor
166
167   if ( fHits) {
168     fHits->Delete() ; 
169     delete fHits ;
170     fHits = 0 ; 
171   }
172   
173   delete fTreeQA ; 
174
175   if ( fQAHitsMulB ) {
176     fQAHitsMulB->Delete() ;
177     delete fQAHitsMulB ; 
178   }
179
180   if ( fQATotEnerB ) {
181     fQATotEnerB->Delete() ;
182     delete fQATotEnerB ; 
183   }
184  
185 }
186
187 //____________________________________________________________________________
188 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
189 {
190   // Add a hit to the hit list.
191   // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
192
193   Int_t hitCounter ;
194   AliPHOSHit *newHit ;
195   AliPHOSHit *curHit ;
196   Bool_t deja = kFALSE ;
197   AliPHOSGeometry * geom = GetGeometry() ; 
198
199   newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
200
201   for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
202     curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
203     if(curHit->GetPrimary() != primary) break ; 
204            // We add hits with the same primary, while GEANT treats primaries succesively 
205     if( *curHit == *newHit ) {
206       *curHit + *newHit ;
207       deja = kTRUE ;
208     }
209   }
210          
211   if ( !deja ) {
212     new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
213     // get the block Id number
214     Int_t relid[4] ;
215     geom->AbsToRelNumbering(Id, relid) ;
216     // and fill the relevant QA checkable (only if in PbW04)
217     if ( relid[1] == 0 ) {
218       fQAHitsMul->Update(1) ; 
219       (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[relid[0]-1]))->Update(1) ;
220     } 
221     fNhits++ ;
222   }
223
224   delete newHit;
225 }
226
227 //____________________________________________________________________________
228 void AliPHOSv1::FinishPrimary() 
229 {
230   // called at the end of each track (primary) by AliRun
231   // hits are reset for each new track
232   // accumulate the total hit-multiplicity
233 //   if ( fQAHitsMul ) 
234 //     fQAHitsMul->Update( fHits->GetEntriesFast() ) ; 
235
236 }
237
238 //____________________________________________________________________________
239 void AliPHOSv1::FinishEvent() 
240 {
241   // called at the end of each event by AliRun
242   // accumulate the hit-multiplicity and total energy per block 
243   // if the values have been updated check it
244
245   if ( fQATotEner ) { 
246     if ( fQATotEner->HasChanged() ) {
247       fQATotEner->CheckMe() ; 
248       fQATotEner->Reset() ; 
249     }
250   }
251   
252   Int_t i ; 
253   if ( fQAHitsMulB && fQATotEnerB ) {
254     for (i = 0 ; i < GetGeometry()->GetNModules() ; i++) {
255       AliPHOSQAIntCheckable * ci = static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]) ;  
256       AliPHOSQAFloatCheckable* cf = static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]) ; 
257       if ( ci->HasChanged() ) { 
258         ci->CheckMe() ;  
259         ci->Reset() ;
260       } 
261       if ( cf->HasChanged() ) { 
262         cf->CheckMe() ; 
263         cf->Reset() ;
264       }
265     } 
266   }
267   
268   // check the total multiplicity 
269   
270   if ( fQAHitsMul ) {
271     if ( fQAHitsMul->HasChanged() ) { 
272       fQAHitsMul->CheckMe() ; 
273       fQAHitsMul->Reset() ; 
274     }
275   } 
276 }
277 //____________________________________________________________________________
278 void AliPHOSv1::StepManager(void)
279 {
280    // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
281
282   Int_t          relid[4] ;           // (box, layer, row, column) indices
283   Int_t          absid    ;           // absolute cell ID number
284   Float_t        xyzte[5]={-1000.,-1000.,-1000.,0.,0.}  ; // position wrt MRS, time and energy deposited
285   TLorentzVector pos      ;           // Lorentz vector of the track current position
286   Int_t          copy     ;
287
288   Int_t tracknumber =  gAlice->CurrentTrack() ; 
289   Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
290   TString name      =  GetGeometry()->GetName() ; 
291
292   Int_t moduleNumber ;
293   
294   if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
295       (gMC->IsTrackEntering() ) &&
296       gMC->TrackCharge() != 0) {      
297     
298     gMC -> TrackPosition(pos);
299     
300     Float_t xyzm[3], xyzd[3] ;
301     Int_t i;
302     for (i=0; i<3; i++) xyzm[i] = pos[i];
303     gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
304     
305     Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
306     xyd[0]  = xyzd[0];
307     xyd[1]  =-xyzd[2];
308     xyd[2]  =-xyzd[1];
309     
310     // Current momentum of the hit's track in the local ref. system
311     TLorentzVector pmom     ;        //momentum of the particle initiated hit
312     gMC -> TrackMomentum(pmom);
313     Float_t pm[3], pd[3];
314     for (i=0; i<3; i++)  
315       pm[i]   = pmom[i];
316     
317     gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
318     pmom[0] = pd[0];
319     pmom[1] =-pd[1];
320     pmom[2] =-pd[2];
321
322     // Digitize the current CPV hit:
323     
324     // 1. find pad response and    
325     gMC->CurrentVolOffID(3,moduleNumber);
326     moduleNumber--;
327     
328     TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
329     CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
330       
331     Float_t xmean = 0;
332     Float_t zmean = 0;
333     Float_t qsum  = 0;
334     Int_t   idigit,ndigits;
335     
336     // 2. go through the current digit list and sum digits in pads
337     
338     ndigits = cpvDigits->GetEntriesFast();
339     for (idigit=0; idigit<ndigits-1; idigit++) {
340       AliPHOSCPVDigit  *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
341       Float_t x1 = cpvDigit1->GetXpad() ;
342       Float_t z1 = cpvDigit1->GetYpad() ;
343       for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
344         AliPHOSCPVDigit  *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
345         Float_t x2 = cpvDigit2->GetXpad() ;
346         Float_t z2 = cpvDigit2->GetYpad() ;
347         if (x1==x2 && z1==z2) {
348           Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
349           cpvDigit2->SetQpad(qsum) ;
350           cpvDigits->RemoveAt(idigit) ;
351         }
352       }
353     }
354     cpvDigits->Compress() ;
355     
356     // 3. add digits to temporary hit list fTmpHits
357     
358     ndigits = cpvDigits->GetEntriesFast();
359     for (idigit=0; idigit<ndigits; idigit++) {
360       AliPHOSCPVDigit  *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
361       relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
362       relid[1] =-1 ;                                            // means CPV
363       relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
364       relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
365       
366       // get the absolute Id number
367       GetGeometry()->RelToAbsNumbering(relid, absid) ; 
368       
369       // add current digit to the temporary hit list
370
371       xyzte[3] = gMC->TrackTime() ;
372       xyzte[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
373       primary = -1;                                             // No need in primary for CPV
374       AddHit(fIshunt, primary, tracknumber, absid, xyzte);
375       
376       if (cpvDigit->GetQpad() > 0.02) {
377         xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
378         zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
379         qsum  += cpvDigit->GetQpad();
380       }
381     }
382     if (cpvDigits) {
383       cpvDigits->Delete();
384       delete cpvDigits;
385       cpvDigits=0;
386     }
387   }
388
389  
390   
391   if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
392
393     gMC->TrackPosition(pos) ;
394     xyzte[0] = pos[0] ;
395     xyzte[1] = pos[1] ;
396     xyzte[2] = pos[2] ;
397
398     Float_t global[3], local[3] ;
399     global[0] = pos[0] ;
400     global[1] = pos[1] ;
401     global[2] = pos[2] ;
402     Float_t lostenergy = gMC->Edep(); 
403     
404     //Put in the TreeK particle entering PHOS and all its parents
405     if ( gMC->IsTrackEntering() ){
406       Float_t xyzd[3] ;
407       gMC -> Gmtod (xyzte, xyzd, 1);    // transform coordinate from master to daughter system    
408       if (xyzd[1] >  GetGeometry()->GetCrystalSize(1)/2-0.002 ||
409           xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2+0.002) {
410         TParticle * part = 0 ; 
411         Int_t parent = gAlice->CurrentTrack() ; 
412         while ( parent != -1 ) {
413           part = gAlice->Particle(parent) ; 
414           part->SetBit(kKeepBit);
415           parent = part->GetFirstMother() ; 
416         }
417       }
418     }
419     if ( lostenergy != 0 ) {  // Track is inside the crystal and deposits some energy 
420       xyzte[3] = gMC->TrackTime() ;     
421       
422       gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
423       
424       Int_t strip ;
425       gMC->CurrentVolOffID(3, strip);
426       Int_t cell ;
427       gMC->CurrentVolOffID(2, cell);
428       
429       Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
430       Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
431       
432       absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + 
433         row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
434       
435       gMC->Gmtod(global, local, 1) ;
436       
437       //Calculates the light yield, the number of photons produced in the
438       //crystal 
439       Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
440                                             exp(-fLightYieldAttenuation *
441                                                 (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
442                                             ) ;
443
444       //Calculates de energy deposited in the crystal  
445       xyzte[4] = fAPDFactor * lightYield  ;
446       
447       // add current hit to the hit list
448       //cout << "AliPHOSv1::StepManager " << primary << " " << tracknumber << endl ; 
449       AddHit(fIshunt, primary,tracknumber, absid, xyzte);
450       
451       // fill the relevant QA Checkables
452       fQATotEner->Update( xyzte[4] ) ;                                             // total energy in PHOS
453       (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[moduleNumber-1]))->Update( xyzte[4] ) ; // energy in this block  
454       
455     } // there is deposited energy
456   } // we are inside a PHOS Xtal
457   
458 }
459
460 //____________________________________________________________________________
461 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
462 {
463   // ------------------------------------------------------------------------
464   // Digitize one CPV hit:
465   // On input take exact 4-momentum p and position zxhit of the hit,
466   // find the pad response around this hit and
467   // put the amplitudes in the pads into array digits
468   //
469   // Author: Yuri Kharlov (after Serguei Sadovsky)
470   // 2 October 2000
471   // ------------------------------------------------------------------------
472
473   const Float_t kCelWr  = GetGeometry()->GetPadSizePhi()/2;  // Distance between wires (2 wires above 1 pad)
474   const Float_t kDetR   = 0.1;     // Relative energy fluctuation in track for 100 e-
475   const Float_t kdEdx   = 4.0;     // Average energy loss in CPV;
476   const Int_t   kNgamz  = 5;       // Ionization size in Z
477   const Int_t   kNgamx  = 9;       // Ionization size in Phi
478   const Float_t kNoise = 0.03;    // charge noise in one pad
479
480   Float_t rnor1,rnor2;
481
482   // Just a reminder on axes notation in the CPV module:
483   // axis Z goes along the beam
484   // axis X goes across the beam in the module plane
485   // axis Y is a normal to the module plane showing from the IP
486
487   Float_t hitX  = zxhit[0];
488   Float_t hitZ  =-zxhit[1];
489   Float_t pX    = p.Px();
490   Float_t pZ    =-p.Pz();
491   Float_t pNorm = p.Py();
492   Float_t eloss = kdEdx;
493
494 //    cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
495
496   Float_t dZY   = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
497   Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
498   gRandom->Rannor(rnor1,rnor2);
499   eloss *= (1 + kDetR*rnor1) *
500            TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
501   Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
502   Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
503   Float_t zhit2 = zhit1 + dZY;
504   Float_t xhit2 = xhit1 + dXY;
505
506   Int_t   iwht1 = (Int_t) (xhit1 / kCelWr);           // wire (x) coordinate "in"
507   Int_t   iwht2 = (Int_t) (xhit2 / kCelWr);           // wire (x) coordinate "out"
508
509   Int_t   nIter;
510   Float_t zxe[3][5];
511   if (iwht1==iwht2) {                      // incline 1-wire hit
512     nIter = 2;
513     zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
514     zxe[1][0] = (iwht1 + 0.5) * kCelWr;
515     zxe[2][0] =  eloss/2;
516     zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
517     zxe[1][1] = (iwht1 + 0.5) * kCelWr;
518     zxe[2][1] =  eloss/2;
519   }
520   else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
521     nIter = 3;
522     Int_t iwht3 = (iwht1 + iwht2) / 2;
523     Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
524     Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
525     Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
526     Float_t xwr13 = (xwht1 + xwht3) / 2;   // center 13
527     Float_t xwr23 = (xwht2 + xwht3) / 2;   // center 23
528     Float_t dxw1  = xhit1 - xwr13;
529     Float_t dxw2  = xhit2 - xwr23;
530     Float_t egm1  = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
531     Float_t egm2  = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
532     Float_t egm3  =           kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
533     zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
534     zxe[1][0] =  xwht1;
535     zxe[2][0] =  eloss * egm1;
536     zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
537     zxe[1][1] =  xwht2;
538     zxe[2][1] =  eloss * egm2;
539     zxe[0][2] =  dXY*(xwht3-xwht1)/dXY + zhit1;
540     zxe[1][2] =  xwht3;
541     zxe[2][2] =  eloss * egm3;
542   }
543   else {                                   // incline 2-wire hit
544     nIter = 2;
545     Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
546     Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
547     Float_t xwr12 = (xwht1 + xwht2) / 2;
548     Float_t dxw1  = xhit1 - xwr12;
549     Float_t dxw2  = xhit2 - xwr12;
550     Float_t egm1  = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
551     Float_t egm2  = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
552     zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
553     zxe[1][0] =  xwht1;
554     zxe[2][0] =  eloss * egm1;
555     zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
556     zxe[1][1] =  xwht2;
557     zxe[2][1] =  eloss * egm2;
558   }
559
560   // Finite size of ionization region
561
562   Int_t nCellZ  = GetGeometry()->GetNumberOfCPVPadsZ();
563   Int_t nCellX  = GetGeometry()->GetNumberOfCPVPadsPhi();
564   Int_t nz3     = (kNgamz+1)/2;
565   Int_t nx3     = (kNgamx+1)/2;
566   cpvDigits->Expand(nIter*kNgamx*kNgamz);
567   TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
568
569   for (Int_t iter=0; iter<nIter; iter++) {
570
571     Float_t zhit = zxe[0][iter];
572     Float_t xhit = zxe[1][iter];
573     Float_t qhit = zxe[2][iter];
574     Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
575     Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
576     if ( zcell<=0      || xcell<=0 ||
577          zcell>=nCellZ || xcell>=nCellX) return;
578     Int_t izcell = (Int_t) zcell;
579     Int_t ixcell = (Int_t) xcell;
580     Float_t zc = zcell - izcell - 0.5;
581     Float_t xc = xcell - ixcell - 0.5;
582     for (Int_t iz=1; iz<=kNgamz; iz++) {
583       Int_t kzg = izcell + iz - nz3;
584       if (kzg<=0 || kzg>nCellZ) continue;
585       Float_t zg = (Float_t)(iz-nz3) - zc;
586       for (Int_t ix=1; ix<=kNgamx; ix++) {
587         Int_t kxg = ixcell + ix - nx3;
588         if (kxg<=0 || kxg>nCellX) continue;
589         Float_t xg = (Float_t)(ix-nx3) - xc;
590         
591         // Now calculate pad response
592         Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
593         qpad += kNoise*rnor2;
594         if (qpad<0) continue;
595         
596         // Fill the array with pad response ID and amplitude
597         new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
598       }
599     }
600   }
601 }
602
603 //____________________________________________________________________________
604 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
605   // ------------------------------------------------------------------------
606   // Calculate the amplitude in one CPV pad using the
607   // cumulative pad response function
608   // Author: Yuri Kharlov (after Serguei Sadovski)
609   // 3 October 2000
610   // ------------------------------------------------------------------------
611
612   Double_t dz = GetGeometry()->GetPadSizeZ()   / 2;
613   Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
614   Double_t z  = zhit * GetGeometry()->GetPadSizeZ();
615   Double_t x  = xhit * GetGeometry()->GetPadSizePhi();
616   Double_t amplitude = qhit *
617     (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
618      CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
619   return (Float_t)amplitude;
620 }
621
622 //____________________________________________________________________________
623 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
624   // ------------------------------------------------------------------------
625   // Cumulative pad response function
626   // It includes several terms from the CF decomposition in electrostatics
627   // Note: this cumulative function is wrong since omits some terms
628   //       but the cell amplitude obtained with it is correct because
629   //       these omitting terms cancel
630   // Author: Yuri Kharlov (after Serguei Sadovski)
631   // 3 October 2000
632   // ------------------------------------------------------------------------
633
634   const Double_t kA=1.0;
635   const Double_t kB=0.7;
636
637   Double_t r2       = x*x + y*y;
638   Double_t xy       = x*y;
639   Double_t cumulPRF = 0;
640   for (Int_t i=0; i<=4; i++) {
641     Double_t b1 = (2*i + 1) * kB;
642     cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
643   }
644   cumulPRF *= kA/(2*TMath::Pi());
645   return cumulPRF;
646 }
647