/* questions to experts: why RemoveDeadPixels should be called before FrompListToDigits ? */ /************************************************************************** * 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. * **************************************************************************/ #include #include #include #include #include "AliITSU.h" #include "AliITSUDigitPix.h" #include "AliITSUHit.h" #include "AliITSUModule.h" #include "AliITSUSensMap.h" #include "AliITSUCalibrationPix.h" #include "AliITSUSegmentationPix.h" #include "AliITSUSimulationPix.h" #include "AliLog.h" #include "AliRun.h" #include "AliMagF.h" #include "AliMathBase.h" #include "AliITSUSimuParam.h" #include "AliITSUSDigit.h" ClassImp(AliITSUSimulationPix) //////////////////////////////////////////////////////////////////////// // 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 // // Version: 0 // Written by Boris Batyunya // December 20 1999 // // Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch // // AliITSUSimulationPix is to do the simulation of pixels // //////////////////////////////////////////////////////////////////////// //______________________________________________________________________ AliITSUSimulationPix::AliITSUSimulationPix() : fTanLorAng(0) ,fStrobe(kTRUE) ,fStrobeLenght(4) ,fStrobePhase(-12.5e-9) { // Default constructor. } //______________________________________________________________________ AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map) :AliITSUSimulation(sim,map) ,fTanLorAng(0) ,fStrobe(kTRUE) ,fStrobeLenght(4) ,fStrobePhase(-12.5e-9) { // standard constructor Init(); } //______________________________________________________________________ AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s) :AliITSUSimulation(s) ,fTanLorAng(s.fTanLorAng) ,fStrobe(s.fStrobe) ,fStrobeLenght(s.fStrobeLenght) ,fStrobePhase(s.fStrobePhase) { // Copy Constructor } //______________________________________________________________________ AliITSUSimulationPix::~AliITSUSimulationPix() { // destructor // only the sens map is owned and it is deleted by ~AliITSUSimulation } //______________________________________________________________________ AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s) { // Assignment operator if (&s == this) return *this; AliITSUSimulation::operator=(s); fStrobe = s.fStrobe; fStrobeLenght = s.fStrobeLenght; fStrobePhase = s.fStrobePhase; return *this; } //______________________________________________________________________ void AliITSUSimulationPix::Init() { // Initilization if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight()); } //______________________________________________________________________ Bool_t AliITSUSimulationPix::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 (weightHole<0) { weightHole=0.; AliWarning("You have asked for negative Hole weight"); AliWarning("I'm going to use only electrons"); } if (weightHole>1) { weightHole=1.; AliWarning("You have asked for weight > 1"); AliWarning("I'm going to use only holes"); } Double_t weightEle=1.-weightHole; AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); if (!fld) AliFatal("The field is not initialized"); Double_t bz = fld->SolenoidField(); fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) + weightEle*fSimuParam->LorentzAngleElectron(bz)); return kTRUE; } //_____________________________________________________________________ void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg) { // This function begins the work of creating S-Digits. if (!(mod->GetNHits())) { AliDebug(1,Form("In event %d module %d there are %d hits returning.", event, mod->GetIndex(),mod->GetNHits())); return; } // if (fStrobe) if (event != fEvent) GenerateStrobePhase(); InitSimulationModule(mod->GetIndex(), event, seg ); // Hits2SDigits(mod); Hits2SDigitsFast(mod); if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels(); if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels(); WriteSDigits(); ClearMap(); } //______________________________________________________________________ void AliITSUSimulationPix::WriteSDigits() { // This function adds each S-Digit to pList static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS"); int nsd = fSensMap->GetEntries(); for (int i=0;iAt(i); // ordered in index if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue; aliITS->AddSumDigit(*sd); } return; } //______________________________________________________________________ void AliITSUSimulationPix::FinishSDigitiseModule() { // This function calls SDigitsToDigits which creates Digits from SDigits FrompListToDigits(); ClearMap(); return; } //______________________________________________________________________ void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg) { // 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 Hits2SDigits // if (fStrobe) if (event != fEvent) GenerateStrobePhase(); InitSimulationModule( mod->GetIndex(), event, seg ); // Hits2SDigits(mod); Hits2SDigitsFast(mod); // if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels(); if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels(); FrompListToDigits(); ClearMap(); } //______________________________________________________________________ void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod) { // Does the charge distributions using Gaussian diffusion charge charing. const Double_t kmictocm = 1.0e-4; // convert microns to cm. const Double_t kcmtomic = 1.e4; const Double_t kBunchLenght = 25e-9; // LHC clock Int_t nhits = mod->GetNHits(); if (!nhits) return; // Int_t h,ix,iz,i; Int_t idtrack; Float_t x,y,z; // keep coordinates float (required by AliSegmentation) 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 t,tp,st,dt=0.2,el,sig,sigx,sigz,fda; Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half Thickness fSimuParam->GetPixSigmaDiffusionAsymmetry(fda); // for (h=0;hGetHit(h)->GetTOF()GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght))) ) continue; // if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); if (st>0.0) { st = (Double_t)((Int_t)(st*kcmtomic)); // 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 (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside el = dt * de / fSimuParam->GetGeVToCharge(); // if (GetDebug(1)) if (el<=0.0) cout<<"el="<GetPixCouplingOption())); } //______________________________________________________________________ void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod) { // Does the charge distributions using Gaussian diffusion charge charing. // Inputs: // AliITSUModule *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}; const Double_t kBunchLenght = 25e-9; // LHC clock TObjArray *hits = mod->GetHits(); Int_t nhits = hits->GetEntriesFast(); if (nhits<=0) return; // Int_t h,ix,iz,i; Int_t idtrack; Float_t x,y,z; // keep coordinates float (required by AliSegmentation) Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0; Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half thickness fSimuParam->GetPixSigmaDiffusionAsymmetry(fda); // for (h=0;hGetHit(h)->GetTOF()GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght))) ) continue; // if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); if (st>0.0) for (i=0;iLocalToDet(x,z,ix,iz))) continue; // outside el = kwi[i]*de/fSimuParam->GetGeVToCharge(); sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); sigx=sig; sigz=sig*fda; if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng; SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); // cout << "sigx sigz " << sigx << " " << sigz << endl; // dom } // end for i // End Integrate over t else { // st == 0.0 deposit it at this point x = x0; y = y0; z = z0; if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside el = de / fSimuParam->GetGeVToCharge(); sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); sigx=sig; sigz=sig*fda; if (fSimuParam->GetPixLorentzDrift()) 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 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster AliITSUSDigit* dg = 0; switch (fSimuParam->GetPixCouplingOption()) { case AliITSUSimuParam::kNewCouplingPix : for (i=nd;i--;) { dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); if (fSensMap->IsDisabled(dg)) continue; SetCoupling(dg,idtrack,h); } case AliITSUSimuParam::kOldCouplingPix: for (i=nd;i--;) { dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); if (fSensMap->IsDisabled(dg)) continue; SetCouplingOld(dg,idtrack,h); } break; default: break; } // end switch if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption())); } //______________________________________________________________________ void AliITSUSimulationPix::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 fSimuParam->GetPixLorentzDrift()=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; // keep coordinates float (required by AliSegmentation) Double_t s,sp,x1,x2,z1,z2; // if (GetDebug(2)) AliInfo(Form("(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. UpdateMapSignal(iz0,ix0,t,hi,el); return; } // end if sp = 1.0/(sig*kRoot2); ixs = TMath::Max(-knx+ix0,0); ixe = TMath::Min(knx+ix0,fSeg->Npx()-1); izs = TMath::Max(-knz+iz0,0); ize = TMath::Min(knz+iz0,fSeg->Npz()-1); for (ix=ixs;ix<=ixe;ix++) for (iz=izs;iz<=ize;iz++) { fSeg->DetToLocal(ix,iz,x,z); // pixel center double dxi = 0.5*kmictocm*fSeg->Dpx(ix); double dzi = 0.5*kmictocm*fSeg->Dpz(iz); x1 = x; z1 = z; x2 = x1 + dxi; // Upper x1 -= dxi; // Lower z2 = z1 + dzi; // Upper z1 -= dzi; // 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 = el*0.25; // Correction based on definision of Erfc s *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2); if (GetDebug(3)) { cout <<"el="<GetMaxIndex(); double minProb = 0.1/maxInd; // int nsd = fSensMap->GetEntries(); Int_t prevID=0,curID=0; TArrayI ordSampleInd(100),ordSample(100); // double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule); fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean); probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold // for (int i=0;iAt(i); // ordered in index if (fSensMap->IsDisabled(sd)) continue; curID = (int)sd->GetUniqueID(); // if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean); prevID = curID+1; } // if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue; if (TMath::Abs(sig)>2147483647.0) { //RS? //PH 2147483647 is the max. integer //PH This apparently is a problem which needs investigation AliWarning(Form("Too big or too small signal value %f",sig)); } fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix); dig.SetCoord1(iz); dig.SetCoord2(ix); dig.SetSignal(1); dig.SetSignalPix((Int_t)sig); int ntr = sd->GetNTracks(); for (int j=0;jGetTrack(j)); dig.SetHit(j,sd->GetHit(j)); } for (int j=ntr;jAddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig); } // if needed, add noisy pixels with id from last real hit to maxID if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean); // } //______________________________________________________________________ Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base) { // create random noisy digits above threshold within id range [minID,maxID[ // see FrompListToDigits for details // static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS"); UInt_t ix,iz; static AliITSUDigitPix dig; static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!! const Int_t knmaxtrk=AliITSdigit::GetNTracks(); // Int_t ncand = 0; int npix = maxID-minID; if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd); int* ordV = ordSample.GetArray(); int* ordI = ordSampleInd.GetArray(); for (int j=0;jGetMapIndex((UInt_t)ordV[ordI[j]],iz,ix); // create noisy digit dig.SetCoord1(iz); dig.SetCoord2(ix); dig.SetSignal(1); dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise)); for (int k=knmaxtrk;k--;) { dig.SetTrack(k,-3); dig.SetHit(k,-1); } aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig); if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix())); } return ncand; } //______________________________________________________________________ void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, 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:
  // old                  existing AliITSUSDigit
  // Int_t ntrack         track incex number
  // Int_t idhit          hit index number
  UInt_t col,row;
  Double_t pulse1,pulse2;
  Double_t couplR=0.0,couplC=0.0;
  Double_t xr=0.;
  //
  fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
  fSimuParam->GetPixCouplingParam(couplC,couplR);
  if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
				col,row,ntrack,idhit,couplC,couplR));
  pulse2 = pulse1 = old->GetSignal();
  if (pulse1GetPixMinElToAdd()) return; // too small signal
  for (Int_t isign=-1;isign<=1;isign+=2) {
    //
    // loop in col direction
    int j1 = int(col) + isign;
    xr = gRandom->Rndm();
    if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
    //
    // loop in row direction
    int j2 = int(row) + isign;
    xr = gRandom->Rndm();
    if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
  } 
  //
}

//______________________________________________________________________
void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, 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:
  // old            existing AliITSUSDigit
  // ntrack         track incex number
  // idhit          hit index number
  // module         module number
  //
  UInt_t col,row;
  Double_t pulse1,pulse2;
  Double_t couplR=0.0,couplC=0.0;
  //
  fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
  fSimuParam->GetPixCouplingParam(couplC,couplR);
  if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
				col,row,ntrack,idhit,couplC,couplR));
 //
 if (old->GetSignal()GetPixMinElToAdd()) return; // too small signal
 for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
   pulse2 = pulse1 = old->GetSignal();
   //
   int j1 = int(col)+isign;
   pulse1 *= couplC;    
   if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1GetPixThreshold(fModule))) pulse1 = old->GetSignal();
   else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
   
   // loop in row direction
   int j2 = int(row) + isign;
   pulse2 *= couplR;
   if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2GetPixThreshold(fModule))) pulse2 = old->GetSignal();
   else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 } // for isign
}

//______________________________________________________________________
void AliITSUSimulationPix::GenerateStrobePhase()
{
  // Generate randomly the strobe
  // phase w.r.t to the LHC clock
  // Done once per event
  const Double_t kBunchLenght = 25e-9; // LHC clock
  fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
    (Double_t)fStrobeLenght*kBunchLenght+
    kBunchLenght/2;
}