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