X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSsimulationSPD.cxx;h=7b32e62287871d7eb883754b396859f6318f2430;hb=c61a7285011faf50a1f39c6a49e0c3fdeeff6260;hp=af2d4af650d7c039621b9236fc38446828f4e9eb;hpb=b0f5e3fcf84ecedb50df35ebe629b44bafc15bc0;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSsimulationSPD.cxx b/ITS/AliITSsimulationSPD.cxx index af2d4af650d..7b32e622878 100644 --- a/ITS/AliITSsimulationSPD.cxx +++ b/ITS/AliITSsimulationSPD.cxx @@ -1,541 +1,984 @@ -#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 "AliRun.h" -#include "AliITSMap.h" // "AliITSMapA2.h" +#include +#include "AliITS.h" +#include "AliITSdigitSPD.h" +#include "AliITShit.h" +#include "AliITSmodule.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 +// Version: 0 +// Written by Boris Batyunya +// December 20 1999 // -//________________________________________________________________________ - -AliITSsimulationSPD::AliITSsimulationSPD(){ - // constructor - fResponse = 0; - fSegmentation = 0; - fHis = 0; - fNoise = 0.; - fBaseline = 0.; -} -//_____________________________________________________________________________ - -AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) { - // constructor - fResponse = resp; - fSegmentation = seg; - - fResponse->GetNoiseParam(fNoise,fBaseline); - - - fMapA2 = new AliITSMapA2(fSegmentation); - - // - fNPixelsZ=fSegmentation->Npz(); - fNPixelsX=fSegmentation->Npx(); +// +// 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 ")); +// AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); +// res->SetTemperature(0.0); +// res->SetDistanceOverVoltage(0.0); + 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())); + 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 + //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; + AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD()); + 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*res->LorentzAngleHole(B[2]) + + WeightEle*res->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; +} +//______________________________________________________________________ +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; + 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)); + 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,ld=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="<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); + 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)); + 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 } -//_____________________________________________________________________________ - -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 spdLenght = fSegmentation->Dz(); - Float_t spdWidth = fSegmentation->Dx(); - - Float_t difCoef = fResponse->DiffCoeff(); - - Float_t zPix0 = 1e+6; - Float_t xPix0 = 1e+6; - Float_t yPix0 = 1e+6; - Float_t yPrev = 1e+6; - - Float_t zPitch = fSegmentation->Dpz(0); - Float_t xPitch = fSegmentation->Dpx(0); - - //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; - - - - // Array of pointers to the label-signal list - - Int_t maxNdigits = fNPixelsX*fNPixelsZ; - 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 - - Int_t layer,idtrack,status,trdown,ndZ,ndX,nsteps,iZi,nZpix,nXpix; - Int_t jzmin,jzmax,jxmin,jxmax,jx,jz,lz,lx,index; - Float_t zArg1,zArg2,xArg1,xArg2; - Float_t zProb1,zProb2,xProb1,xProb2; - Float_t dZCharge,dXCharge; - Double_t signal; - Float_t projDif; // RMS of xDif and zDif - Float_t dProjn; // RMS of dXn and dZn - Float_t depEnergy; // The deposited energy from this hit - Float_t zPix; // hit position in microns - Float_t xPix; // hit position in microns - Float_t yPix; // hit position in microns - Float_t zPixn; // hit position in microns - Float_t xPixn; // hit position in microns - Float_t charge,dCharge;// charge in e- - Float_t drPath; - Float_t tAng,dZ,dX,dXn,dZn; - Float_t sigmaDif; - Float_t dZright,dZleft; - Float_t dXright,dXleft; - Float_t dZprev,dZnext,dXprev,dXnext; - - static Bool_t first=kTRUE; - Int_t lasttrack = -2,hit; - for(hit=0;hitAt(hit); - layer = iHit->GetLayer(); - - // work with the idtrack=entry number in the TreeH - // Int_t idtrack=mod->GetHitTrackIndex(ii); - // or store straight away the particle position in the array - // of particles : - idtrack = iHit->GetTrack(); - - //if(module==11 || module==157) { - // Int_t track = iHit->fTrack; - // Int_t primary = gAlice->GetPrimary(track); - // Int_t parent = iHit->GetParticle()->GetFirstMother(); - // printf("module,hit,track,primary,parent %d %d %d %d %d \n", - // md,hit,track,primary,parent); - //} - - // Get hit z and x(r*phi) cordinates for each module (detector) - // in local system. - - zPix = kconv*iHit->GetZL(); // Geant cm to microns - xPix = kconv*iHit->GetXL(); // Geant cm to micron - yPix = kconv*iHit->GetYL(); // Geant cm to micron - - // Get track status - status = iHit->GetTrackStatus(); - if(status != 66 && status != 68) { - printf("!!!!! no order status %d\n",status); - } - - 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; - } - // enter Si only - if (layer == 2 && status == 66 && yPix < -71.) { - yPix0 = yPix; - } - depEnergy = iHit->GetIonization(); - // skip if the input point to Si - if(depEnergy <= 0.) 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 ------- - - charge = depEnergy*kEntoEl; // charge in e- - drPath = 0.; - tAng = 0.; - 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 - - - projDif = sqrt(xDif*xDif + zDif*zDif); - ndZ = (Int_t)TMath::Abs(zDif/zPitch) + 1; - ndX = (Int_t)TMath::Abs(xDif/xPitch) + 1; - - // number of the steps along the track: - 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; - dCharge = charge/nsteps; // charge in e- for one step - dZ = zDif/nsteps; - 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++) { - dZn = iZi*dZ; - dXn = iZi*dX; - zPixn = zPix0 + dZn; - xPixn = xPix0 + dXn; - - if(TMath::Abs(projDif) >= 5.) { - 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) { - 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 + spdLenght/2.); - xPixn = (xPixn + spdWidth/2.); - fSegmentation->GetCellIxz(xPixn,zPixn,nXpix,nZpix); - zPitch = fSegmentation->Dpz(nZpix); - // set the window for the integration - jzmin = 1; - jzmax = 3; - if(nZpix == 1) jzmin =2; - if(nZpix == fNPixelsZ) jzmax = 2; - - jxmin = 1; - jxmax = 3; - if(nXpix == 1) jxmin =2; - if(nXpix == fNPixelsX) jxmax = 2; - - zPix = nZpix; - dZright = zPitch*(zPix - zPixn); - dZleft = zPitch - dZright; - xPix = nXpix; - dXright = xPitch*(xPix - xPixn); - dXleft = xPitch - dXright; - dZprev = 0.; - dZnext = 0.; - dXprev = 0.; - 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; - } - // lz changes from 1 to the fNofPixels(270) - lz = nZpix + jz -2; - - zArg1 = dZprev/sigmaDif; - zArg2 = dZnext/sigmaDif; - zProb1 = TMath::Erfc(zArg1); - zProb2 = TMath::Erfc(zArg2); - dZCharge =0.5*(zProb1-zProb2)*dCharge; - - //printf("dZCharge %f \n",dZCharge); - - // ----------- 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; - } - lx = nXpix + jx -2; - - xArg1 = dXprev/sigmaDif; - xArg2 = dXnext/sigmaDif; - xProb1 = TMath::Erfc(xArg1); - xProb2 = TMath::Erfc(xArg2); - dXCharge =0.5*(xProb1-xProb2)*dZCharge; - - //printf("dXCharge %f \n",dXCharge); - - if(dXCharge > 1.) { - index = lz-1; - if (first) { - indexRange[0]=indexRange[1]=index; - indexRange[2]=indexRange[3]=lx-1; - first=kFALSE; - } - - indexRange[0]=TMath::Min(indexRange[0],lz-1); - indexRange[1]=TMath::Max(indexRange[1],lz-1); - indexRange[2]=TMath::Min(indexRange[2],lx-1); - indexRange[3]=TMath::Max(indexRange[3],lx-1); - // build the list of digits for this module - signal=fMapA2->GetSignal(index,lx-1); - signal+=dXCharge; - fMapA2->SetHit(index,lx-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; //ch - - if (lasttrack != idtrack || hit==(nhits-1)) { - GetList(idtrack,pList,indexRange); - first=kTRUE; - } - lasttrack=idtrack; - } // 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,Float_t **pList,Int_t *indexRange) { - // loop over nonzero digits - - Int_t ix,iz,globalIndex; - Float_t signal; - Float_t highest,middle,lowest; - - for(iz=indexRange[0];izGetSignal(iz,ix); - - globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 0! - if(!pList[globalIndex]){ - - // - // Create new list (6 elements - 3 signals and 3 tracks + total sig) - // - - pList[globalIndex] = new Float_t [6]; - - // set list to -2 - - *pList[globalIndex] = -2.; - *(pList[globalIndex]+1) = -2.; - *(pList[globalIndex]+2) = -2.; - *(pList[globalIndex]+3) = 0.; - *(pList[globalIndex]+4) = 0.; - *(pList[globalIndex]+5) = 0.; - - - *pList[globalIndex] = (float)label; - *(pList[globalIndex]+3) = signal; - } - else{ - - // check the signal magnitude - - highest = *(pList[globalIndex]+3); - middle = *(pList[globalIndex]+4); - 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; - } - else if (signal>middle){ - *(pList[globalIndex]+5) = middle; - *(pList[globalIndex]+4) = signal; - - *(pList[globalIndex]+2) = *(pList[globalIndex]+1); - *(pList[globalIndex]+1) = label; - } - else{ - *(pList[globalIndex]+5) = signal; - *(pList[globalIndex]+2) = label; - } +} +//______________________________________________________________________ +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; + 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::ChargeToSignal(Float_t **pList) { - // charge to signal - - AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS"); - - - TRandom *random = new TRandom(); - Float_t threshold = (float)fResponse->MinVal(); - - Int_t digits[3], tracks[3],gi,j1; - Float_t charges[3]; - Float_t electronics; - Float_t signal,phys; - Int_t iz,ix; - for(iz=0;izGaus(); - signal = (float)fMapA2->GetSignal(iz,ix); - signal += electronics; - if (signal > threshold) { - digits[0]=iz; - digits[1]=ix; - digits[2]=1; - gi =iz*fNPixelsX+ix; // global index - for(j1=0;j1<3;j1++){ - tracks[j1] = (Int_t)(*(pList[gi]+j1)); - charges[j1] = 0; - } - phys=0; - aliITS->AddDigit(0,phys,digits,tracks,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()); + 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::CreateHistograms() { - // CreateHistograms - - Int_t i; - for(i=0;iAppend(candnum); - (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps", - fNPixelsX,0.,(Float_t) fNPixelsX); - delete spdname; - } - +//______________________________________________________________________ +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() { - // +//______________________________________________________________________ +void AliITSsimulationSPD::ResetHistograms(){ // Reset histograms for this detector - // - Int_t i; - for(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
 }