X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PHOS%2FAliPHOSv1.cxx;h=b99e522710377f6031b0aae1b98de7515809c64a;hp=0e2464fab7375e8c5d85f047c783a0b42d575abf;hb=6bcdf1d1c0ec8e749cf34d632a4a8facd385645f;hpb=8dfa469d2de8f9b361ba598614b73e01876e5b1c diff --git a/PHOS/AliPHOSv1.cxx b/PHOS/AliPHOSv1.cxx index 0e2464fab73..b99e5227103 100644 --- a/PHOS/AliPHOSv1.cxx +++ b/PHOS/AliPHOSv1.cxx @@ -14,257 +14,550 @@ **************************************************************************/ /* $Id$ */ + +/* History of cvs commits: + * + * $Log$ + * Revision 1.111 2007/07/24 09:41:19 morsch + * AliStack included for kKeepBit. + * + * Revision 1.110 2007/03/10 08:58:52 kharlov + * Protection for noCPV geometry + * + * Revision 1.109 2007/03/01 11:37:37 kharlov + * Strip units changed from 8x1 to 8x2 (T.Pocheptsov) + * + * Revision 1.108 2007/02/02 09:40:50 alibrary + * Includes required by ROOT head + * + * Revision 1.107 2007/02/01 10:34:47 hristov + * Removing warnings on Solaris x86 + * + * Revision 1.106 2006/11/14 17:11:15 hristov + * 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 + * + * Revision 1.105 2006/09/13 07:31:01 kharlov + * Effective C++ corrections (T.Pocheptsov) + * + * Revision 1.104 2005/05/28 14:19:05 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ -// Version of AliPHOSv0 which allows for keeping all hits in TreeH -// -//*-- Author: Gines MARTINEZ (SUBATECH) +// Implementation version v1 of PHOS Manager class +//--- +//--- +// Layout EMC + CPV has name IHEP: +// Produces hits for CPV, cumulated hits +//--- +//--- +//*-- Author: Yves Schutz (SUBATECH) // --- ROOT system --- - -#include "TBRIK.h" -#include "TNode.h" -#include "TRandom.h" +#include +#include +#include // --- Standard library --- -#include -#include -#include -#include // --- AliRoot header files --- - -#include "AliPHOSv1.h" +#include "AliPHOSCPVDigit.h" +#include "AliPHOSGeometry.h" #include "AliPHOSHit.h" -#include "AliPHOSDigit.h" -#include "AliPHOSReconstructioner.h" +#include "AliPHOSv1.h" #include "AliRun.h" -#include "AliConst.h" +#include "AliMC.h" +#include "AliStack.h" +#include "AliPHOSSimParam.h" ClassImp(AliPHOSv1) //____________________________________________________________________________ -AliPHOSv1::AliPHOSv1() +AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20) { - // default ctor - fNTmpHits = 0 ; - fTmpHits = 0 ; + //Def ctor. } //____________________________________________________________________________ AliPHOSv1::AliPHOSv1(const char *name, const char *title): - AliPHOSv0(name,title) + AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20) { - // ctor - fHits= new TClonesArray("AliPHOSHit",1000) ; + // + // We store hits : + // - fHits (the "normal" one), which retains the hits associated with + // the current primary particle being tracked + // (this array is reset after each primary has been tracked). + // + + + + // We do not want to save in TreeH the raw hits + // But save the cumulated hits instead (need to create the branch myself) + // It is put in the Digit Tree because the TreeH is filled after each primary + // and the TreeD at the end of the event (branch is set in FinishEvent() ). + + fHits= new TClonesArray("AliPHOSHit",1000) ; +// fCPVDigits("AliPHOSCPVDigit",20); + gAlice->GetMCApp()->AddHitList(fHits) ; + + fNhits = 0 ; + + fIshunt = 2 ; // All hits are associated with primary particles } //____________________________________________________________________________ AliPHOSv1::~AliPHOSv1() { // dtor - if ( fTmpHits) { - fTmpHits->Delete() ; - delete fTmpHits ; - fTmpHits = 0 ; - } - - if ( fEmcRecPoints ) { - fEmcRecPoints->Delete() ; - delete fEmcRecPoints ; - fEmcRecPoints = 0 ; - } - - if ( fPpsdRecPoints ) { - fPpsdRecPoints->Delete() ; - delete fPpsdRecPoints ; - fPpsdRecPoints = 0 ; - } - - if ( fTrackSegments ) { - fTrackSegments->Delete() ; - delete fTrackSegments ; - fTrackSegments = 0 ; - } + if ( fHits) { + fHits->Delete() ; + delete fHits ; + fHits = 0 ; + } } //____________________________________________________________________________ -void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits) +void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits) { // Add a hit to the hit list. - // In this version of AliPHOSv0, a PHOS hit is real geant - // hits in a single crystal or in a single PPSD gas cell - - // cout << "Primary is " << primary << endl; - //cout << "Tracknumber is " << tracknumber << endl; - //cout << "Vol Id is " << Id << endl; - //cout << "hits is " << hits[0] << " " << hits[1] << " " << hits[2] << " " << hits[3] <= 0 && !deja ; hitCounter-- ) { + curHit = dynamic_cast((*fHits)[hitCounter]) ; + if(curHit->GetPrimary() != primary) break ; + // We add hits with the same primary, while GEANT treats primaries succesively + if( *curHit == *newHit ) { + *curHit + *newHit ; + deja = kTRUE ; + } + } + + if ( !deja ) { + new((*fHits)[fNhits]) AliPHOSHit(*newHit) ; + // get the block Id number + Int_t relid[4] ; + geom->AbsToRelNumbering(Id, relid) ; + + fNhits++ ; + } + + delete newHit; +} - // fHits->Print(""); - - newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ; - - // We DO want to save in TreeH the raw hits - // TClonesArray &lhits = *fHits; - // cout << " Adding a hit to fHits TCloneArray number " << fNhits << endl ; - new(ltmphits[fNhits]) AliPHOSHit(*newHit) ; - fNhits++ ; +//____________________________________________________________________________ +void AliPHOSv1::FinishPrimary() +{ + // called at the end of each track (primary) by AliRun + // hits are reset for each new track + // accumulate the total hit-multiplicity - // cout << " Added a hit to fHits TCloneArray number " << fNhits << endl ; - // We do not want to save in TreeH the raw hits - // new(lhits[fNhits]) AliPHOSHit(*newHit) ; - // fNhits++ ; +} - // Please note that the fTmpHits array must survive up to the - // end of the events, so it does not appear e.g. in ResetHits() ( - // which is called at the end of each primary). +//____________________________________________________________________________ +void AliPHOSv1::FinishEvent() +{ + // called at the end of each event by AliRun + // accumulate the hit-multiplicity and total energy per block + // if the values have been updated check it + + AliDetector::FinishEvent(); +} +//____________________________________________________________________________ +void AliPHOSv1::StepManager(void) +{ + // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell - delete newHit; + Int_t relid[4] ; // (box, layer, row, column) indices + Int_t absid ; // absolute cell ID number + Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited + TLorentzVector pos ; // Lorentz vector of the track current position + Int_t copy ; -} + Int_t moduleNumber ; + + static Int_t idPCPQ = -1; + if (strstr(fTitle.Data(),"noCPV") == 0) + idPCPQ = gMC->VolId("PCPQ"); + if( gMC->CurrentVolID(copy) == idPCPQ && + (gMC->IsTrackEntering() ) && + gMC->TrackCharge() != 0) { + + gMC -> TrackPosition(pos); + Float_t xyzm[3], xyzd[3] ; + Int_t i; + for (i=0; i<3; i++) xyzm[i] = pos[i]; + gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system -//___________________________________________________________________________ -void AliPHOSv1::FinishEvent() -{ - // Makes the digits from the sum of summed hit in a single crystal or PPSD gas cell - // Adds to the energy the electronic noise - // Keeps digits with energy above fDigitThreshold + Float_t xyd[3]={0,0,0} ; //local position of the entering + xyd[0] = xyzd[0]; + xyd[1] =-xyzd[2]; + xyd[2] =-xyzd[1]; + + // Current momentum of the hit's track in the local ref. system + TLorentzVector pmom ; //momentum of the particle initiated hit + gMC -> TrackMomentum(pmom); + Float_t pm[3], pd[3]; + for (i=0; i<3; i++) + pm[i] = pmom[i]; + + gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system + pmom[0] = pd[0]; + pmom[1] =-pd[1]; + pmom[2] =-pd[2]; - // Save the cumulated hits instead of raw hits (need to create the branch myself) - // It is put in the Digit Tree because the TreeH is filled after each primary - // and the TreeD at the end of the event. - - Int_t i ; - Int_t relid[4]; - Int_t j ; - TClonesArray &lDigits = *fDigits ; - AliPHOSHit * hit ; - AliPHOSDigit * newdigit ; - AliPHOSDigit * curdigit ; - Bool_t deja = kFALSE ; - - for ( i = 0 ; i < fNhits ; i++ ) { - hit = (AliPHOSHit*)fHits->At(i) ; - newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ; - deja =kFALSE ; - for ( j = 0 ; j < fNdigits ; j++) { - curdigit = (AliPHOSDigit*) lDigits[j] ; - if ( *curdigit == *newdigit) { - *curdigit = *curdigit + *newdigit ; - deja = kTRUE ; + // Digitize the current CPV hit: + + // 1. find pad response and + gMC->CurrentVolOffID(3,moduleNumber); + moduleNumber--; + +// TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit + CPVDigitize(pmom,xyd,&fCPVDigits); + + Float_t xmean = 0; + Float_t zmean = 0; + Float_t qsum = 0; + Int_t idigit,ndigits; + + // 2. go through the current digit list and sum digits in pads + + ndigits = fCPVDigits.GetEntriesFast(); + for (idigit=0; idigit(fCPVDigits.UncheckedAt(idigit)); + Float_t x1 = cpvDigit1->GetXpad() ; + Float_t z1 = cpvDigit1->GetYpad() ; + for (Int_t jdigit=idigit+1; jdigit(fCPVDigits.UncheckedAt(jdigit)); + Float_t x2 = cpvDigit2->GetXpad() ; + Float_t z2 = cpvDigit2->GetYpad() ; + if (x1==x2 && z1==z2) { + Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ; + cpvDigit2->SetQpad(qsumpad) ; + fCPVDigits.RemoveAt(idigit) ; + } } } - if ( !deja ) { - new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ; - fNdigits++ ; + fCPVDigits.Compress() ; + + // 3. add digits to temporary hit list fTmpHits + + ndigits = fCPVDigits.GetEntriesFast(); + for (idigit=0; idigit(fCPVDigits.UncheckedAt(idigit)); + relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number + relid[1] =-1 ; // means CPV + relid[2] = cpvDigit->GetXpad() ; // column number of a pad + relid[3] = cpvDigit->GetYpad() ; // row number of a pad + + // get the absolute Id number + GetGeometry()->RelToAbsNumbering(relid, absid) ; + + // add current digit to the temporary hit list + + xyzte[3] = gMC->TrackTime() ; + xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad + + Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); + AddHit(fIshunt, primary, absid, xyzte); + + if (cpvDigit->GetQpad() > 0.02) { + xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5); + zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5); + qsum += cpvDigit->GetQpad(); + } } - - delete newdigit ; - } - - // Noise induced by the PIN diode of the PbWO crystals - - Float_t energyandnoise ; - for ( i = 0 ; i < fNdigits ; i++ ) { - newdigit = (AliPHOSDigit * ) fDigits->At(i) ; - fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ; + fCPVDigits.Clear(); + } - if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals) - energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ; + + static Int_t idPXTL = gMC->VolId("PXTL"); + if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal - if (energyandnoise < 0 ) - energyandnoise = 0 ; + gMC->TrackPosition(pos) ; + xyzte[0] = pos[0] ; + xyzte[1] = pos[1] ; + xyzte[2] = pos[2] ; - if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list - fDigits->RemoveAt(i) ; + Float_t lostenergy = gMC->Edep(); + + //Put in the TreeK particle entering PHOS and all its parents + if ( gMC->IsTrackEntering() ){ + Float_t xyzd[3] ; + gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system + if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface + Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ; + TParticle * part = gAlice->GetMCApp()->Particle(parent) ; + Float_t vert[3],vertd[3] ; + vert[0]=part->Vx() ; + vert[1]=part->Vy() ; + vert[2]=part->Vz() ; + gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system + if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS + //0.1 to get rid of numerical errors + part->SetBit(kKeepBit); + while ( parent != -1 ) { + part = gAlice->GetMCApp()->Particle(parent) ; + part->SetBit(kKeepBit); + parent = part->GetFirstMother() ; + } + } + } } - } + if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy + xyzte[3] = gMC->TrackTime() ; + + gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ; + + Int_t strip ; + gMC->CurrentVolOffID(3, strip); + Int_t cell ; + gMC->CurrentVolOffID(2, cell); + + //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1. + //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29) + //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ; + Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ; + Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ; + + // Absid for 8x2-strips. Looks nice :) + absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + + row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0); + + + //Calculates the light yield, the number of photons produced in the + //crystal + //There is no dependence of reponce on distance from energy deposition to APD + Float_t lightYield = gRandom->Poisson(AliPHOSSimParam::GetInstance()->GetLightFactor() * lostenergy) ; + + //Calculates de energy deposited in the crystal + xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ; + + Int_t primary ; + if(fIshunt == 2){ + primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ; + TParticle * part = gAlice->GetMCApp()->Particle(primary) ; + while ( !part->TestBit(kKeepBit) ) { + primary = part->GetFirstMother() ; + if(primary == -1){ + primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); + break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side + //surface of the crystal. In this case it may have no primary at all. + //We can not easily separate this case from the case when this is part of the shower, + //developed in the neighboring crystal. + } + part = gAlice->GetMCApp()->Particle(primary) ; + } + } + else{ + primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); + } + + // add current hit to the hit list + // Info("StepManager","%d %d", primary, tracknumber) ; + AddHit(fIshunt, primary, absid, xyzte); + + } // there is deposited energy + } // we are inside a PHOS Xtal - fDigits->Compress() ; - - fNdigits = fDigits->GetEntries() ; - for (i = 0 ; i < fNdigits ; i++) { - newdigit = (AliPHOSDigit *) fDigits->At(i) ; - newdigit->SetIndexInList(i) ; - } - } -void AliPHOSv1::StepManager(void) +//____________________________________________________________________________ +void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits) { - // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell - - Int_t relid[4] ; // (box, layer, row, column) indices - Float_t xyze[4] ; // position wrt MRS and energy deposited - TLorentzVector pos ; - Int_t copy ; - - Int_t tracknumber = gAlice->CurrentTrack() ; - Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() ); - - TString name = fGeom->GetName() ; - if ( name == "GPS2" ) { // the CPV is a PPSD - if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell - { - gMC->TrackPosition(pos) ; - xyze[0] = pos[0] ; - xyze[1] = pos[1] ; - xyze[2] = pos[2] ; - xyze[3] = gMC->Edep() ; - - if ( xyze[3] != 0 ) { // there is deposited energy - gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number - gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number - // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper - // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower - gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell - gMC->CurrentVolID(relid[3]) ; // get the column number - - // get the absolute Id number - - Int_t absid ; - fGeom->RelToAbsNumbering(relid, absid) ; - - // add current hit to the hit list - AddHit(fIshunt, primary, tracknumber, absid, xyze); - - } // there is deposited energy - } // We are inside the gas of the CPV - } // GPS2 configuration - - 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 + // ------------------------------------------------------------------------ + // 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 kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad) + const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e- + const Float_t kdEdx = 4.0; // Average energy loss in CPV; + const Int_t kNgamz = 5; // Ionization size in Z + const Int_t kNgamx = 9; // Ionization size in Phi + const Float_t kNoise = 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 + + 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 eloss = kdEdx; + + Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness(); + Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness(); + gRandom->Rannor(rnor1,rnor2); + eloss *= (1 + kDetR*rnor1) * + TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2))); + Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2; + Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2; + Float_t zhit2 = zhit1 + dZY; + Float_t xhit2 = xhit1 + dXY; + + Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in" + Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // 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) * kCelWr; + zxe[2][0] = eloss/2; + zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2; + zxe[1][1] = (iwht1 + 0.5) * kCelWr; + zxe[2][1] = eloss/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) * kCelWr; // wire 1 + Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2 + Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // 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) + kCelWr ); + Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr ); + Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr ); + zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2; + zxe[1][0] = xwht1; + zxe[2][0] = eloss * egm1; + zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2; + zxe[1][1] = xwht2; + zxe[2][1] = eloss * egm2; + zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1; + zxe[1][2] = xwht3; + zxe[2][2] = eloss * egm3; + } + else { // incline 2-wire hit + nIter = 2; + Float_t xwht1 = (iwht1 + 0.5) * kCelWr; + Float_t xwht2 = (iwht2 + 0.5) * kCelWr; + 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] = eloss * egm1; + zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2; + zxe[1][1] = xwht2; + zxe[2][1] = eloss * egm2; + } - // get the absolute Id number + // Finite size of ionization region + + Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ(); + Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi(); + Int_t nz3 = (kNgamz+1)/2; + Int_t nx3 = (kNgamx+1)/2; + cpvDigits->Expand(nIter*kNgamx*kNgamz); + TClonesArray &ldigits = *(static_cast(cpvDigits)); + + for (Int_t iter=0; iterGetPadSizeZ(); + Float_t xcell = xhit / GetGeometry()->GetPadSizePhi(); + if ( zcell<=0 || xcell<=0 || + zcell>=nCellZ || xcell>=nCellX) return; + Int_t izcell = (Int_t) zcell; + Int_t ixcell = (Int_t) xcell; + Float_t zc = zcell - izcell - 0.5; + Float_t xc = xcell - ixcell - 0.5; + for (Int_t iz=1; iz<=kNgamz; iz++) { + Int_t kzg = izcell + iz - nz3; + if (kzg<=0 || kzg>nCellZ) continue; + Float_t zg = (Float_t)(iz-nz3) - zc; + for (Int_t ix=1; ix<=kNgamx; ix++) { + Int_t kxg = ixcell + ix - nx3; + if (kxg<=0 || kxg>nCellX) continue; + Float_t xg = (Float_t)(ix-nx3) - xc; + + // Now calculate pad response + Float_t qpad = CPVPadResponseFunction(qhit,zg,xg); + qpad += kNoise*rnor2; + if (qpad<0) continue; + + // Fill the array with pad response ID and amplitude + new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad); + } + } + } +} - Int_t absid ; - fGeom->RelToAbsNumbering(relid, absid) ; - - // add current hit to the hit list +//____________________________________________________________________________ +Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) { + // ------------------------------------------------------------------------ + // Calculate the amplitude in one CPV pad using the + // cumulative pad response function + // Author: Yuri Kharlov (after Serguei Sadovski) + // 3 October 2000 + // ------------------------------------------------------------------------ + + Double_t dz = GetGeometry()->GetPadSizeZ() / 2; + Double_t dx = GetGeometry()->GetPadSizePhi() / 2; + Double_t z = zhit * GetGeometry()->GetPadSizeZ(); + Double_t x = xhit * GetGeometry()->GetPadSizePhi(); + Double_t amplitude = qhit * + (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) - + CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx)); + return (Float_t)amplitude; +} - AddHit(fIshunt, primary,tracknumber, absid, xyze); - - } // there is deposited energy - } // we are inside a PHOS Xtal +//____________________________________________________________________________ +Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) { + // ------------------------------------------------------------------------ + // Cumulative pad response function + // It includes several terms from the CF decomposition in electrostatics + // Note: this cumulative function is wrong since omits some terms + // but the cell amplitude obtained with it is correct because + // these omitting terms cancel + // Author: Yuri Kharlov (after Serguei Sadovski) + // 3 October 2000 + // ------------------------------------------------------------------------ + + const Double_t kA=1.0; + const Double_t kB=0.7; + + Double_t r2 = x*x + y*y; + Double_t xy = x*y; + Double_t cumulPRF = 0; + for (Int_t i=0; i<=4; i++) { + Double_t b1 = (2*i + 1) * kB; + cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) ); + } + cumulPRF *= kA/(2*TMath::Pi()); + return cumulPRF; }