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