X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSsimulationSPD.cxx;h=a7e3ee97ed09357512e5360639554c417d6576b8;hb=8809be5b5e578480560070566463e843d4d9f9b3;hp=bec9375d1f24311f2d39b0d0cd09ed259884f857;hpb=03898a572903e84d489de082f1784fc61512613d;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSPD.cxx b/ITS/AliITSsimulationSPD.cxx index bec9375d1f2..a7e3ee97ed0 100644 --- a/ITS/AliITSsimulationSPD.cxx +++ b/ITS/AliITSsimulationSPD.cxx @@ -1,625 +1,987 @@ -#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 "AliMagF.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 +// October 2007: GetCalibrationObjects() removed // -// 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(), +fLorentz(kFALSE), +fTanLorAng(0){ + // 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(), +fLorentz(kFALSE), +fTanLorAng(0){ + // 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(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; +//______________________________________________________________________ +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(kTRUE); + if (fLorentz) SetTanLorAngle(); } +//______________________________________________________________________ +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 *mf = gAlice->Field(); + Float_t pos[3]={0.,0.,0.}; + Float_t B[3]={0.,0.,0.}; + mf->Field(pos,B); + fTanLorAng = TMath::Tan(WeightHole*simpar->LorentzAngleHole(B[2]) + + WeightEle*simpar->LorentzAngleElectron(B[2])); + fTanLorAng*=-1.; // this only for the old geometry + // comment the upper line for the new geometry + 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){ + // 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; - - //cout<<"len,wid,dy,nx,nz,pitchx,pitchz ="<Dy()<<","<At(hit); - Int_t layer = iHit->GetLayer(); - Float_t yPix0 = -73; - if(layer == 1) yPix0 = -77; - - 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(); // the momenta at the - Float_t py = iHit->GetPYL(); // each GEANT step - Float_t pz = iHit->GetPZL(); - Float_t ptot = 1000*sqrt(px*px+py*py+pz*pz); - */ - - Float_t pmod = iHit->GetParticle()->P(); // total momentum at the - // vertex - pmod *= 1000; - - - 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(); - //cout<<"hit,status,y ="< spdLength/2) { - //cout<<"!!!1 z outside ="<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()); + SetEventNumber(event); + InitSimulationModule( GetModuleNumber() , 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 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. + + AliDebug(1,Form("(mod=%p,,)",mod)); + // HitToSDigit(mod); + InitSimulationModule( mod->GetIndex(), event ); + 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,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 + 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 + 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,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,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 *= 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+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 *= TMath::Erfc(spx*x1) - TMath::Erfc(spx*x2); + if(GetDebug(3)){ + cout <<"el="< " << moduleNr << endl; // dom +// cout << "nr of dead " << nrDead << endl; // dom + for (Int_t i=0; iDeleteHit(calObj->GetBadColAt(i),calObj->GetBadRowAt(i)); +// cout << "dead index " << i << endl; // dom +// cout << "col row --> " << calObj->GetDeadColAt(i) << " " << calObj->GetDeadRowAt(i) << endl; // dom } - - 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; - } +} +//______________________________________________________________________ +void AliITSsimulationSPD::FrompListToDigits(){ + // add noise and electronics, perform the zero suppression and add the + // digit to the list + // Inputs: + // none. + // Outputs: + // none. + // Return: + // none. + static AliITS *aliITS = (AliITS*)gAlice->GetModule("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); + } // for ix/iz } - - -//--------------------------------------------- -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"); - - - 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]) { - //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]; - } - } - delete [] pList; - +//______________________________________________________________________ +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::CreateHistograms() -{ - // create 1D histograms for tests - - printf("SPD - create histograms\n"); - - fHis=new TObjArray(fNPixelsZ); - for (Int_t i=0;iFill(ix,v); } - -//____________________________________________ - -void AliITSsimulationSPD::ResetHistograms() -{ - // +//______________________________________________________________________ +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
+}