+ if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
+ gMC->TrackPosition(pos) ;
+ xyze[0] = pos[0] ;
+ xyze[1] = pos[1] ;
+ xyze[2] = pos[2] ;
+ xyze[3] = gMC->Edep() ;
+
+ if ( xyze[3] != 0 ) {
+ gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
+ relid[1] = 0 ; // means PBW04
+ gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
+ gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
+
+ // get the absolute Id number
+
+ fGeom->RelToAbsNumbering(relid, absid) ;
+
+ // add current hit to the hit list
+
+ AddHit(fIshunt, primary,tracknumber, absid, xyze);
+
+ } // there is deposited energy
+ } // we are inside a PHOS Xtal
+}
+
+//____________________________________________________________________________
+void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
+{
+ // ------------------------------------------------------------------------
+ // Digitize one CPV hit:
+ // On input take exact 4-momentum p and position zxhit of the hit,
+ // find the pad response around this hit and
+ // put the amplitudes in the pads into array digits
+ //
+ // Author: Yuri Kharlov (after Serguei Sadovsky)
+ // 2 October 2000
+ // ------------------------------------------------------------------------
+
+ const Float_t celWr = fGeom->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
+ const Float_t detR = 0.1; // Relative energy fluctuation in track for 100 e-
+ const Float_t dEdx = 4.0; // Average energy loss in CPV;
+ const Int_t ngamz = 5; // Ionization size in Z
+ const Int_t ngamx = 9; // Ionization size in Phi
+ const Float_t qNoise = 0.03; // charge noise in one pad
+
+ Float_t rnor1,rnor2;
+
+ // Just a reminder on axes notation in the CPV module:
+ // axis Z goes along the beam
+ // axis X goes across the beam in the module plane
+ // axis Y is a normal to the module plane showing from the IP
+
+// cout << __PRETTY_FUNCTION__ << ": YVK : Start digitization\n";
+
+ Float_t hitX = zxhit[0];
+ Float_t hitZ =-zxhit[1];
+ Float_t pX = p.Px();
+ Float_t pZ =-p.Pz();
+ Float_t pNorm = p.Py();
+ Float_t E = dEdx;
+
+// cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
+
+ Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
+ Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
+ gRandom->Rannor(rnor1,rnor2);
+ E *= (1 + detR*rnor1) *
+ TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
+ Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
+ Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
+ Float_t zhit2 = zhit1 + dZY;
+ Float_t xhit2 = xhit1 + dXY;
+
+// cout << "CPVDigitize: YVK : "<<xhit1<<" "<<xhit2<<" | "<<zhit1<<" "<<zhit2<<" | "
+// << (Int_t)(xhit1/fGeom->GetPadSizePhi())<<" "<<(Int_t)(zhit1/fGeom->GetPadSizeZ())<<" "<<endl;
+
+ Int_t iwht1 = (Int_t) (xhit1 / celWr); // wire (x) coordinate "in"
+ Int_t iwht2 = (Int_t) (xhit2 / celWr); // wire (x) coordinate "out"
+
+ Int_t nIter;
+ Float_t zxe[3][5];
+ if (iwht1==iwht2) { // incline 1-wire hit
+ nIter = 2;
+ zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
+ zxe[1][0] = (iwht1 + 0.5) * celWr;
+ zxe[2][0] = E/2;
+ zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
+ zxe[1][1] = (iwht1 + 0.5) * celWr;
+ zxe[2][1] = E/2;
+ }
+ else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
+ nIter = 3;
+ Int_t iwht3 = (iwht1 + iwht2) / 2;
+ Float_t xwht1 = (iwht1 + 0.5) * celWr; // wire 1
+ Float_t xwht2 = (iwht2 + 0.5) * celWr; // wire 2
+ Float_t xwht3 = (iwht3 + 0.5) * celWr; // wire 3
+ Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
+ Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
+ Float_t dxw1 = xhit1 - xwr13;
+ Float_t dxw2 = xhit2 - xwr23;
+ Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
+ Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
+ Float_t egm3 = celWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
+ zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
+ zxe[1][0] = xwht1;
+ zxe[2][0] = E * egm1;
+ zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
+ zxe[1][1] = xwht2;
+ zxe[2][1] = E * egm2;
+ zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
+ zxe[1][2] = xwht3;
+ zxe[2][2] = E * egm3;
+ }
+ else { // incline 2-wire hit
+ nIter = 2;
+ Float_t xwht1 = (iwht1 + 0.5) * celWr;
+ Float_t xwht2 = (iwht2 + 0.5) * celWr;
+ Float_t xwr12 = (xwht1 + xwht2) / 2;
+ Float_t dxw1 = xhit1 - xwr12;
+ Float_t dxw2 = xhit2 - xwr12;
+ Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
+ Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
+ zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
+ zxe[1][0] = xwht1;
+ zxe[2][0] = E * egm1;
+ zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
+ zxe[1][1] = xwht2;
+ zxe[2][1] = E * egm2;
+ }