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