/* 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
*
#include "AliPHOSv1.h"
#include "AliRun.h"
#include "AliMC.h"
+#include "AliStack.h"
+#include "AliPHOSSimParam.h"
ClassImp(AliPHOSv1)
//____________________________________________________________________________
-AliPHOSv1::AliPHOSv1():
- fLightYieldMean(0.),
- fIntrinsicPINEfficiency(0.),
- fLightYieldAttenuation(0.),
- fRecalibrationFactor(0.),
- fElectronsPerGeV(0.),
- fAPDGain(0.),
- fLightFactor(0.),
- fAPDFactor(0.)
+AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20)
{
//Def ctor.
}
//____________________________________________________________________________
AliPHOSv1::AliPHOSv1(const char *name, const char *title):
- AliPHOSv0(name,title),
- fLightYieldMean(0.),
- fIntrinsicPINEfficiency(0.),
- fLightYieldAttenuation(0.),
- fRecalibrationFactor(0.),
- fElectronsPerGeV(0.),
- fAPDGain(0.),
- fLightFactor(0.),
- fAPDFactor(0.)
+ AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20)
{
//
// We store hits :
// 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
-
- //Photoelectron statistics:
- // The light yield is a poissonian distribution of the number of
- // photons created in the PbWo4 crystal, calculated using following formula
- // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
- // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
- // LightYieldMean is parameter calculated to be over 47000 photons per GeV
- // APDEfficiency is 0.02655
- // k_0 is 0.0045 from Valery Antonenko
- // The number of electrons created in the APD is
- // NumberOfElectrons = APDGain * LightYield
- // The APD Gain is 300
- fLightYieldMean = 47000;
- fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
- fLightYieldAttenuation = 0.0045 ;
- fRecalibrationFactor = 13.418/ fLightYieldMean ;
- fElectronsPerGeV = 2.77e+8 ;
- fAPDGain = 300. ;
- fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
- fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
}
//____________________________________________________________________________
newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
- curHit =
dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
+ curHit =
static_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
if(curHit->GetPrimary() != primary) break ;
// We add hits with the same primary, while GEANT treats primaries succesively
if( *curHit == *newHit ) {
TLorentzVector pos ; // Lorentz vector of the track current position
Int_t copy ;
- TString name = GetGeometry()->GetName() ;
-
Int_t moduleNumber ;
- static Int_t idPCPQ = gMC->VolId("PCPQ");
- if( gMC->CurrentVolID(copy) == idPCPQ &&
- (gMC->IsTrackEntering() ) &&
- gMC->TrackCharge() != 0) {
-
- gMC -> TrackPosition(pos);
+ static Int_t idPCPQ = -1;
+ if (strstr(fTitle.Data(),"noCPV") == 0)
+ idPCPQ = TVirtualMC::GetMC()->VolId("PCPQ");
+
+ if( TVirtualMC::GetMC()->CurrentVolID(copy) == idPCPQ &&
+ (TVirtualMC::GetMC()->IsTrackEntering() ) &&
+ TVirtualMC::GetMC()->TrackCharge() != 0) {
+ TVirtualMC::GetMC() -> 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
-
+ TVirtualMC::GetMC() -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
+
+
Float_t xyd[3]={0,0,0} ; //local position of the entering
xyd[0] = xyzd[0];
xyd[1] =-xyzd[2];
// Current momentum of the hit's track in the local ref. system
TLorentzVector pmom ; //momentum of the particle initiated hit
- gMC -> TrackMomentum(pmom);
+ TVirtualMC::GetMC() -> 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
+ TVirtualMC::GetMC() -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
pmom[0] = pd[0];
pmom[1] =-pd[1];
pmom[2] =-pd[2];
// Digitize the current CPV hit:
// 1. find pad response and
- gMC->CurrentVolOffID(3,moduleNumber);
+ TVirtualMC::GetMC()->CurrentVolOffID(3,moduleNumber);
moduleNumber--;
- TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
- CPVDigitize(pmom,xyd,cpvDigits);
+// TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
+ CPVDigitize(pmom,xyd,&fCPVDigits);
Float_t xmean = 0;
Float_t zmean = 0;
// 2. go through the current digit list and sum digits in pads
- ndigits = cpvDigits->GetEntriesFast();
+ ndigits = fCPVDigits.GetEntriesFast();
for (idigit=0; idigit<ndigits-1; idigit++) {
- AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+ AliPHOSCPVDigit *cpvDigit1 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
Float_t x1 = cpvDigit1->GetXpad() ;
Float_t z1 = cpvDigit1->GetYpad() ;
for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
- AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
+ AliPHOSCPVDigit *cpvDigit2 = static_cast<AliPHOSCPVDigit*>(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) ;
- cpvDigits->RemoveAt(idigit) ;
+ fCPVDigits.RemoveAt(idigit) ;
}
}
}
- cpvDigits->Compress() ;
+ fCPVDigits.Compress() ;
// 3. add digits to temporary hit list fTmpHits
- ndigits = cpvDigits->GetEntriesFast();
+ ndigits = fCPVDigits.GetEntriesFast();
for (idigit=0; idigit<ndigits; idigit++) {
- AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+ AliPHOSCPVDigit *cpvDigit = static_cast<AliPHOSCPVDigit*>(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
// add current digit to the temporary hit list
- xyzte[3] = gMC->TrackTime() ;
+ xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
qsum += cpvDigit->GetQpad();
}
}
- if (cpvDigits) {
- cpvDigits->Delete();
- delete cpvDigits;
- cpvDigits=0;
- }
+ fCPVDigits.Clear();
}
- static Int_t idPXTL = gMC->VolId("PXTL");
- if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
+ static Int_t idPXTL = TVirtualMC::GetMC()->VolId("PXTL");
+ if(TVirtualMC::GetMC()->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
- gMC->TrackPosition(pos) ;
+ TVirtualMC::GetMC()->TrackPosition(pos) ;
xyzte[0] = pos[0] ;
xyzte[1] = pos[1] ;
xyzte[2] = pos[2] ;
- Float_t global[3], local[3] ;
- global[0] = pos[0] ;
- global[1] = pos[1] ;
- global[2] = pos[2] ;
- Float_t lostenergy = gMC->Edep();
+ Float_t lostenergy = TVirtualMC::GetMC()->Edep();
//Put in the TreeK particle entering PHOS and all its parents
- if ( gMC->IsTrackEntering() ){
+ if ( TVirtualMC::GetMC()->IsTrackEntering() ){
Float_t xyzd[3] ;
- gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
+ TVirtualMC::GetMC() -> 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) ;
vert[0]=part->Vx() ;
vert[1]=part->Vy() ;
vert[2]=part->Vz() ;
- gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
+ TVirtualMC::GetMC() -> 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);
}
}
if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
- xyzte[3] = gMC->TrackTime() ;
+ xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
- gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
+ TVirtualMC::GetMC()->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
Int_t strip ;
- gMC->CurrentVolOffID(3, strip);
+ TVirtualMC::GetMC()->CurrentVolOffID(3, strip);
Int_t cell ;
- gMC->CurrentVolOffID(2, cell);
-
- Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
- Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
-
+ TVirtualMC::GetMC()->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 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
-
- gMC->Gmtod(global, local, 1) ;
+ 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
- Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
- exp(-fLightYieldAttenuation *
- (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
- ) ;
+ //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] = fAPDFactor * lightYield ;
+ xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ;
Int_t primary ;
if(fIshunt == 2){
part = gAlice->GetMCApp()->Particle(primary) ;
}
}
- else
+ else{
primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
-
-
+ }
// add current hit to the hit list
// Info("StepManager","%d %d", primary, tracknumber) ;
Float_t pNorm = p.Py();
Float_t eloss = kdEdx;
-//Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
-
Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
gRandom->Rannor(rnor1,rnor2);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff