X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONResponseV0.cxx;h=d47dd9aa05f37ede279846c69a03b4d3080db6d1;hb=e475373795a3db0246740b19449e2577400ff03b;hp=496f7061360087113f132395639a2aba724ee203;hpb=f7b62f08b840171d1d540bacf94bba090b04c91d;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONResponseV0.cxx b/MUON/AliMUONResponseV0.cxx index 496f7061360..d47dd9aa05f 100644 --- a/MUON/AliMUONResponseV0.cxx +++ b/MUON/AliMUONResponseV0.cxx @@ -15,45 +15,167 @@ /* $Id$ */ +//----------------------------------------------------------------------------- +// Class AliMUONResponseV0 +// -------------------------- +// Implementation of +// Mathieson response +//----------------------------------------------------------------------------- + #include "AliMUONResponseV0.h" -#include "AliSegmentation.h" +#include "AliMUON.h" +#include "AliMUONConstants.h" +#include "AliMUONDigit.h" +#include "AliMUONGeometryTransformer.h" +#include "AliMUONHit.h" + +#include "AliMpArea.h" +#include "AliMpDEManager.h" +#include "AliMpVPadIterator.h" +#include "AliMpSegmentation.h" +#include "AliMpVSegmentation.h" +#include "AliMpCathodType.h" + +#include "AliRun.h" +#include "AliLog.h" + +#include "Riostream.h" +#include "TVector2.h" #include #include - +/// \cond CLASSIMP ClassImp(AliMUONResponseV0) +/// \endcond +AliMUON* muon() +{ + return static_cast(gAlice->GetModule("MUON")); +} + +void Global2Local(Int_t detElemId, Double_t xg, Double_t yg, Double_t zg, + Double_t& xl, Double_t& yl, Double_t& zl) +{ + /// ideally should be : + /// Double_t x,y,z; + /// AliMUONGeometry::Global2Local(detElemId,xg,yg,zg,x,y,z); + /// but while waiting for this geometry singleton, let's go through + /// AliMUON still. + + const AliMUONGeometryTransformer* transformer = muon()->GetGeometryTransformer(); + transformer->Global2Local(detElemId,xg,yg,zg,xl,yl,zl); +} + +//__________________________________________________________________________ +AliMUONResponseV0::AliMUONResponseV0() + : AliMUONResponse(), + fChargeSlope(0.0), + fChargeSpreadX(0.0), + fChargeSpreadY(0.0), + fSigmaIntegration(0.0), + fMaxAdc(0), + fSaturation(0), + fZeroSuppression(0), + fChargeCorrel(0.0), + fMathieson(new AliMUONMathieson), + fChargeThreshold(1e-4) +{ + /// Normal constructor + AliDebug(1,Form("Default ctor")); +} + +//__________________________________________________________________________ +AliMUONResponseV0::AliMUONResponseV0(const AliMUONResponseV0& other) +: AliMUONResponse(), +fChargeSlope(0.0), +fChargeSpreadX(0.0), +fChargeSpreadY(0.0), +fSigmaIntegration(0.0), +fMaxAdc(0), +fSaturation(0), +fZeroSuppression(0), +fChargeCorrel(0.0), +fMathieson(0), +fChargeThreshold(1e-4) +{ + /// copy ctor + other.CopyTo(*this); +} + +//__________________________________________________________________________ +AliMUONResponseV0& +AliMUONResponseV0::operator=(const AliMUONResponseV0& other) +{ + /// Assignment operator + other.CopyTo(*this); + return *this; +} + +//__________________________________________________________________________ +AliMUONResponseV0::~AliMUONResponseV0() +{ +/// Destructor + + AliDebug(1,""); + delete fMathieson; +} + +//______________________________________________________________________________ +void +AliMUONResponseV0::CopyTo(AliMUONResponseV0& other) const +{ + /// Copy *this to other + other.fChargeSlope=fChargeSlope; + other.fChargeSpreadX=fChargeSpreadX; + other.fChargeSpreadY=fChargeSpreadY; + other.fSigmaIntegration=fSigmaIntegration; + other.fMaxAdc=fMaxAdc; + other.fSaturation=fSaturation; + other.fZeroSuppression=fZeroSuppression; + other.fChargeCorrel=fChargeCorrel; + delete other.fMathieson; + other.fMathieson = new AliMUONMathieson(*fMathieson); + other.fChargeThreshold=fChargeThreshold; +} + +//______________________________________________________________________________ +void +AliMUONResponseV0::Print(Option_t*) const +{ +/// Printing + + cout << " ChargeSlope=" << fChargeSlope + << " ChargeSpreadX,Y=" << fChargeSpreadX + << fChargeSpreadY + << " ChargeCorrelation=" << fChargeCorrel + << endl; +} + //__________________________________________________________________________ void AliMUONResponseV0::SetSqrtKx3AndDeriveKx2Kx4(Float_t SqrtKx3) { - // Set to "SqrtKx3" the Mathieson parameter K3 ("fSqrtKx3") - // in the X direction, perpendicular to the wires, - // and derive the Mathieson parameters K2 ("fKx2") and K4 ("fKx4") - // in the same direction - fSqrtKx3 = SqrtKx3; - fKx2 = TMath::Pi() / 2. * (1. - 0.5 * fSqrtKx3); - Float_t cx1 = fKx2 * fSqrtKx3 / 4. / TMath::ATan(Double_t(fSqrtKx3)); - fKx4 = cx1 / fKx2 / fSqrtKx3; + /// Set to "SqrtKx3" the Mathieson parameter K3 ("fSqrtKx3") + /// in the X direction, perpendicular to the wires, + /// and derive the Mathieson parameters K2 ("fKx2") and K4 ("fKx4") + /// in the same direction + fMathieson->SetSqrtKx3AndDeriveKx2Kx4(SqrtKx3); } //__________________________________________________________________________ void AliMUONResponseV0::SetSqrtKy3AndDeriveKy2Ky4(Float_t SqrtKy3) { - // Set to "SqrtKy3" the Mathieson parameter K3 ("fSqrtKy3") - // in the Y direction, along the wires, - // and derive the Mathieson parameters K2 ("fKy2") and K4 ("fKy4") - // in the same direction - fSqrtKy3 = SqrtKy3; - fKy2 = TMath::Pi() / 2. * (1. - 0.5 * fSqrtKy3); - Float_t cy1 = fKy2 * fSqrtKy3 / 4. / TMath::ATan(Double_t(fSqrtKy3)); - fKy4 = cy1 / fKy2 / fSqrtKy3; + /// Set to "SqrtKy3" the Mathieson parameter K3 ("fSqrtKy3") + /// in the Y direction, along the wires, + /// and derive the Mathieson parameters K2 ("fKy2") and K4 ("fKy4") + /// in the same direction + fMathieson->SetSqrtKy3AndDeriveKy2Ky4(SqrtKy3); } - -Float_t AliMUONResponseV0::IntPH(Float_t eloss) + //__________________________________________________________________________ +Float_t AliMUONResponseV0::IntPH(Float_t eloss) const { - // Calculate charge from given ionization energy loss + /// Calculate charge from given ionization energy loss Int_t nel; - nel= Int_t(eloss*1.e9/32.); + nel= Int_t(eloss*1.e9/27.4); Float_t charge=0; if (nel == 0) nel=1; for (Int_t i=1;i<=nel;i++) { @@ -63,54 +185,118 @@ Float_t AliMUONResponseV0::IntPH(Float_t eloss) } return charge; } -// ------------------------------------------- -Float_t AliMUONResponseV0::IntXY(AliSegmentation * segmentation) +//_____________________________________________________________________________ +Float_t +AliMUONResponseV0::GetAnod(Float_t x) const { -// Calculate charge on current pad according to Mathieson distribution -// - const Float_t kInversePitch = 1/fPitch; -// -// Integration limits defined by segmentation model -// - Float_t xi1, xi2, yi1, yi2; - segmentation->IntegrationLimits(xi1,xi2,yi1,yi2); - xi1=xi1*kInversePitch; - xi2=xi2*kInversePitch; - yi1=yi1*kInversePitch; - yi2=yi2*kInversePitch; -// -// The Mathieson function - Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1); - Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2); - - Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1); - Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2); + /// Return wire coordinate closest to x. - - return Float_t(4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))* - fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1))); + Int_t n = Int_t(x/Pitch()); + Float_t wire = (x>0) ? n+0.5 : n-0.5; + return Pitch()*wire; } -Int_t AliMUONResponseV0::DigitResponse(Int_t digit, AliMUONTransientDigit* /*where*/) +//______________________________________________________________________________ +void +AliMUONResponseV0::DisIntegrate(const AliMUONHit& hit, TList& digits) { - // add white noise and do zero-suppression and signal truncation -// Float_t meanNoise = gRandom->Gaus(1, 0.2); - // correct noise for slat chambers; - // one more field to add to AliMUONResponseV0 to allow different noises ???? - Float_t meanNoise = gRandom->Gaus(1.5, 0.2); - Float_t noise = gRandom->Gaus(0, meanNoise); - digit+=(Int_t)noise; - if ( digit <= ZeroSuppression()) digit = 0; - if ( digit > MaxAdc()) digit=MaxAdc(); - return digit; -} - - - - - + /// Go from 1 hit to a list of digits. + /// The energy deposition of that hit is first converted into charge + /// (in IntPH() method), and then this charge is dispatched on several + /// pads, according to the Mathieson distribution. + + digits.Clear(); + + Int_t detElemId = hit.DetElemId(); + + // Width of the integration area + Double_t dx = SigmaIntegration()*ChargeSpreadX(); + Double_t dy = SigmaIntegration()*ChargeSpreadY(); + + // Use that (dx,dy) to specify the area upon which + // we will iterate to spread charge into. + Double_t x,y,z; + Global2Local(detElemId,hit.X(),hit.Y(),hit.Z(),x,y,z); + x = GetAnod(x); + TVector2 hitPosition(x,y); + AliMpArea area(hitPosition,TVector2(dx,dy)); + + // Get pulse height from energy loss. + Float_t qtot = IntPH(hit.Eloss()); + + // Get the charge correlation between cathodes. + Float_t currentCorrel = TMath::Exp(gRandom->Gaus(0.0,ChargeCorrel()/2.0)); + for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++cath ) + { + Float_t qcath = qtot * ( cath == 0 ? currentCorrel : 1.0/currentCorrel); + + // Get an iterator to loop over pads, within the given area. + const AliMpVSegmentation* seg = + AliMpSegmentation::Instance() + ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cath)); + + AliMpVPadIterator* it = seg->CreateIterator(area); + + if (!it) + { + AliError(Form("Could not get iterator for detElemId %d",detElemId)); + return; + } + + // Start loop over pads. + it->First(); + + if ( it->IsDone() ) + { + // Exceptional case : iterator is built, but is invalid from the start. + AliMpPad pad = seg->PadByPosition(area.Position(),kFALSE); + if ( pad.IsValid() ) + { + AliWarning(Form("Got an invalid iterator bug (area.Position() is within " + " DE but the iterator is void) for detElemId %d cath %d", + detElemId,cath)); + } + else + { + AliError(Form("Got an invalid iterator bug for detElemId %d cath %d." + "Might be a bad hit ? area.Position()=(%e,%e) " + "Dimensions()=(%e,%e)", + detElemId,cath,area.Position().X(),area.Position().Y(), + area.Dimensions().X(),area.Dimensions().Y())); + } + delete it; + return; + } + + while ( !it->IsDone() ) + { + // For each pad given by the iterator, compute the charge of that + // pad, according to the Mathieson distribution. + AliMpPad pad = it->CurrentItem(); + TVector2 lowerLeft(hitPosition-pad.Position()-pad.Dimensions()); + TVector2 upperRight(lowerLeft + pad.Dimensions()*2.0); + Float_t qp = TMath::Abs(fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(), + upperRight.X(),upperRight.Y())); + + Int_t icharge = Int_t(qp*qcath); + + if ( qp > fChargeThreshold ) + { + // If we're above threshold, then we create a digit, + // and fill it with relevant information, including electronics. + AliMUONDigit* d = new AliMUONDigit(detElemId,pad.GetLocation().GetFirst(), + pad.GetLocation().GetSecond(),cath); + d->SetPadXY(pad.GetIndices().GetFirst(),pad.GetIndices().GetSecond()); + d->SetCharge(icharge); + digits.Add(d); + } + it->Next(); + } + delete it; + } +}