X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSsimulationSSD.cxx;h=b28c28803109989f5f936ffa553b04f95ed25b7d;hb=0854b067d05f31ae31cefa0f144c4cf62bf501d8;hp=f15ee050b7590e3c204403798ca558d6e0f96c52;hpb=14a74335211f412ee0e8c98f83faa3157a733234;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSSD.cxx b/ITS/AliITSsimulationSSD.cxx index f15ee050b75..b28c2880310 100644 --- a/ITS/AliITSsimulationSSD.cxx +++ b/ITS/AliITSsimulationSSD.cxx @@ -19,6 +19,9 @@ #include #include #include +#include + +#include #include "AliITSmodule.h" #include "AliITSMapA2.h" #include "AliITSpList.h" @@ -29,16 +32,21 @@ #include "AliITShit.h" #include "AliITSdigitSSD.h" #include "AliRun.h" +#include "AliMagF.h" +#include "AliITSgeom.h" #include "AliITSsimulationSSD.h" #include "AliITSTableSSD.h" -//#include "AliITSresponseSSD.h" +#include +#include "AliMathBase.h" +using std::endl; +using std::cout; ClassImp(AliITSsimulationSSD) //////////////////////////////////////////////////////////////////////// // // // Author: Enrico Fragiacomo // // enrico.fragiacomo@ts.infn.it // -// Last revised: march 2006 // +// Last revised: june 2008 // // // // AliITSsimulationSSD is the simulation of SSD. // //////////////////////////////////////////////////////////////////////// @@ -49,7 +57,12 @@ AliITSsimulationSSD::AliITSsimulationSSD():AliITSsimulation(), fMapA2(0), fIonE(0.0), fDifConst(), -fDriftVel(){ +fDriftVel(), +fTimeResponse(NULL), +fLorentz(kFALSE), +fTanLorAngP(0), +fTanLorAngN(0) +{ //default Constructor //Inputs: // none. @@ -65,7 +78,12 @@ AliITSsimulation(dettyp), fMapA2(0), fIonE(0.0), fDifConst(), -fDriftVel(){ +fDriftVel(), +fTimeResponse(NULL), +fLorentz(kFALSE), +fTanLorAngP(0), +fTanLorAngN(0) +{ // Constructor // Input: // AliITSDetTypeSim Pointer to the SSD dettype to be used @@ -74,6 +92,7 @@ fDriftVel(){ // Return // A standard constructed AliITSsimulationSSD class + fTimeResponse = new TF1("ftimeresponse",".5*x*exp(1.-.5*x)"); Init(); } //---------------------------------------------------------------------- @@ -87,13 +106,38 @@ void AliITSsimulationSSD::Init(){ // Return // none. AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); SetDriftVelocity(); // use default values in .h file SetIonizeE(); // use default values in .h file SetDiffConst(); // use default values in .h file fpList = new AliITSpList(2,GetNStrips()); fMapA2 = new AliITSMapA2(seg); + SetLorentzDrift(simpar->GetSSDLorentzDrift()); + if (fLorentz) SetTanLorAngle(); } + +//______________________________________________________________________ +Bool_t AliITSsimulationSSD::SetTanLorAngle() { + // This function set the Tangent of the Lorentz angles. + // output: Bool_t : kTRUE in case of success + // + + if(!fDetType) { + AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set "); + return kFALSE;} + + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); + if (!fld) AliFatal("The field is not initialized"); + Double_t bz = fld->SolenoidField(); + + fTanLorAngN = TMath::Tan( simpar->LorentzAngleElectron(bz) ); + fTanLorAngP = TMath::Tan( simpar->LorentzAngleHole(bz) ); + + return kTRUE; +} + //______________________________________________________________________ AliITSsimulationSSD& AliITSsimulationSSD::operator=( const AliITSsimulationSSD &s){ @@ -108,8 +152,13 @@ AliITSsimulationSSD& AliITSsimulationSSD::operator=( this->fDifConst[1] = s.fDifConst[1]; this->fDriftVel[0] = s.fDriftVel[0]; this->fDriftVel[1] = s.fDriftVel[1]; + this->fTimeResponse = s.fTimeResponse; + this->fLorentz = s.fLorentz; + this->fTanLorAngP = s.fTanLorAngP; + this->fTanLorAngN = s.fTanLorAngN; return *this; } +/* //______________________________________________________________________ AliITSsimulation& AliITSsimulationSSD::operator=( const AliITSsimulation &s){ @@ -121,17 +170,30 @@ AliITSsimulation& AliITSsimulationSSD::operator=( return *this; } +*/ //______________________________________________________________________ AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source): - AliITSsimulation(source){ + AliITSsimulation(source), +fMapA2(source.fMapA2), +fIonE(source.fIonE), +fDifConst(), +fDriftVel(), +fTimeResponse(source.fTimeResponse), +fLorentz(source.fLorentz), +fTanLorAngP(source.fTanLorAngP), +fTanLorAngN(source.fTanLorAngN) +{ // copy constructor - - *this = source; + fDifConst[0] = source.fDifConst[0]; + fDifConst[1] = source.fDifConst[1]; + fDriftVel[0] = source.fDriftVel[0]; + fDriftVel[1] = source.fDriftVel[1]; } //______________________________________________________________________ AliITSsimulationSSD::~AliITSsimulationSSD() { // destructor delete fMapA2; + delete fTimeResponse; //delete fDCS; } //______________________________________________________________________ @@ -208,36 +270,49 @@ void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod, Double_t x1=0.0, y1=0.0, z1=0.0; Double_t de=0.0; Int_t module = mod->GetIndex(); + Double_t tof = 0.; + AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2); - TObjArray *hits = mod->GetHits(); - Int_t nhits = hits->GetEntriesFast(); - if (nhits<=0) return; - AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips()); - module = mod->GetIndex(); - if ( mod->GetLayer() == 6 ) seg->SetLayer(6); - if ( mod->GetLayer() == 5 ) seg->SetLayer(5); - for(Int_t i=0; iGetHit(i)->GetXL() << " "<GetHit(i)->GetYL(); - cout << " " << mod->GetHit(i)->GetZL(); - cout << endl; - } // end if - if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) { - HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav); - if (lasttrack != idtrack || i==(nhits-1)) { - GetList(idtrack,i,module,pList,tav); - } // end if - lasttrack=idtrack; + TObjArray *hits = mod->GetHits(); + Int_t nhits = hits->GetEntriesFast(); + if (nhits<=0) return; + AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips()); + module = mod->GetIndex(); + if ( mod->GetLayer() == 6 ) seg->SetLayer(6); + if ( mod->GetLayer() == 5 ) seg->SetLayer(5); + + for(Int_t i=0; iGetHit(i)->GetXL() << " "<GetHit(i)->GetYL(); + cout << " " << mod->GetHit(i)->GetZL(); + cout << endl; + } // end if + if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) { + + // Scale down dE/dx according to the hit's TOF wrt to the trigger + // Necessary for pileup simulation + // EF - 21/04/09 + tof = mod->GetHit(i)->GetTOF(); + tof *= 1.E+6; // convert time in microsecond + if(tof<2.) de = de * fTimeResponse->Eval(-1.*tof+2.); + else de = 0.; + // + + HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav); + if (lasttrack != idtrack || i==(nhits-1)) { + GetList(idtrack,i,module,pList,tav); } // end if - } // end loop over hits - delete tav; tav=0; - return; + lasttrack=idtrack; + } // end if + } // end loop over hits + delete tav; tav=0; + return; } //---------------------------------------------------------------------- void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, @@ -249,8 +324,8 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2); // Turns hits in SSD module into one or more digits. - Float_t tang[2] = {0.0,0.0}; - seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo + //Float_t tang[2] = {0.0,0.0}; + //seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo Double_t x, y, z; Double_t dex=0.0, dey=0.0, dez=0.0; Double_t pairs; // pair generation energy per step. @@ -258,66 +333,151 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, Double_t tdrift[2] = {0.,0.}; // time of drift Double_t w; Double_t inf[2], sup[2], par0[2]; - + + // Set up corrections for Lorentz drift (ExB) + Double_t tanLorAngP = fTanLorAngP; + Double_t tanLorAngN = fTanLorAngN; + if(seg->GetLayer()==6) { + tanLorAngP = -1.*fTanLorAngP; + tanLorAngN = -1.*fTanLorAngN; + } + // Steps in the module are determined "manually" (i.e. No Geant) // NumOfSteps divide path between entering and exiting hits in steps Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez); // Enery loss is equally distributed among steps de = de/numOfSteps; pairs = de/GetIonizeE(); // e-h pairs generated + + //----------------------------------------------------- + // stepping + //----------------------------------------------------- for(Int_t j=0; j (seg->Dy()/2+10)*1.0E-4 ) { // check if particle is within the detector Warning("HitToDigit", - "hit out of detector y0=%e,y=%e,dey=%e,j =%e module=%d", - y0,y,dey,j,module); + "hit out of detector y0=%e,y=%e,dey=%e,j =%d module=%d, exceed=%e", + y0,y,dey,j,module, y-(seg->Dy()/2+10)*1.0E-4); return; } // end if z = z0 + (j+0.5)*dez; + if(GetDebug(4)) cout <<"HitToDigit "<GetLayer()==6) { + y=-y; // Lay6 module has sensor up-side-down!!! + } + + Int_t k; + //--------------------------------------------------------- + // Pside + //------------------------------------------------------------ + k=0; + + // w is the coord. perpendicular to the strips + // Float_t xp=x*1.e+4,zp=z*1.e+4; // microns + Float_t xp=x,zp=z; + + // correction for the Lorentz's angle + if(fLorentz) { + Float_t deltaxp = (y+(seg->Dy()*1.0E-4)/2)*tanLorAngP; + xp+=deltaxp; + } + + seg->GetPadTxz(xp,zp); + // calculate drift time // y is the minimum path tdrift[0] = (y+(seg->Dy()*1.0E-4)/2)/GetDriftVelocity(0); + + w = xp; // P side strip number + + if((w<(-0.5)) || (w>(GetNStrips()-0.5))) { + // this check rejects hits in regions not covered by strips + // 0.5 takes into account boundaries + if(GetDebug(4)) cout << "Dead SSD region, x,z="<Dy()*1.0E-4)/2-y)*tanLorAngN; + xp+=deltaxn; + } + + + seg->GetPadTxz(xp,zp); + tdrift[1] = ((seg->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1); - for(Int_t k=0; k<2; k++) { // both sides remember: 0=Pside 1=Nside - - tang[k]=TMath::Tan(tang[k]); - - // w is the coord. perpendicular to the strips - Float_t xp=x*1.e+4,zp=z*1.e+4; // microns - seg->GetPadTxz(xp,zp); - if(k==0) w = xp; // P side strip number - else w = zp; // N side strip number - - if((w<(-0.5)) || (w>(GetNStrips()-0.5))) { - // this check rejects hits in regions not covered by strips - // 0.5 takes into account boundaries - if(GetDebug(4)) cout << "x,z="<Gaus(0,res->GetNoiseP().At(ix)); + noise = (Double_t) gRandom->Gaus(0,res->GetNoiseP(ix)); // need to calibrate noise // NOTE. noise from the calibration database comes uncalibrated, // it needs to be calibrated in order to be added // to the signal. It will be decalibrated later on together with the noise - noise *= res->GetGainP(ix); + noise *= (Double_t) res->GetGainP(ix); // noise comes in ADC channels from the calibration database // It needs to be converted back to electronVolts - noise /= res->GetDEvToADC(1.); + noise /= res->GetSSDDEvToADC(1.); // Finally, noise is added to the signal signal = noise + fMapA2->GetSignal(0,ix);//get signal from map @@ -352,9 +512,9 @@ void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){ // Nside for(ix=0;ixGaus(0,res->GetNoiseN().At(ix));// give noise to signal - noise *= res->GetGainN(ix); - noise /= res->GetDEvToADC(1.); + noise = (Double_t) gRandom->Gaus(0,res->GetNoiseN(ix));// give noise to signal + noise *= (Double_t) res->GetGainN(ix); + noise /= res->GetSSDDEvToADC(1.); signal = noise + fMapA2->GetSignal(1,ix);//get signal from map fMapA2->SetHit(1,ix,signal); // give back signal to map if(signal>0.0) pList->AddNoise(1,ix,module,noise); @@ -366,37 +526,38 @@ void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) { // Apply the effect of electronic coupling between channels Int_t ix; Double_t signal=0; - AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module); - + //AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module); + AliITSSimuParam* res = fDetType->GetSimuParam(); + Double_t *contrLeft = new Double_t[GetNStrips()]; Double_t *contrRight = new Double_t[GetNStrips()]; // P side coupling for(ix=0;ix0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetCouplingPL(); + if(ix>0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetSSDCouplingPL(); else contrLeft[ix] = 0.0; - if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetCouplingPR(); + if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetSSDCouplingPR(); else contrRight[ix] = 0.0; } // loop over strips for(ix=0;ixGetCouplingPL() * fMapA2->GetSignal(0,ix) - - res->GetCouplingPR() * fMapA2->GetSignal(0,ix); + signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingPL() * fMapA2->GetSignal(0,ix) + - res->GetSSDCouplingPR() * fMapA2->GetSignal(0,ix); fMapA2->AddSignal(0,ix,signal); if(signal>0.0) pList->AddNoise(0,ix,module,signal); } // loop over strips // N side coupling for(ix=0;ix0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetCouplingNL(); + if(ix>0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetSSDCouplingNL(); else contrLeft[ix] = 0.0; - if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetCouplingNR(); + if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetSSDCouplingNR(); else contrRight[ix] = 0.0; } // loop over strips for(ix=0;ixGetCouplingNL() * fMapA2->GetSignal(0,ix) - - res->GetCouplingNR() * fMapA2->GetSignal(0,ix); + signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingNL() * fMapA2->GetSignal(0,ix) + - res->GetSSDCouplingNR() * fMapA2->GetSignal(0,ix); fMapA2->AddSignal(1,ix,signal); if(signal>0.0) pList->AddNoise(1,ix,module,signal); } // loop over strips @@ -410,20 +571,14 @@ void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) { void AliITSsimulationSSD::ApplyDeadChannels(Int_t module) { // Kill dead channels setting gain to zero - Int_t deadentries; - AliITSCalibrationSSD* res = (AliITSCalibrationSSD*)GetCalibrationModel(module); - deadentries = res->GetDeadPChannelsList().GetSize(); - //cout<AddGainP(res->GetDeadPChannelsList().At(i),0.0); - } + for(Int_t i=0;iGetDeadNChannelsList().GetSize(); - for(Int_t i=0; iAddGainN(res->GetDeadNChannelsList().At(i),0.0); - } + if(res->IsPChannelBad(i)) res->SetGainP(i,0.0); + if(res->IsNChannelBad(i)) res->SetGainN(i,0.0); + + } // loop over strips } @@ -434,7 +589,7 @@ Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) { Float_t sigm2 = sqrt2*s; Float_t integral; - integral = 0.5 * TMath::Erf( (x - av) / sigm2); + integral = 0.5 * AliMathBase::ErfFast( (x - av) / sigm2); return integral; } //______________________________________________________________________ @@ -564,7 +719,7 @@ void AliITSsimulationSSD::GetList(Int_t label,Int_t hit,Int_t mod, tav->Clear(); } //---------------------------------------------------------------------- -void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) { +void AliITSsimulationSSD::ChargeToSignal(Int_t module,const AliITSpList *pList) { // charge to signal static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS"); Float_t threshold = 0.; @@ -576,6 +731,7 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) { Float_t charges[3] = {0.0,0.0,0.0}; Float_t signal; AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); for(Int_t k=0;k<2;k++){ // both sides (0=Pside, 1=Nside) for(Int_t ix=0;ix gain=0 if( ((k==0)&&(res->GetGainP(ix)==0)) || ((k==1)&&(res->GetGainN(ix)==0))) continue; + signal = fMapA2->GetSignal(k,ix); // signal has to be uncalibrated // In real life, gains are supposed to be calculated from calibration runs, // stored in the calibration DB and used in the reconstruction @@ -591,15 +748,18 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) { else signal /= res->GetGainN(ix); // signal is converted in unit of ADC - signal = res->GetDEvToADC(fMapA2->GetSignal(k,ix)); - if(signal>4096.) signal = 4096.;//if exceeding, accumulate last one + signal = res->GetSSDDEvToADC(signal); + if(signal>4095.) signal = 4095.;//if exceeding, accumulate last one // threshold for zero suppression is set on the basis of the noise // A good value is 3*sigma_noise - if(k==0) threshold = res->GetNoiseP().At(ix); - else threshold = res->GetNoiseN().At(ix); - threshold *= res->GetZSThreshold(); // threshold at 3 sigma noise + if(k==0) threshold = res->GetNoiseP(ix); + else threshold = res->GetNoiseN(ix); + + threshold *= simpar->GetSSDZSThreshold(); // threshold at 3 sigma noise + if(signal < threshold) continue; + //cout<>(istream &os,AliITSsimulationSSD &source){ - -