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