X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSsimulationSPD.cxx;h=d2308f7bc1edc4e152cb0f85f99d360f3321d74f;hb=790a8a943dfd39ce2939493300fe69954712c249;hp=0befd6bfb1bebd4aaed814e0cf07ab4958f911f5;hpb=9f0330014970465b2e5e924cac1813479fb2c489;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSPD.cxx b/ITS/AliITSsimulationSPD.cxx index 0befd6bfb1b..d2308f7bc1e 100644 --- a/ITS/AliITSsimulationSPD.cxx +++ b/ITS/AliITSsimulationSPD.cxx @@ -1,634 +1,969 @@ -#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 "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 "AliCDBEntry.h" +#include "AliCDBLocal.h" +//#define DEBUG 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 // -// AliITSsimulationSPD is 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(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(); +// Version: 0 +// Written by Boris Batyunya +// December 20 1999 +// +// +// AliITSsimulationSPD is to do the simulation of SPDs. +// +//////////////////////////////////////////////////////////////////////// +//______________________________________________________________________ +AliITSsimulationSPD::AliITSsimulationSPD(): +AliITSsimulation(), +fHis(0), +fSPDname(), +fCoupling(){ + // Default constructor. + // Inputs: + // none. + // Outputs: + // none. + // Return: + // A default constructed AliITSsimulationSPD class. + + AliDebug(1,Form("Calling default constructor")); +// Init(); } - -//_____________________________________________________________________________ - -AliITSsimulationSPD::~AliITSsimulationSPD() { - // destructor - - delete fMapA2; - - if (fHis) { - fHis->Delete(); - delete fHis; - } +//______________________________________________________________________ +AliITSsimulationSPD::AliITSsimulationSPD(AliITSDetTypeSim *dettyp): +AliITSsimulation(dettyp), +fHis(0), +fSPDname(), +fCoupling(){ + // 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 ")); +// AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); +// res->SetTemperature(0.0); +// res->SetDistanceOverVoltage(0.0); + Init(); } +//______________________________________________________________________ +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())); + AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + Double_t bias = res->GetBiasVoltage(); +// cout << "Bias Voltage --> " << bias << endl; // dom + res->SetDistanceOverVoltage(kmictocm*seg->Dy(),bias); +// set kind of coupling ("old" or "new") + char opt[20]; + res->GetCouplingOption(opt); + char *old = strstr(opt,"old"); + if (old) { + fCoupling=2; + } else { + fCoupling=1; + } // end if + + // Get the calibration objects for each module(ladder) + GetCalibrationObjects(0); //RunNr 0 hard coded for now - -//__________________________________________________________________________ -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; } +//______________________________________________________________________ +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){ + // 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 +} +//______________________________________________________________________ +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; + 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; +} - Float_t spdLength = fSegmentation->Dz(); - Float_t spdWidth = fSegmentation->Dx(); +//______________________________________________________________________ +void AliITSsimulationSPD::GetCalibrationObjects(Int_t RunNr) { + // Gets the calibration objects for each module (ladder) + // Inputs: + // RunNr: hard coded to RunNr=0 for now + // Outputs: + // none. + // Return: + // none. - Float_t difCoef, dum; - fResponse->DiffCoeff(difCoef,dum); + AliCDBManager* man = AliCDBManager::Instance(); - Float_t zPix0 = 1e+6; - Float_t xPix0 = 1e+6; - Float_t yPix0 = 1e+6; - Float_t yPrev = 1e+6; - Float_t zP0 = 100.; - Float_t xP0 = 100.; + AliCDBEntry *entrySPD=0; + entrySPD = man->Get("ITS/Calib/CalibSPD", RunNr); - Float_t zPitch = fSegmentation->Dpz(0); - Float_t xPitch = fSegmentation->Dpx(0); + if(!entrySPD){ + AliFatal("Cannot find SPD calibration entry in default storage!"); + return; + } - //cout << "pitch per z: " << zPitch << endl; - //cout << "pitch per r*phi: " << xPitch << endl; - - TObjArray *fHits = mod->GetHits(); - Int_t nhits = fHits->GetEntriesFast(); - if (!nhits) return; - // cout << "module, nhits ="<At(hit); - Int_t layer = iHit->GetLayer(); - - // 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 : - - 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(TMath::Abs(zPix) > spdLength/2.) { - printf("!! Zpix outside = %f\n",zPix); - if(status == 66) zP0=100; - continue; - } - - - if (TMath::Abs(xPix) > spdWidth/2.) { - printf("!! Xpix outside = %f\n",xPix); - if (status == 66) xP0=100; - continue; - } - - Float_t zP = (zPix + spdLength/2.)/1000.; - Float_t xP = (xPix + spdWidth/2.)/1000.; - - Int_t trdown = 0; - - // enter Si or after event in Si - if (status == 66 ) { - zPix0 = zPix; - xPix0 = xPix; - yPrev = yPix; - } - // enter Si only - if (layer == 1 && status == 66 && yPix > 71.) { - yPix0 = yPix; - zP0 = zP; - xP0 = xP; - } - // enter Si only - if (layer == 2 && status == 66 && yPix < -71.) { - yPix0 = yPix; - zP0 = zP; - xP0 = xP; - } - Float_t depEnergy = iHit->GetIonization(); - // skip if the input point to Si - if(depEnergy <= 0.) continue; - // skip if the input point is outside of Si, but the next - // point is inside of Si - if(zP0 > 90 || xP0 > 90) continue; - // if track returns to the opposite direction: - if (layer == 1 && yPix > yPrev) { - yPix0 = yPrev; - trdown = 1; - } - if (layer == 2 && yPix < yPrev) { - yPix0 = 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 = yPix - yPix0; - - if(TMath::Abs(ydif) < 0.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(TMath::Abs(projDif) > 5.0) tang = ydif/projDif; - Float_t dCharge = charge/nsteps; // charge in e- for one step - Float_t dZ = zdif/nsteps; - Float_t dX = xdif/nsteps; - - if (TMath::Abs(projDif) < 5.0 ) { - drPath = ydif*1.e-4; - drPath = TMath::Abs(drPath); // drift path in cm - sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm - } - - 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(TMath::Abs(projDif) >= 5.) { - Float_t dProjn = sqrt(dZn*dZn+dXn*dXn); - if(trdown == 0) { - drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm - drPath = TMath::Abs(drPath); - } - if(trdown == 1) { - Float_t dProjn = projDif/nsteps; - drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4; - 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 - - - //set protection - for(int k=0;k<4;k++) { - if (indexRange[k] < 0) indexRange[k]=0; + TObjArray *respSPD = (TObjArray *)entrySPD->GetObject(); + if ((! respSPD)) { + AliFatal("Cannot get data from SPD database entry!"); + return; + } + for (Int_t mod=0; mod<240; mod++) { + fCalObj[mod] = (AliITSCalibrationSPD*) respSPD->At(mod); } - - for(Int_t iz=indexRange[0];izGetSignal(iz,ix); - - if (!signal) continue; - - Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0! - if(!pList[globalIndex]){ - - // - // Create new list (9 elements - 3 signals and 3 tracks + 3 hits) - // - - pList[globalIndex] = new Float_t [9]; - - // set list to -3 - - *pList[globalIndex] = -3.; - *(pList[globalIndex]+1) = -3.; - *(pList[globalIndex]+2) = -3.; - *(pList[globalIndex]+3) = 0.; - *(pList[globalIndex]+4) = 0.; - *(pList[globalIndex]+5) = 0.; - *(pList[globalIndex]+6) = -1.; - *(pList[globalIndex]+7) = -1.; - *(pList[globalIndex]+8) = -1.; - - - *pList[globalIndex] = (float)label; - *(pList[globalIndex]+3) = signal; - *(pList[globalIndex]+6) = (float)idhit; - } - else{ - - // check the signal magnitude - - Float_t highest = *(pList[globalIndex]+3); - Float_t middle = *(pList[globalIndex]+4); - Float_t lowest = *(pList[globalIndex]+5); - - signal -= (highest+middle+lowest); - - // - // compare the new signal with already existing list - // - - if(signalhighest){ - *(pList[globalIndex]+5) = middle; - *(pList[globalIndex]+4) = highest; - *(pList[globalIndex]+3) = signal; - - *(pList[globalIndex]+2) = *(pList[globalIndex]+1); - *(pList[globalIndex]+1) = *pList[globalIndex]; - *pList[globalIndex] = label; - - *(pList[globalIndex]+8) = *(pList[globalIndex]+7); - *(pList[globalIndex]+7) = *(pList[globalIndex]+6); - *(pList[globalIndex]+6) = idhit; - } - else if (signal>middle){ - *(pList[globalIndex]+5) = middle; - *(pList[globalIndex]+4) = signal; - - *(pList[globalIndex]+2) = *(pList[globalIndex]+1); - *(pList[globalIndex]+1) = label; - - *(pList[globalIndex]+8) = *(pList[globalIndex]+7); - *(pList[globalIndex]+7) = idhit; - } - else{ - *(pList[globalIndex]+5) = signal; - *(pList[globalIndex]+2) = label; - *(pList[globalIndex]+8) = idhit; - } - } - } // end of loop pixels in x - } // end of loop over pixels in z - - } - -//--------------------------------------------- -void AliITSsimulationSPD::ChargeToSignal(Float_t **pList) -{ - // add noise and electronics, perform the zero suppression and add the - // digit to the list - - AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS"); - - - TRandom *random = new TRandom(); - Float_t threshold = (float)fResponse->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]) { - 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::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)); + SetModuleNumber(module); + SetEventNumber(event); + 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()); + SetEventNumber(event); + // HitToSDigit(mod); + HitToSDigitFast(mod); + RemoveDeadPixels(mod); +// 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){ + // 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. + + AliDebug(1,Form("(mod=%p,,)",mod)); + // HitToSDigit(mod); + HitToSDigitFast(mod); + RemoveDeadPixels(mod); +// 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. + 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; + Double_t x,y,z,t,tp,st,dt=0.2,el,sig,sigx,sigz,fda; + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); + Double_t thick = 0.5*kmictocm*seg->Dy(); // Half Thickness + res->GetSigmaDiffusionAsymmetry(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 + 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)); + if(GetDebug(1)){ + if(el<=0.0) cout<<"el="<LocalToDet(x,z,ix,iz))) continue; // outside + el = res->GeVToCharge((Double_t)de); + sig = res->SigmaDiffusion1D(TMath::Abs(thick + y)); + // SpreadCharge(x,z,ix,iz,el,sig,idtrack,h); + sigx=sig; + sigz=sig*fda; + SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,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,idtrack,h); + } // 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,idtrack,h); + } // end for i + break; + } // end switch + if(GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling); +} +//______________________________________________________________________ +void AliITSsimulationSPD::HitToSDigitFast(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 Int_t kn10=10; + const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1, + 4.325316833e-1,4.869532643e-1,5.130467358e-1, + 5.674683167e-1,6.602952159e-1,7.833023029e-1, + 9.255628306e-1}; + const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1, + 7.472567455e-2,3.333567215e-2,3.333567215e-2, + 7.472567455e-2,1.095431813e-1,1.346333597e-1, + 1.477621124e-1}; + 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; + Double_t x,y,z,t,st,el,sig,sigx,sigz,fda; + AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0); + AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); + Double_t thick = 0.5*kmictocm*seg->Dy(); // Half thickness + res->GetSigmaDiffusionAsymmetry(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 + 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 = res->GeVToCharge((Double_t)(dt*de)); + // el = 1./kn10*res->GeVToCharge((Double_t)de); + el = kwi[i]*res->GeVToCharge((Double_t)de); + if(GetDebug(1)){ + if(el<=0.0) cout<<"el="<LocalToDet(x,z,ix,iz))) continue; // outside + el = res->GeVToCharge((Double_t)de); + sig = res->SigmaDiffusion1D(TMath::Abs(thick + y)); + //SpreadCharge(x,z,ix,iz,el,sig,idtrack,h); + sigx=sig; + sigz=sig*fda; + SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,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,idtrack,h); + } // 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,idtrack,h); + } // 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,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) + // 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 y0 y 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) + // 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; // Distance from where track traveled + x2 -= x0; // Distance from where track traveled + z1 -= z0; // Distance from where track traveled + z2 -= z0; // Distance from where track traveled + s = 0.25; // Correction based on definision of Erfc + s *= TMath::Erfc(sp*x1) - TMath::Erfc(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; // Distance from where track traveled + x2 -= x0; // Distance from where track traveled + z1 -= z0; // Distance from where track traveled + z2 -= z0; // Distance from where track traveled + s = 0.25; // Correction based on definision of Erfc + s *= TMath::Erfc(spx*x1) - TMath::Erfc(spx*x2); + if(GetDebug(3)){ + cout <<"el="<GetModule("ITS"); + Int_t j,ix,iz; + Double_t electronics; + Double_t sig; + const Int_t knmaxtrk=AliITSdigit::GetNTracks(); + static AliITSdigitSPD dig; + AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); + if(GetDebug(1)) Info("FrompListToDigits","()"); + for(iz=0; izRndm() >= eff) continue; +// } // end if + // END parametrize the efficiency + // + electronics = res->ApplyBaselineAndNoise(); + 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); + } // for ix/iz } - -//____________________________________________ - -void AliITSsimulationSPD::ResetHistograms() -{ - // +//______________________________________________________________________ +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()); + TString 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 - // - for ( int i=0;iReset(); - } - + // 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::SetCoupling(Int_t col, Int_t row, Int_t ntrack, + Int_t idhit) { + // 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
+    //    Int_t ntrack         track incex number
+    //    Int_t idhit          hit index number
+    // 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,ntrack=%d,idhit=%d) "
+                         "Calling SetCoupling couplC=%e couplR=%e",
+                         col,row,ntrack,idhit,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{
+                UpdateMapSignal(row,j1,ntrack,idhit,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{
+                UpdateMapSignal(j2,col,ntrack,idhit,pulse2);
+                //  flag = 0;
+                flag = 1; // only first next!!
+            } // end if
+        } while(flag == 0);
+    } // for isign
+}
+//______________________________________________________________________
+void AliITSsimulationSPD::SetCouplingOld(Int_t col, Int_t row,
+                Int_t ntrack,Int_t idhit) {
+    //  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 ntrack         track incex number
+    //    Int_t idhit          hit index number
+    //    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,ntrack=%d,idhit=%d) "
+                         "Calling SetCoupling couplC=%e couplR=%e",
+                         col,row,ntrack,idhit,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{
+                UpdateMapSignal(row,j1,ntrack,idhit,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{
+                UpdateMapSignal(j2,col,ntrack,idhit,pulse2);
+                // flag = 0;
+                flag = 1; // only first next!!
+            } // end if
+        } while(flag == 0);
+    } // for isign
+}