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