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