X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=ITS%2FAliITSsimulationSPD.cxx;h=ac65a13e6345d7eded82390dbd4553d5ab3e9b6e;hb=bf27d64a0f4d49bc005c0fc4ad812e1f92578c67;hp=20b2c930efe5a5307a43d8790880cb89533e9db0;hpb=5a200b64ec3030b25e8ca6bb4f75ea92db634eb4;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSPD.cxx b/ITS/AliITSsimulationSPD.cxx index 20b2c930efe..ac65a13e634 100644 --- a/ITS/AliITSsimulationSPD.cxx +++ b/ITS/AliITSsimulationSPD.cxx @@ -1,607 +1,1030 @@ -#include -#include +/************************************************************************** +* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +* * +* Author: The ALICE Off-line Project. * +* Contributors are mentioned in the code where appropriate. * +* * +* Permission to use, copy, modify and distribute this software and its * +* documentation strictly for non-commercial purposes is hereby granted * +* without fee, provided that the above copyright notice appears in all * +* copies and that both the copyright notice and this permission notice * +* appear in the supporting documentation. The authors make no claims * +* about the suitability of this software for any purpose. It is * +* provided "as is" without express or implied warranty. * +**************************************************************************/ + +/* +$Id$ +*/ + +#include +#include #include -#include #include -#include - - -#include "AliRun.h" #include "AliITS.h" +#include "AliITSdigitSPD.h" #include "AliITShit.h" -#include "AliITSdigit.h" #include "AliITSmodule.h" -#include "AliITSMapA2.h" +#include "AliITSpList.h" +#include "AliITSCalibrationSPD.h" +#include "AliITSsegmentationSPD.h" #include "AliITSsimulationSPD.h" -#include "AliITSsegmentation.h" -#include "AliITSresponse.h" - - +#include "AliLog.h" +#include "AliRun.h" +#include "AliMagF.h" +#include "AliMathBase.h" +//#define DEBUG +using std::endl; +using std::cout; ClassImp(AliITSsimulationSPD) //////////////////////////////////////////////////////////////////////// -// Version: 0 -// Written by Boris Batyunya -// December 20 1999 +// Version: 1 +// Modified by D. Elia, G.E. Bruno, H. Tydesjo +// Fast diffusion code by Bjorn S. Nilsen +// March-April 2006 +// October 2007: GetCalibrationObjects() removed // -// AliITSsimulationSPD is the simulation of SPDs -//________________________________________________________________________ - +// Version: 0 +// Written by Boris Batyunya +// December 20 1999 +// +// +// AliITSsimulationSPD is to do the simulation of SPDs. +// +//////////////////////////////////////////////////////////////////////// -AliITSsimulationSPD::AliITSsimulationSPD() -{ - // constructor - fResponse = 0; - fSegmentation = 0; - fMapA2=0; - fHis = 0; - fNoise=0.; - fBaseline=0.; - fNPixelsZ=0; - fNPixelsX=0; +//______________________________________________________________________ +AliITSsimulationSPD::AliITSsimulationSPD(): +AliITSsimulation(), +fHis(0), +fSPDname(), +fCoupling(), +fLorentz(kFALSE), +fTanLorAng(0), +fStrobe(kTRUE), +fStrobeLenght(4), +fStrobePhase(-12.5e-9){ + // Default constructor. + // Inputs: + // none. + // Outputs: + // none. + // Return: + // A default constructed AliITSsimulationSPD class. + + AliDebug(1,Form("Calling default constructor")); +// Init(); } - - -//_____________________________________________________________________________ - -AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) { - // standard constructor - - fHis = 0; - fResponse = resp; - fSegmentation = seg; - - fResponse->GetNoiseParam(fNoise,fBaseline); - - fMapA2 = new AliITSMapA2(fSegmentation); - - // - - fNPixelsZ=fSegmentation->Npz(); - fNPixelsX=fSegmentation->Npx(); - +//______________________________________________________________________ +AliITSsimulationSPD::AliITSsimulationSPD(AliITSDetTypeSim *dettyp): +AliITSsimulation(dettyp), +fHis(0), +fSPDname(), +fCoupling(), +fLorentz(kFALSE), +fTanLorAng(0), +fStrobe(kTRUE), +fStrobeLenght(4), +fStrobePhase(-12.5e-9){ + // standard constructor + // Inputs: + // AliITSsegmentation *seg A pointer to the segmentation class + // to be used for this simulation + // AliITSCalibration *resp A pointer to the responce class to + // be used for this simulation + // Outputs: + // none. + // Return: + // A default constructed AliITSsimulationSPD class. + + AliDebug(1,Form("Calling standard constructor ")); + Init(); } - -//_____________________________________________________________________________ - -AliITSsimulationSPD::~AliITSsimulationSPD() { - // destructor - - delete fMapA2; - - if (fHis) { - fHis->Delete(); - delete fHis; - } +//______________________________________________________________________ +void AliITSsimulationSPD::Init(){ + // Initilization + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + const Double_t kmictocm = 1.0e-4; // convert microns to cm. + + SetModuleNumber(0); + SetEventNumber(0); + SetMap(new AliITSpList(GetNPixelsZ(),GetNPixelsX())); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + Double_t bias = simpar->GetSPDBiasVoltage(); +// cout << "Bias Voltage --> " << bias << endl; // dom + simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),bias); +// set kind of coupling ("old" or "new") + char opt[20]; + simpar->GetSPDCouplingOption(opt); + char *old = strstr(opt,"old"); + if (old) { + fCoupling=2; + } else { + fCoupling=1; + } // end if + SetLorentzDrift(simpar->GetSPDLorentzDrift()); + if (fLorentz) SetTanLorAngle(simpar->GetSPDLorentzHoleWeight()); + //SetStrobeGeneration(kFALSE); + if (fStrobe) GenerateStrobePhase(); } - - -//__________________________________________________________________________ -AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){ - // Copy Constructor - if(&source == this) return; - this->fMapA2 = source.fMapA2; - this->fNoise = source.fNoise; - this->fBaseline = source.fBaseline; - this->fNPixelsX = source.fNPixelsX; - this->fNPixelsZ = source.fNPixelsZ; - this->fHis = source.fHis; - return; +//______________________________________________________________________ +Bool_t AliITSsimulationSPD::SetTanLorAngle(Double_t WeightHole) { + // This function set the Tangent of the Lorentz angle. + // A weighted average is used for electrons and holes + // Input: Double_t WeightHole: wheight for hole: it should be in the range [0,1] + // output: Bool_t : kTRUE in case of success + // + if(!fDetType) { + AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set "); + return kFALSE;} + if(WeightHole<0) { + WeightHole=0.; + AliWarning("AliITSsimulationSPD::SetTanLorAngle: You have asked for negative Hole weight"); + AliWarning("AliITSsimulationSPD::SetTanLorAngle: I'm going to use only electrons"); + } + if(WeightHole>1) { + WeightHole=1.; + AliWarning("AliITSsimulationSPD::SetTanLorAngle: You have asked for weight > 1"); + AliWarning("AliITSsimulationSPD::SetTanLorAngle: I'm going to use only holes"); + } + Double_t WeightEle=1.-WeightHole; + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); + if (!fld) AliFatal("The field is not initialized"); + Double_t bz = fld->SolenoidField(); + fTanLorAng = TMath::Tan(WeightHole*simpar->LorentzAngleHole(bz) + + WeightEle*simpar->LorentzAngleElectron(bz)); + return kTRUE; } +//______________________________________________________________________ +AliITSsimulationSPD::~AliITSsimulationSPD(){ + // destructor + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + + if (fHis) { + fHis->Delete(); + delete fHis; + } // end if fHis +} +//______________________________________________________________________ +AliITSsimulationSPD::AliITSsimulationSPD(const + AliITSsimulationSPD + &s) : AliITSsimulation(s), +fHis(s.fHis), +fSPDname(s.fSPDname), +fCoupling(s.fCoupling), +fLorentz(s.fLorentz), +fTanLorAng(s.fTanLorAng), +fStrobe(s.fStrobe), +fStrobeLenght(s.fStrobeLenght), +fStrobePhase(s.fStrobePhase){ + // Copy Constructor + // Inputs: + // AliITSsimulationSPD &s The original class for which + // this class is a copy of + // Outputs: + // none. + // Return: -//_________________________________________________________________________ -AliITSsimulationSPD& - AliITSsimulationSPD::operator=(const AliITSsimulationSPD &source) { - // Assignment operator - if(&source == this) return *this; - this->fMapA2 = source.fMapA2; - this->fNoise = source.fNoise; - this->fBaseline = source.fBaseline; - this->fNPixelsX = source.fNPixelsX; - this->fNPixelsZ = source.fNPixelsZ; - this->fHis = source.fHis; - return *this; - } -//_____________________________________________________________________________ - -void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy) -{ - // digitize module - - const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons - // for 3.6 eV/pair - const Float_t kconv = 10000.; // cm -> microns - - Float_t spdLength = fSegmentation->Dz(); - Float_t spdWidth = fSegmentation->Dx(); - - Float_t difCoef, dum; - fResponse->DiffCoeff(difCoef,dum); - - Float_t zPix0 = 1e+6; - Float_t xPix0 = 1e+6; - Float_t yPrev = 1e+6; - - Float_t zPitch = fSegmentation->Dpz(0); - Float_t xPitch = fSegmentation->Dpx(0); - - TObjArray *fHits = mod->GetHits(); - Int_t nhits = fHits->GetEntriesFast(); - if (!nhits) return; - - // Array of pointers to the label-signal list - - Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;; - Float_t **pList = new Float_t* [maxNDigits]; - memset(pList,0,sizeof(Float_t*)*maxNDigits); - Int_t indexRange[4] = {0,0,0,0}; - - // Fill detector maps with GEANT hits - // loop over hits in the module - static Bool_t first; - Int_t lasttrack=-2; - Int_t hit, iZi, jz, jx; - for (hit=0;hitAt(hit); - Int_t layer = iHit->GetLayer(); - Float_t yPix0 = -73; - if(layer == 1) yPix0 = -77; - - // work with the idtrack=entry number in the TreeH - Int_t idhit,idtrack; - mod->GetHitTrackAndHitIndex(hit,idtrack,idhit); - //Int_t idtrack=mod->GetHitTrackIndex(hit); - // or store straight away the particle position in the array - // of particles : - - //b.b. if(iHit->StatusEntering()) idhit=hit; - Int_t itrack = iHit->GetTrack(); - Int_t dray = 0; - - if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE; - - // Int_t parent = iHit->GetParticle()->GetFirstMother(); - Int_t partcode = iHit->GetParticle()->GetPdgCode(); - -// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+ -// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda - - Float_t px = iHit->GetPXL(); - Float_t py = iHit->GetPYL(); - Float_t pz = iHit->GetPZL(); - Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz); - - - if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- - // at p < 6 MeV/c - - - // Get hit z and x(r*phi) cordinates for each module (detector) - // in local system. - - Float_t zPix = kconv*iHit->GetZL(); - Float_t xPix = kconv*iHit->GetXL(); - Float_t yPix = kconv*iHit->GetYL(); - - // Get track status - Int_t status = iHit->GetTrackStatus(); - - // Check boundaries - if(zPix > spdLength/2) zPix = spdLength/2 - 10; - if(zPix < 0 && zPix < -spdLength/2) zPix = -spdLength/2 + 10; - if(xPix > spdWidth/2) xPix = spdWidth/2 - 10; - if(xPix < 0 && xPix < -spdWidth/2) zPix = -spdWidth/2 + 10; - - Int_t trdown = 0; - - // enter Si or after event in Si - if (status == 66 ) { - zPix0 = zPix; - xPix0 = xPix; - yPrev = yPix; - } - - Float_t depEnergy = iHit->GetIonization(); - // skip if the input point to Si - - if(depEnergy <= 0.) continue; - - // if track returns to the opposite direction: - if (yPix < yPrev) { - trdown = 1; - } - - - // take into account the holes diffusion inside the Silicon - // the straight line between the entrance and exit points in Si is - // divided into the several steps; the diffusion is considered - // for each end point of step and charge - // is distributed between the pixels through the diffusion. - - - // ---------- the diffusion in Z (beam) direction ------- - - Float_t charge = depEnergy*kEnToEl; // charge in e- - Float_t drPath = 0.; - Float_t tang = 0.; - Float_t sigmaDif = 0.; - Float_t zdif = zPix - zPix0; - Float_t xdif = xPix - xPix0; - Float_t ydif = TMath::Abs(yPix - yPrev); - Float_t ydif0 = TMath::Abs(yPrev - yPix0); - - if(ydif < 1) continue; // ydif is not zero - - Float_t projDif = sqrt(xdif*xdif + zdif*zdif); - Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1; - Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1; - - // number of the steps along the track: - Int_t nsteps = ndZ; - if(ndX > ndZ) nsteps = ndX; - if(nsteps < 6) nsteps = 6; // minimum number of the steps - - if (projDif < 5 ) { - drPath = (yPix-yPix0)*1.e-4; - drPath = TMath::Abs(drPath); // drift path in cm - sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm - sigmaDif = sigmaDif*kconv; // sigma diffusion in microns - nsteps = 1; - } - - if(projDif > 5) tang = ydif/projDif; - Float_t dCharge = charge/nsteps; // charge in e- for one step - Float_t dZ = zdif/nsteps; - Float_t dX = xdif/nsteps; - - for (iZi = 1;iZi <= nsteps;iZi++) { - Float_t dZn = iZi*dZ; - Float_t dXn = iZi*dX; - Float_t zPixn = zPix0 + dZn; - Float_t xPixn = xPix0 + dXn; - - if(projDif >= 5) { - Float_t dProjn = sqrt(dZn*dZn+dXn*dXn); - drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm - if(trdown == 0) { - drPath = TMath::Abs(drPath) + ydif0*1.e-4; - } - if(trdown == 1) { - drPath = ydif0*1.e-4 - TMath::Abs(drPath); - drPath = TMath::Abs(drPath); - } - sigmaDif = difCoef*sqrt(drPath); - sigmaDif = sigmaDif*kconv; // sigma diffusion in microns - } - - zPixn = (zPixn + spdLength/2.); - xPixn = (xPixn + spdWidth/2.); - Int_t nZpix, nXpix; - fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix); - zPitch = fSegmentation->Dpz(nZpix); - fSegmentation->GetPadTxz(xPixn,zPixn); - // set the window for the integration - Int_t jzmin = 1; - Int_t jzmax = 3; - if(nZpix == 1) jzmin =2; - if(nZpix == fNPixelsZ) jzmax = 2; - - Int_t jxmin = 1; - Int_t jxmax = 3; - if(nXpix == 1) jxmin =2; - if(nXpix == fNPixelsX) jxmax = 2; - - Float_t zpix = nZpix; - Float_t dZright = zPitch*(zpix - zPixn); - Float_t dZleft = zPitch - dZright; - - Float_t xpix = nXpix; - Float_t dXright = xPitch*(xpix - xPixn); - Float_t dXleft = xPitch - dXright; - - Float_t dZprev = 0.; - Float_t dZnext = 0.; - Float_t dXprev = 0.; - Float_t dXnext = 0.; - - for(jz=jzmin; jz <=jzmax; jz++) { - if(jz == 1) { - dZprev = -zPitch - dZleft; - dZnext = -dZleft; - } - if(jz == 2) { - dZprev = -dZleft; - dZnext = dZright; - } - if(jz == 3) { - dZprev = dZright; - dZnext = dZright + zPitch; - } - // kz changes from 1 to the fNofPixels(270) - Int_t kz = nZpix + jz -2; - - Float_t zArg1 = dZprev/sigmaDif; - Float_t zArg2 = dZnext/sigmaDif; - Float_t zProb1 = TMath::Erfc(zArg1); - Float_t zProb2 = TMath::Erfc(zArg2); - Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge; - - - // ----------- holes diffusion in X(r*phi) direction -------- - - if(dZCharge > 1.) { - for(jx=jxmin; jx <=jxmax; jx++) { - if(jx == 1) { - dXprev = -xPitch - dXleft; - dXnext = -dXleft; - } - if(jx == 2) { - dXprev = -dXleft; - dXnext = dXright; - } - if(jx == 3) { - dXprev = dXright; - dXnext = dXright + xPitch; - } - Int_t kx = nXpix + jx -2; - - Float_t xArg1 = dXprev/sigmaDif; - Float_t xArg2 = dXnext/sigmaDif; - Float_t xProb1 = TMath::Erfc(xArg1); - Float_t xProb2 = TMath::Erfc(xArg2); - Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge; - - if(dXCharge > 1.) { - Int_t index = kz-1; - - if (first) { - indexRange[0]=indexRange[1]=index; - indexRange[2]=indexRange[3]=kx-1; - first=kFALSE; - } - - indexRange[0]=TMath::Min(indexRange[0],kz-1); - indexRange[1]=TMath::Max(indexRange[1],kz-1); - indexRange[2]=TMath::Min(indexRange[2],kx-1); - indexRange[3]=TMath::Max(indexRange[3],kx-1); - - // build the list of digits for this module - Double_t signal=fMapA2->GetSignal(index,kx-1); - signal+=dXCharge; - fMapA2->SetHit(index,kx-1,(double)signal); - } // dXCharge > 1 e- - } // jx loop - } // dZCharge > 1 e- - } // jz loop - } // iZi loop - - if (status == 65) { // the step is inside of Si - zPix0 = zPix; - xPix0 = xPix; - } - yPrev = yPix; - - if(dray == 0) { - GetList(itrack,idhit,pList,indexRange); - } - - lasttrack=itrack; - } // hit loop inside the module - - - // introduce the electronics effects and do zero-suppression - ChargeToSignal(pList); - - // clean memory - - fMapA2->ClearMap(); - - -} - -//--------------------------------------------- -void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange) -{ - // lop over nonzero digits - +} +//______________________________________________________________________ +AliITSsimulationSPD& AliITSsimulationSPD::operator=(const + AliITSsimulationSPD &s){ + // Assignment operator + // Inputs: + // AliITSsimulationSPD &s The original class for which + // this class is a copy of + // Outputs: + // none. + // Return: + + if(&s == this) return *this; + this->fHis = s.fHis; + fCoupling = s.fCoupling; + fSPDname = s.fSPDname; + fLorentz = s.fLorentz; + fTanLorAng = s.fTanLorAng; + fStrobe = s.fStrobe; + fStrobeLenght = s.fStrobeLenght; + fStrobePhase = s.fStrobePhase; + return *this; +} +/* +//______________________________________________________________________ +AliITSsimulation& AliITSsimulationSPD::operator=(const + AliITSsimulation &s){ + // Assignment operator + // Inputs: + // AliITSsimulationSPD &s The original class for which + // this class is a copy of + // Outputs: + // none. + // Return: + + if(&s == this) return *this; + Error("AliITSsimulationSPD","Not allowed to make a = with " + "AliITSsimulationSPD","Using default creater instead"); + + return *this; +} +*/ +//______________________________________________________________________ +void AliITSsimulationSPD::InitSimulationModule(Int_t module, Int_t event){ + // This function creates maps to build the list of tracks for each + // summable digit. Inputs defined by base class. + // Inputs: + // Int_t module // Module number to be simulated + // Int_t event // Event number to be simulated + // Outputs: + // none + // Returns: + // none + + AliDebug(1,Form("(module=%d,event=%d)",module,event)); + const Double_t kmictocm = 1.0e-4; // convert microns to cm. + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + SetModuleNumber(module); + SetEventNumber(event); + simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),simpar->GetSPDBiasVoltage(module)); + ClearMap(); +} +//_____________________________________________________________________ +void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod,Int_t, + Int_t event){ + // This function begins the work of creating S-Digits. Inputs defined + // by base class. + // Inputs: + // AliITSmodule *mod // module + // Int_t // not used + // Int_t event // Event number + // Outputs: + // none + // Return: + // test // test returns kTRUE if the module contained hits + // // test returns kFALSE if it did not contain hits + + AliDebug(1,Form("(mod=%p, ,event=%d)",mod,event)); + if(!(mod->GetNhits())){ + AliDebug(1,Form("In event %d module %d there are %d hits returning.", + event, mod->GetIndex(),mod->GetNhits())); + return;// if module has no hits don't create Sdigits + } // end if + SetModuleNumber(mod->GetIndex()); + if (fStrobe) if(event != GetEventNumber()) GenerateStrobePhase(); + SetEventNumber(event); + InitSimulationModule( GetModuleNumber() , event ); + // HitToSDigit(mod); + HitToSDigitFast(mod); + if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag()) AddNoisyPixels(); + if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels(); + +// cout << "After Remove in SDigitiseModule !!!!!" << endl; // dom +// cout << "Module " << mod->GetIndex() << " Event " << event << endl; // dom + WriteSDigits(); + ClearMap(); +} +//______________________________________________________________________ +void AliITSsimulationSPD::WriteSDigits(){ + // This function adds each S-Digit to pList + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none + Int_t ix, nix, iz, niz; + static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS"); + + AliDebug(1,Form("Writing SDigits for module %d",GetModuleNumber())); +// cout << "WriteSDigits for module " << GetModuleNumber() << endl; // dom + GetMap()->GetMaxMapIndex(niz, nix); + for(iz=0; izGetSignalOnly(iz,ix)>0.0){ +// cout << " Signal gt 0 iz ix " << iz << ix << " Module " << GetModuleNumber() << endl; // dom + aliITS->AddSumDigit(*(GetMap()->GetpListItem(iz,ix))); + if(AliDebugLevel()>0) { + AliDebug(1,Form("%d, %d",iz,ix)); + cout << *(GetMap()->GetpListItem(iz,ix)) << endl; + } // end if GetDebug + } // end if GetMap()->GetSignalOnly(iz,ix)>0.0 + } // end for iz,ix + return; +} +//______________________________________________________________________ +void AliITSsimulationSPD::FinishSDigitiseModule(){ + // This function calls SDigitsToDigits which creates Digits from SDigits + // Inputs: + // none + // Outputs: + // none + // Return + // none + + AliDebug(1,"()"); +// cout << "FinishSDigitiseModule for module " << GetModuleNumber() << endl; // dom + FrompListToDigits(); // Charge To Signal both adds noise and + ClearMap(); + return; +} +//______________________________________________________________________ +void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod,Int_t, + Int_t event){ + // This function creates Digits straight from the hits and then adds + // electronic noise to the digits before adding them to pList + // Each of the input variables is passed along to HitToSDigit + // Inputs: + // AliITSmodule *mod module + // Int_t Dummy. + // Int_t Dummy + // Outputs: + // none. + // Return: + // none. + + if (fStrobe) if(event != GetEventNumber()) GenerateStrobePhase(); + AliDebug(1,Form("(mod=%p,,0)",mod)); + // HitToSDigit(mod); + InitSimulationModule( mod->GetIndex(), event ); + HitToSDigitFast(mod); - //set protection - for(int k=0;k<4;k++) { - if (indexRange[k] < 0) indexRange[k]=0; + if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag()) AddNoisyPixels(); + if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels(); +// cout << "After Remove in DigitiseModule in module " << mod->GetIndex() << endl; // dom + FrompListToDigits(); + ClearMap(); +} +//______________________________________________________________________ +void AliITSsimulationSPD::HitToSDigit(AliITSmodule *mod){ + // Does the charge distributions using Gaussian diffusion charge charing. + // Inputs: + // AliITSmodule *mod Pointer to this module + // Output: + // none. + // Return: + // none. + const Double_t kmictocm = 1.0e-4; // convert microns to cm. + const Double_t kBunchLenght = 25e-9; // LHC clock + TObjArray *hits = mod->GetHits(); + Int_t nhits = hits->GetEntriesFast(); + Int_t h,ix,iz,i; + Int_t idtrack; + Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0; + Double_t x,y,z,t,tp,st,dt=0.2,el,sig,sigx,sigz,fda; + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + AliITSSimuParam *simpar = fDetType->GetSimuParam(); + Double_t thick = 0.5*kmictocm*seg->Dy(); // Half Thickness + simpar->GetSPDSigmaDiffusionAsymmetry(fda); // + + AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling)); + if(nhits<=0) return; + for(h=0;h0) { + AliDebug(1,Form("Hits, %d", h)); + cout << *(mod->GetHit(h)) << endl; + } // end if GetDebug + // Check if the hit is inside readout window + if (fStrobe) + if ((mod->GetHit(h)->GetTOF() < fStrobePhase) || + (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue; + if(!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; + st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); + if(st>0.0){ + st = (Double_t)((Int_t)(st/kmictocm)); // number of microns + if(st<=1.0) st = 1.0; + dt = 1.0/st; + for(t=0.0;t<1.0;t+=dt){ // Integrate over t + tp = t+0.5*dt; + x = x0+x1*tp; + y = y0+y1*tp; + z = z0+z1*tp; + if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside + //el = res->GeVToCharge((Double_t)(dt*de)); + el = dt * de / simpar->GetGeVToCharge(); + if(GetDebug(1)){ + if(el<=0.0) cout<<"el="<GetHits(); + Int_t nhits = hits->GetEntriesFast(); + Int_t h,ix,iz,i; + Int_t idtrack; + Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0; + Double_t x,y,z,t,st,el,sig,sigx,sigz,fda; + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + AliITSSimuParam* simpar = fDetType->GetSimuParam(); + Double_t thick = 0.5*kmictocm*seg->Dy(); // Half thickness + simpar->GetSPDSigmaDiffusionAsymmetry(fda); +// cout << "Half Thickness " << thick << endl; // dom +// cout << "Diffusion asymm " << fda << endl; // dom + + AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling)); + if(nhits<=0) return; + for(h=0;h0) { + AliDebug(1,Form("Hits, %d", h)); + cout << *(mod->GetHit(h)) << endl; + } // end if GetDebug + // Check if the hit is inside readout window + if (fStrobe) + if ((mod->GetHit(h)->GetTOF() < fStrobePhase) || + (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue; + if(!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; + st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); + if(st>0.0) for(i=0;iLocalToDet(x,z,ix,iz))) continue; // outside + el = kwi[i]*de/simpar->GetGeVToCharge(); + if(GetDebug(1)){ + if(el<=0.0) cout<<"el="<LocalToDet(x,z,ix,iz))) continue; // outside + el = de / simpar->GetGeVToCharge(); + sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y)); + sigx=sig; + sigz=sig*fda; + if (fLorentz) ld=(y+thick)*fTanLorAng; + SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); + } // end if st>0.0 + + } // Loop over all hits h + + // Coupling + switch (fCoupling) { + default: + break; + case 1: // case 3: + for(i=0;iGetEntries();i++) + if(GetMap()->GetpListItem(i)==0) continue; + else{ + GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix); + SetCoupling(iz,ix); + } // end for i + break; + case 2: // case 4: + for(i=0;iGetEntries();i++) + if(GetMap()->GetpListItem(i)==0) continue; + else{ + GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix); + SetCouplingOld(iz,ix); + } // end for i + break; + } // end switch + if(GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling); +} +//______________________________________________________________________ +void AliITSsimulationSPD::SpreadCharge(Double_t x0,Double_t z0, + Int_t ix0,Int_t iz0, + Double_t el,Double_t sig,Double_t ld, + Int_t t,Int_t hi){ + // Spreads the charge over neighboring cells. Assume charge is distributed + // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg) + // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig) + // if fLorentz=kTRUE, then x0=x0+ld (Lorentz drift taken into account) + // Defined this way, the integral over all x and z is el. + // Inputs: + // Double_t x0 x position of point where charge is liberated + // Double_t z0 z position of point where charge is liberated + // Int_t ix0 row of cell corresponding to point x0 + // Int_t iz0 columb of cell corresponding to point z0 + // Double_t el number of electrons liberated in this step + // Double_t sig Sigma difusion for this step (y0 dependent) + // Double_t ld lorentz drift in x for this step (y0 dependent) + // Int_t t track number + // Int_t ti hit track index number + // Int_t hi hit "hit" index number + // Outputs: + // none. + // Return: + // none. + const Int_t knx = 3,knz = 2; + const Double_t kRoot2 = 1.414213562; // Sqrt(2). + const Double_t kmictocm = 1.0e-4; // convert microns to cm. + Int_t ix,iz,ixs,ixe,izs,ize; + Float_t x,z; + Double_t x1,x2,z1,z2,s,sp; + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + + + if(GetDebug(4)) Info("SpreadCharge","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e," + "sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi); + if(sig<=0.0) { // if sig<=0 No diffusion to simulate. + GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el); + if(GetDebug(2)){ + cout << "sig<=0.0=" << sig << endl; + } // end if GetDebug + return; + } // end if + sp = 1.0/(sig*kRoot2); + if(GetDebug(2)){ + cout << "sig=" << sig << " sp=" << sp << endl; + } // end if GetDebug + ixs = TMath::Max(-knx+ix0,0); + ixe = TMath::Min(knx+ix0,seg->Npx()-1); + izs = TMath::Max(-knz+iz0,0); + ize = TMath::Min(knz+iz0,seg->Npz()-1); + for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){ + seg->DetToLocal(ix,iz,x,z); // pixel center + x1 = x; + z1 = z; + x2 = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper + x1 -= 0.5*kmictocm*seg->Dpx(ix); // Lower + z2 = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper + z1 -= 0.5*kmictocm*seg->Dpz(iz); // Lower + x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) + x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) + z1 -= z0; // Distance from where track traveled + z2 -= z0; // Distance from where track traveled + s = 0.25; // Correction based on definision of Erfc + s *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2); + if(GetDebug(3)){ + cout <<"el="<AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el); + if(GetDebug(2)){ + cout << "sigx<=0.0=" << sigx << endl; + cout << "sigz<=0.0=" << sigz << endl; + } // end if GetDebug + return; + } // end if + spx = 1.0/(sigx*kRoot2); spz = 1.0/(sigz*kRoot2); + if(GetDebug(2)){ + cout << "sigx=" << sigx << " spx=" << spx << endl; + cout << "sigz=" << sigz << " spz=" << spz << endl; + } // end if GetDebug + ixs = TMath::Max(-knx+ix0,0); + ixe = TMath::Min(knx+ix0,seg->Npx()-1); + izs = TMath::Max(-knz+iz0,0); + ize = TMath::Min(knz+iz0,seg->Npz()-1); + for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){ + seg->DetToLocal(ix,iz,x,z); // pixel center + x1 = x; + z1 = z; + x2 = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper + x1 -= 0.5*kmictocm*seg->Dpx(ix); // Lower + z2 = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper + z1 -= 0.5*kmictocm*seg->Dpz(iz); // Lower + x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) + x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) + z1 -= z0; // Distance from where track traveled + z2 -= z0; // Distance from where track traveled + s = 0.25; // Correction based on definision of Erfc + s *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2); + if(GetDebug(3)){ + cout <<"el="<MinVal(); - - Int_t digits[3], tracks[3], hits[3],gi,j1; - Float_t charges[3]; - Float_t electronics; - Float_t signal,phys; - for(Int_t iz=0;izGaus(); - signal = (float)fMapA2->GetSignal(iz,ix); - signal += electronics; - gi =iz*fNPixelsX+ix; // global index - if (signal > threshold) { - digits[0]=iz; - digits[1]=ix; - digits[2]=1; - for(j1=0;j1<3;j1++){ - if (pList[gi]) { - //b.b. tracks[j1]=-3; - tracks[j1] = (Int_t)(*(pList[gi]+j1)); - hits[j1] = (Int_t)(*(pList[gi]+j1+6)); - }else { - tracks[j1]=-2; //noise - hits[j1] = -1; - } - charges[j1] = 0; - } - - if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) { - tracks[1] = -3; - hits[1] = -1; - tracks[2] = -3; - hits[2] = -1; - } - if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) { - tracks[1] = -3; - hits[1] = -1; - } - if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) { - tracks[2] = -3; - hits[2] = -1; - } - if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) { - tracks[2] = -3; - hits[2] = -1; - } - - phys=0; - aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges); - } - if(pList[gi]) delete [] pList[gi]; - } +//______________________________________________________________________ +void AliITSsimulationSPD::AddNoisyPixels() { + // Adds noisy pixels on each module (ladder) + // This should be called before going from sdigits to digits (FrompListToDigits) + Int_t mod = GetModuleNumber(); + AliITSCalibrationSPD* calObj = (AliITSCalibrationSPD*) fDetType->GetSPDNoisyModel(mod); + + Int_t nrNoisy = calObj->GetNrBad(); + for (Int_t i=0; iAddNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), mod, 10*GetThreshold()); } - delete [] pList; - } - - -//____________________________________________ - -void AliITSsimulationSPD::CreateHistograms() -{ - // create 1D histograms for tests - - printf("SPD - create histograms\n"); - - fHis=new TObjArray(fNPixelsZ); - TString spdName("spd_"); - for (Int_t i=0;iGetModule("ITS"); + Int_t j,ix,iz; + Double_t electronics; + Double_t sig; + const Int_t knmaxtrk=AliITSdigit::GetNTracks(); + static AliITSdigitSPD dig; + AliITSSimuParam *simpar = fDetType->GetSimuParam(); + if(GetDebug(1)) Info("FrompListToDigits","()"); + for(iz=0; izRndm() >= eff) continue; +// } // end if +// END parametrize the efficiency +// + electronics = simpar->ApplySPDBaselineAndNoise(); + UpdateMapNoise(ix,iz,electronics); + // + // Apply Threshold and write Digits. + sig = GetMap()->GetSignalOnly(iz,ix); + FillHistograms(ix,iz,sig+electronics); + if(GetDebug(3)){ + cout<threshold("<GetNEntries()) { + dig.SetTrack(j,GetMap()->GetTrack(iz,ix,j)); + dig.SetHit(j,GetMap()->GetHit(iz,ix,j)); + }else { // Default values + dig.SetTrack(j,-3); + dig.SetHit(j,-1); + } // end if GetMap() + } // end for j + if(GetDebug(3)){ + cout<GetpListItem(iz,ix))<AddSimDigit(0,&dig); + // simulate fo signal response for this pixel hit: + fDetType->ProcessSPDDigitForFastOr(fModule, dig.GetCoord1(), dig.GetCoord2()); + } // for ix/iz +} +//______________________________________________________________________ +void AliITSsimulationSPD::CreateHistograms(){ + // create 1D histograms for tests + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + + if(GetDebug(1)) Info("CreateHistograms","create histograms"); + + fHis = new TObjArray(GetNPixelsZ()); + fSPDname="spd_"; + for(Int_t i=0;iAddAt(new TH1F(fSPDname.Data(),"SPD maps", + GetNPixelsX(),0.,(Double_t)GetNPixelsX()),i); + } // end for i +} +//______________________________________________________________________ +void AliITSsimulationSPD::FillHistograms(Int_t ix,Int_t iz,Double_t v){ + // Fill the histogram + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + + if(!GetHistArray()) return; // Only fill if setup. + if(GetDebug(2)) Info("FillHistograms","fill histograms"); + GetHistogram(iz)->Fill(ix,v); +} +//______________________________________________________________________ +void AliITSsimulationSPD::ResetHistograms(){ + // Reset histograms for this detector + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + + if(!GetHistArray()) return; // Only fill if setup. + if(GetDebug(2)) Info("FillHistograms","fill histograms"); + for ( int i=0;iAt(i)) ((TH1F*)fHis->At(i))->Reset(); + } // end for i } -//____________________________________________ - -void AliITSsimulationSPD::ResetHistograms() +//______________________________________________________________________ +void AliITSsimulationSPD::SetCoupling(Int_t col, Int_t row) { + // Take into account the coupling between adiacent pixels. + // The parameters probcol and probrow are the probability of the + // signal in one pixel shared in the two adjacent pixels along + // the column and row direction, respectively. + // Note pList is goten via GetMap() and module is not need any more. + // Otherwise it is identical to that coded by Tiziano Virgili (BSN). + //Begin_Html + /* + + +
+ + . + + 
+   */
+   //End_Html
+   // Inputs:
+   //    Int_t col            z cell index
+   //    Int_t row            x cell index
+   // Outputs:
+   //    none.
+   // Return:
+   //     none.
+   Int_t j1,j2,flag=0;
+   Double_t pulse1,pulse2;
+   Double_t couplR=0.0,couplC=0.0;
+   Double_t xr=0.;
+
+   GetCouplings(couplC,couplR);
+   if(GetDebug(3)) Info("SetCoupling","(col=%d,row=%d) "
+                        "Calling SetCoupling couplC=%e couplR=%e",
+                        col,row,couplC,couplR);
+   j1 = col;
+   j2 = row;
+   pulse1 = GetMap()->GetSignalOnly(col,row);
+   pulse2 = pulse1;
+   for (Int_t isign=-1;isign<=1;isign+=2){// loop in col direction
+       do{
+           j1 += isign;
+           xr = gRandom->Rndm();
+           if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplC)){
+               j1 = col;
+               flag = 1;
+           }else{
+               UpdateMapNoise(row,j1,pulse1);
+               //  flag = 0;
+               flag = 1; // only first next!!
+           } // end if
+       } while(flag == 0);
+       // loop in row direction
+       do{
+           j2 += isign;
+           xr = gRandom->Rndm();
+           if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplR)){
+               j2 = row;
+               flag = 1;
+           }else{
+               UpdateMapNoise(j2,col,pulse2);
+               //  flag = 0;
+               flag = 1; // only first next!!
+           } // end if
+       } while(flag == 0);
+   } // for isign
+}
+//______________________________________________________________________
+void AliITSsimulationSPD::SetCouplingOld(Int_t col, Int_t row) {
+   //  Take into account the coupling between adiacent pixels.
+   //  The parameters probcol and probrow are the fractions of the
+   //  signal in one pixel shared in the two adjacent pixels along
+   //  the column and row direction, respectively.
+   //Begin_Html
+   /*
+     
+
+
+ + . + + 
+   */
+   //End_Html
+   // Inputs:
+   //    Int_t col            z cell index
+   //    Int_t row            x cell index
+   //    Int_t module         module number
+   // Outputs:
+   //    none.
+   // Return:
+   //     none.
+   Int_t j1,j2,flag=0;
+   Double_t pulse1,pulse2;
+   Double_t couplR=0.0,couplC=0.0;
+
+   GetCouplings(couplC,couplR);
+
+   //  Debugging ...
+//    cout << "Threshold --> " << GetThreshold() << endl;  // dom
+//    cout << "Couplings --> " << couplC << " " << couplR << endl;  // dom
+
+   if(GetDebug(3)) Info("SetCouplingOld","(col=%d,row=%d) "
+                        "Calling SetCoupling couplC=%e couplR=%e",
+                        col,row,couplC,couplR);
+   for (Int_t isign=-1;isign<=1;isign+=2){// loop in col direction
+   pulse1 = GetMap()->GetSignalOnly(col,row);
+   pulse2 = pulse1;
+   j1 = col;
+   j2 = row;
+       do{
+           j1 += isign;
+           pulse1 *= couplC;
+           if ((j1<0)||(j1>GetNPixelsZ()-1)||(pulse1GetSignalOnly(col,row);
+               j1 = col;
+               flag = 1;
+           }else{
+               UpdateMapNoise(row,j1,pulse1);
+               // flag = 0;
+               flag = 1;  // only first next !!
+           } // end if
+       } while(flag == 0);
+       // loop in row direction
+       do{
+           j2 += isign;
+           pulse2 *= couplR;
+           if((j2<0)||(j2>(GetNPixelsX()-1))||(pulse2GetSignalOnly(col,row);
+               j2 = row;
+               flag = 1;
+           }else{
+               UpdateMapNoise(j2,col,pulse2);
+               // flag = 0;
+               flag = 1; // only first next!!
+           } // end if
+       } while(flag == 0);
+   } // for isign
+}
+//______________________________________________________________________
+void AliITSsimulationSPD::GenerateStrobePhase()
 {
-    //
-    // Reset histograms for this detector
-    //
-
-    for ( int i=0;iReset();
-    }
-
+ // Generate randomly the strobe
+ // phase w.r.t to the LHC clock
+ // Done once per event
+   const Double_t kBunchLenght = 25e-9; // LHC clock
+   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
+                  (Double_t)fStrobeLenght*kBunchLenght+
+                  kBunchLenght/2;
 }
 
-
-
-
-
-
-
-
-