X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONResponseV0.cxx;h=5bc82aea728017280d5d068f973aa26f4de1f086;hb=b7342e4d71169a9f54c9c0fad7c6ecb0c5893e0e;hp=4a7f347c404aba1d098d486827b0d2c5e17dcf18;hpb=a9e2aefa97f1153d6f61e580a32d396156706b7b;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONResponseV0.cxx b/MUON/AliMUONResponseV0.cxx index 4a7f347c404..5bc82aea728 100644 --- a/MUON/AliMUONResponseV0.cxx +++ b/MUON/AliMUONResponseV0.cxx @@ -13,77 +13,346 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.1.2.1 2000/06/09 21:33:35 morsch -AliMUONResponse code from AliMUONSegResV0.cxx +/* $Id$ */ -*/ +//----------------------------------------------------------------------------- +// Class AliMUONResponseV0 +// -------------------------- +// Implementation of +// Mathieson response +//----------------------------------------------------------------------------- #include "AliMUONResponseV0.h" -#include "AliMUONSegmentation.h" +#include "AliMUON.h" +#include "AliMUONConstants.h" +#include "AliMUONDigit.h" +#include "AliMUONGeometryTransformer.h" +#include "AliMUONHit.h" +#include "AliMUONConstants.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 +using std::endl; +using std::cout; +/// \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), + fIsTailEffect(kFALSE) +{ + /// 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), +fIsTailEffect(kFALSE) +{ + /// copy ctor + other.CopyTo(*this); +} + +//__________________________________________________________________________ +AliMUONResponseV0& +AliMUONResponseV0::operator=(const AliMUONResponseV0& other) +{ + /// Assignment operator + other.CopyTo(*this); + return *this; +} -ClassImp(AliMUONResponseV0) -Float_t AliMUONResponseV0::IntPH(Float_t eloss) +//__________________________________________________________________________ +AliMUONResponseV0::~AliMUONResponseV0() { - // Calculate charge from given ionization energy loss +/// 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 + 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 + fMathieson->SetSqrtKy3AndDeriveKy2Ky4(SqrtKy3); +} + //__________________________________________________________________________ +Float_t AliMUONResponseV0::IntPH(Float_t eloss) const +{ + /// 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++) { - charge -= fChargeSlope*TMath::Log(gRandom->Rndm()); + Float_t arg=0.; + while(!arg) arg = gRandom->Rndm(); + charge -= fChargeSlope*TMath::Log(arg); } return charge; } -// ------------------------------------------- -Float_t AliMUONResponseV0::IntXY(AliMUONSegmentation * 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) +//______________________________________________________________________________ +void +AliMUONResponseV0::DisIntegrate(const AliMUONHit& hit, TList& digits, Float_t timeDif) { - // add white noise and do zero-suppression and signal truncation - Float_t meanNoise = gRandom->Gaus(1, 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(); + Double_t hitX = hit.X() ; + Double_t hitY = hit.Y() ; + Double_t hitZ = hit.Z() ; + // Width of the integration area + Double_t dx = SigmaIntegration()*ChargeSpreadX(); + Double_t dy = SigmaIntegration()*ChargeSpreadY(); + + //Modify to take the tailing effect. + if(fIsTailEffect){ + Double_t locX,locY,locZ,globXCenter,globYCenter,globZ; + Int_t para = 5; // This parameter is a natural number(excluding zero), higher the value less is the tailing effect + Double_t termA = 1.0; + Double_t termB = 1.0; + if(para>0){ + for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++cath ) + { + // Get an iterator to loop over pads, within the given area. + const AliMpVSegmentation* seg = + AliMpSegmentation::Instance() + ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cath)); + AliMp::PlaneType plane = seg->PlaneType(); + + if(plane == AliMp::kBendingPlane) { + Global2Local(detElemId,hitX,hitY,hitZ,locX,locY,locZ); + AliMpPad pad = seg->PadByPosition(locX,locY,kFALSE); + if(pad.IsValid()){ + Double_t locYCenter = pad.GetPositionY(); + Double_t locXCenter = pad.GetPositionX(); + const AliMUONGeometryTransformer* transformer = muon()->GetGeometryTransformer(); + transformer->Local2Global(detElemId,locXCenter,locYCenter,locZ,globXCenter,globYCenter,globZ); + for(Int_t itime = 0; itime 0 condn + }// if tail effect + // Use that (dx,dy) to specify the area upon which + // we will iterate to spread charge into. + Double_t x,y,z; + Global2Local(detElemId,hitX,hitY,hitZ,x,y,z); + x = GetAnod(x); + AliMpArea area(x,y,dx,dy); + + // Get pulse height from energy loss. + Float_t qtot = IntPH(hit.Eloss()); + + // If from a pileup event we apply a reduction factor to the charge + if (timeDif!=0){ + qtot = AliMUONConstants::ReducedQTot(qtot,timeDif); + } + + // Scale the charge to it'll (roughly) be in fC + qtot *= AliMUONConstants::DefaultADC2MV()*AliMUONConstants::DefaultA0()*AliMUONConstants::DefaultCapa(); + + // 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.GetPositionX(),area.GetPositionY(), + kFALSE); + if ( pad.IsValid() ) + { + AliDebug(1, Form("Got an invalid iterator bug (area.Position() is within " + " DE but the iterator is void) for detElemId %d cath %d", + detElemId,cath)); + } + else + { + AliDebug(1, 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.GetPositionX(),area.GetPositionY(), + area.GetDimensionX(),area.GetDimensionY())); + } + 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(TVector2(x,y)-TVector2(pad.GetPositionX(),pad.GetPositionY())- + TVector2(pad.GetDimensionX(),pad.GetDimensionY())); + TVector2 upperRight(lowerLeft + TVector2(pad.GetDimensionX(),pad.GetDimensionY())*2.0); + Float_t qp = TMath::Abs(fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(), + upperRight.X(),upperRight.Y())); + + if ( qp > fChargeThreshold && + qp*qcath > AliMUONConstants::DefaultADC2MV()*AliMUONConstants::DefaultA0()*AliMUONConstants::DefaultCapa() ) + { + // If we're above threshold, then we create a digit, + // and fill it with relevant information, including electronics. + + // note that the second condition above is to be backward compatible (when + // the sdigitizer was making a cut on Int_t(qp*qcath) > 0 and qcath was in ADC, not in fC) + + AliMUONDigit* d = new AliMUONDigit(detElemId,pad.GetManuId(), + pad.GetManuChannel(),cath); + d->SetPadXY(pad.GetIx(),pad.GetIy()); + d->SetCharge(qp*qcath); + digits.Add(d); + } + it->Next(); + } + delete it; + } +}