2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 #include <TGeoGlobalMagField.h>
22 #include "AliITSUDigitPix.h"
23 #include "AliITSUHit.h"
24 #include "AliITSUChip.h"
25 #include "AliITSUSensMap.h"
26 #include "AliITSUCalibrationPix.h"
27 #include "AliITSUSegmentationPix.h"
28 #include "AliITSUSimulationPix.h"
32 #include "AliMathBase.h"
33 #include "AliITSUSimuParam.h"
34 #include "AliITSUSDigit.h"
35 #include "AliITSUParamList.h"
39 using namespace TMath;
41 ClassImp(AliITSUSimulationPix)
42 ////////////////////////////////////////////////////////////////////////
44 // Modified by D. Elia, G.E. Bruno, H. Tydesjo
45 // Fast diffusion code by Bjorn S. Nilsen
47 // October 2007: GetCalibrationObjects() removed
50 // Written by Boris Batyunya
53 // Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
55 // AliITSUSimulationPix is to do the simulation of pixels
57 // 2013 Feb: Added MonoPix response and nois calculation al la MIMOSA32 (levente.molnar@cern.ch)
59 ////////////////////////////////////////////////////////////////////////
61 //______________________________________________________________________
62 AliITSUSimulationPix::AliITSUSimulationPix()
64 ,fGlobalChargeScale(1.0)
69 // Default constructor.
70 SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
73 //______________________________________________________________________
74 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
75 :AliITSUSimulation(sim,map)
77 ,fGlobalChargeScale(1.0)
82 // standard constructor
83 SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
87 //______________________________________________________________________
88 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
90 ,fTanLorAng(s.fTanLorAng)
91 ,fGlobalChargeScale(s.fGlobalChargeScale)
92 ,fSpread2DHisto(s.fSpread2DHisto)
93 ,fSpreadFun(s.fSpreadFun)
94 ,fROTimeFun(s.fROTimeFun)
100 //______________________________________________________________________
101 AliITSUSimulationPix::~AliITSUSimulationPix()
104 // only the sens map is owned and it is deleted by ~AliITSUSimulation
107 //______________________________________________________________________
108 AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
110 // Assignment operator
111 if (&s == this) return *this;
112 AliITSUSimulation::operator=(s);
113 fSpread2DHisto = s.fSpread2DHisto;
115 fGlobalChargeScale = s.fGlobalChargeScale;
116 fSpreadFun = s.fSpreadFun;
117 fROTimeFun = s.fROTimeFun;
122 //______________________________________________________________________
123 void AliITSUSimulationPix::Init()
126 if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
129 //______________________________________________________________________
130 Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
132 // This function set the Tangent of the Lorentz angle.
133 // A weighted average is used for electrons and holes
134 // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
135 // output: Bool_t : kTRUE in case of success
139 AliWarning("You have asked for negative Hole weight");
140 AliWarning("I'm going to use only electrons");
144 AliWarning("You have asked for weight > 1");
145 AliWarning("I'm going to use only holes");
147 Double_t weightEle=1.-weightHole;
148 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
149 if (!fld) AliFatal("The field is not initialized");
150 Double_t bz = fld->SolenoidField();
151 fTanLorAng = Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
152 weightEle*fSimuParam->LorentzAngleElectron(bz));
156 //_____________________________________________________________________
157 void AliITSUSimulationPix::SDigitiseChip()
159 // This function begins the work of creating S-Digits.
161 AliDebug(10,Form("In event %d chip %d there are %d hits", fEvent, fChip->GetIndex(),fChip->GetNHits()));
162 if (fChip->GetNHits()) {
163 if(fResponseParam->GetParameter(kDigitalSim) == 0 ) Hits2SDigitsFast(); // analogue chip response simulation
164 else Hits2SDigitsFastDigital(); // digital chip response
166 if (!fSensMap->GetEntries()) return;
172 //______________________________________________________________________
173 void AliITSUSimulationPix::WriteSDigits()
175 // This function adds each S-Digit to pList
176 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
177 int nsd = fSensMap->GetEntries();
180 for (int i=0;i<nsd;i++) {
181 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
182 if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
183 aliITS->AddSumDigit(*sd);
188 //______________________________________________________________________
189 void AliITSUSimulationPix::FinishSDigitiseChip()
191 // This function calls SDigitsToDigits which creates Digits from SDigits
197 //______________________________________________________________________
198 void AliITSUSimulationPix::DigitiseChip()
200 // This function creates Digits straight from the hits and then adds
201 // electronic noise to the digits before adding them to pList
202 // Each of the input variables is passed along to Hits2SDigits
204 if(fResponseParam->GetParameter(kDigitalSim) == 0 ) Hits2SDigitsFast(); // analogue chip response simulation
205 else Hits2SDigitsFastDigital(); // digital chip response
206 FinishSDigitiseChip();
209 //______________________________________________________________________
210 void AliITSUSimulationPix::Hits2SDigits()
212 // Does the charge distributions using Gaussian diffusion charge charing.
213 Int_t nhits = fChip->GetNHits();
218 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
219 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
220 Double_t t,tp,st,dt=0.2,el;
221 Double_t thick = 0.5*fSeg->Dy(); // Half Thickness
224 for (h=0;h<nhits;h++) {
226 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
227 st = Sqrt(x1*x1+y1*y1+z1*z1);
229 st = (Double_t)((Int_t)(st*1e4)); // number of microns
230 if (st<=1.0) st = 1.0;
231 dt = 1.0/st; // RS TODO: do we need 1 micron steps?
232 double dy = dt*thick;
234 for (t=0.0;t<1.0;t+=dt) { // Integrate over t
239 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
240 el = fGlobalChargeScale * dt * de / fSimuParam->GetGeVToCharge();
242 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
243 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
245 } else { // st == 0.0 deposit it at this point
249 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
250 el = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
251 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
252 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
254 } // Loop over all hits h
257 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
258 AliITSUSDigit* dg = 0;
259 switch (fSimuParam->GetPixCouplingOption()) {
260 case AliITSUSimuParam::kNoCouplingPix :
262 case AliITSUSimuParam::kNewCouplingPix :
264 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
265 if (fSensMap->IsDisabled(dg)) continue;
269 case AliITSUSimuParam::kOldCouplingPix:
271 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
272 if (fSensMap->IsDisabled(dg)) continue;
280 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
283 //______________________________________________________________________
284 void AliITSUSimulationPix::Hits2SDigitsFast()
286 // Does the charge distributions using Gaussian diffusion charge charing. // Inputs:
287 // AliITSUChip *mod Pointer to this chip
289 TObjArray *hits = fChip->GetHits();
290 Int_t nhits = hits->GetEntriesFast();
291 if (nhits<=0) return;
295 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
296 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
297 Double_t step,st,el,de=0.0;
298 Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
299 Double_t thick = fSeg->Dy();
301 for (h=0;h<nhits;h++) {
303 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
305 st = Sqrt(x1*x1+y1*y1+z1*z1);
307 int np = int(1.5*st/minDim); //RStmp: inject the points in such a way that there is ~1.5 point per cell
308 np = TMath::Max(1.0*np,fResponseParam->GetParameter(kSpreadFunMinSteps));
309 AliDebug(10,Form(" Number of charge injection steps is set to %d ",np));
310 double dstep = 1./np;
311 double dy = dstep*thick;
314 for (i=0;i<np;i++) { //RStmp Integrate over t
315 // for (i=0;i<kn10;i++) { // Integrate over t
316 step += dstep; // RStmp kti[i];
320 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
321 el = fGlobalChargeScale * dstep * de/fSimuParam->GetGeVToCharge();
322 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
323 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
324 } // end for i // End Integrate over t
326 else { // st == 0.0 deposit it at this point
330 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
331 el = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
332 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
333 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
336 } // Loop over all hits h
339 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
340 AliITSUSDigit* dg = 0;
341 switch (fSimuParam->GetPixCouplingOption()) {
342 case AliITSUSimuParam::kNoCouplingPix :
344 case AliITSUSimuParam::kNewCouplingPix :
346 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
347 if (fSensMap->IsDisabled(dg)) continue;
350 case AliITSUSimuParam::kOldCouplingPix:
352 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
353 if (fSensMap->IsDisabled(dg)) continue;
360 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
363 //______________________________________________________________________
364 void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
365 Double_t el, Double_t tof, Int_t tID, Int_t hID)
367 // Spreads the charge over neighboring cells. Assume charge is distributed
368 // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
369 // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
370 // Defined this way, the integral over all x and z is el.
372 // Double_t x0 x position of point where charge is liberated (local)
373 // Double_t z0 z position of point where charge is liberated (local)
374 // Double_t dy distance from the entrance surface (diffusion sigma may depend on it)
375 // Int_t ix0 row of cell corresponding to point x0
376 // Int_t iz0 columb of cell corresponding to point z0
377 // Double_t el number of electrons liberated in this step
378 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
379 // Double_t sigz Sigma difusion along z for this step (z0 dependent)
380 // Int_t tID track number
381 // Int_t hID hit "hit" index number
383 Int_t ix,iz,ixs,ixe,izs,ize;
384 Float_t x,z; // keep coordinates float (required by AliSegmentation)
385 Float_t xdioshift = 0 , zdioshift = 0 ;
386 Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
388 if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,dy=%e, ix0=%d,iz0=%d,el=%e,tID=%d,hID=%d)",x0,z0,dy,ix0,iz0,el,tID,hID));
390 Double_t &x1 = dtIn[kCellX1];
391 Double_t &x2 = dtIn[kCellX2];
392 Double_t &z1 = dtIn[kCellZ1];
393 Double_t &z2 = dtIn[kCellZ2];
395 int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
396 int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
398 dtIn[kCellYDepth] = dy;
399 ixs = Max(-nx+ix0,0);
400 ixe = Min( nx+ix0,fSeg->Npx()-1);
401 izs = Max(-nz+iz0,0);
402 ize = Min( nz+iz0,fSeg->Npz()-1);
403 for (ix=ixs;ix<=ixe;ix++)
404 for (iz=izs;iz<=ize;iz++) {
406 // Check if the hit is inside readout window
407 int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
408 if (Abs(cycleRO)>kMaxROCycleAccept) continue;
410 fSeg->DetToLocal(ix,iz,x,z); // pixel center
411 xdioshift = zdioshift = 0;
412 double dxi = fSeg->Dpx(ix);
413 double dzi = fSeg->Dpz(iz);
414 CalcDiodeShiftInPixel(ix,iz,xdioshift,zdioshift); // Check and apply diode shift if needed
419 // printf("DShift: %d %d -> %.4f %.4f\n",ix,iz,xdioshift,zdioshift);
420 x1 = (x + xdioshift) - x0; // calculate distance of cell boundaries from injection center
421 z1 = (z + zdioshift) - z0;
422 x2 = x1 + dxi; // Upper
424 z2 = z1 + dzi; // Upper
426 s = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
427 if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s,cycleRO);
432 //______________________________________________________________________
433 Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
435 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
436 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
437 // The spread function is assumed to be double gaussian in 2D
438 // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
441 double intg1 = GausInt2D(fResponseParam->GetParameter(kG2SigX0), // sigX
442 dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX0), // x1-xmean
443 dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX0), // x2-xmean
444 fResponseParam->GetParameter(kG2SigZ0), // sigZ
445 dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ0), // z1-zmean
446 dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ0)); // z2-zmean
448 double intg2 = GausInt2D(fResponseParam->GetParameter(kG2SigX1), // sigX
449 dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX1), // x1-xmean
450 dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX1), // x2-xmean
451 fResponseParam->GetParameter(kG2SigZ1), // sigZ
452 dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ1), // z1-zmean
453 dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ1)); // z2-zmean
454 double scl = fResponseParam->GetParameter(kG2ScaleG2);
455 return (intg1+intg2*scl)/(1+scl);
460 //______________________________________________________________________
461 Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
463 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
464 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
465 // The spread function integral is taken from fSpread2DHisto extracted from the sensor response parameters
466 // list in the method SetResponseParam. The histo must return the fraction of charge integrates in the
467 // cell whose center is on the distance X=(dtIn[kCellX1]+dtIn[kCellX2])/2 and Z=(dtIn[kCellZ1]+dtIn[kCellZ2])/2
468 // from the injection point.
470 Double_t qpixfrac = 0;
471 Double_t xintp = 1e4*(dtIn[kCellX1]+dtIn[kCellX2])/2.0;
472 Double_t zintp = 1e4*(dtIn[kCellZ1]+dtIn[kCellZ2])/2.0;
474 qpixfrac = fSpread2DHisto->Interpolate(xintp,zintp); //the PSF map is given in um but the dtIn is in cm so we need to convert it
479 //______________________________________________________________________
480 Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
482 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
483 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
484 // The spread function is assumed to be gaussian in 2D
485 // Parameters should be: mean0,sigma0
486 return GausInt2D(fResponseParam->GetParameter(kG1SigX), // sigX
487 dtIn[kCellX1]-fResponseParam->GetParameter(kG1MeanX), // x1-xmean
488 dtIn[kCellX2]-fResponseParam->GetParameter(kG1MeanX), // x2-xmean
489 fResponseParam->GetParameter(kG1SigZ), // sigZ
490 dtIn[kCellZ1]-fResponseParam->GetParameter(kG1MeanZ), // z1-zmean
491 dtIn[kCellZ2]-fResponseParam->GetParameter(kG1MeanZ)); // z2-zmean
495 //______________________________________________________________________
496 void AliITSUSimulationPix::RemoveDeadPixels()
498 // Removes dead pixels on each chip (ladder)
499 // This should be called before going from sdigits to digits (i.e. from FrompListToDigits)
501 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
504 if (calObj->IsBad()) {ClearMap(); return;} // whole chip is masked
506 // prepare the list of r/o cycles seen
507 Char_t cyclesSeen[2*kMaxROCycleAccept+1];
509 for (int i=(2*kMaxROCycleAccept+1);i--;) if (fCyclesID[i]) cyclesSeen[ncycles++]=i-kMaxROCycleAccept;
511 // remove single bad pixels one by one
512 int nsingle = calObj->GetNrBadSingle();
515 for (int i=nsingle;i--;) {
516 calObj->GetBadPixelSingle(i,row,col);
517 for (int icl=ncycles;icl--;) fSensMap->DeleteItem(col,row,cyclesSeen[icl]);
519 int nsd = fSensMap->GetEntriesUnsorted();
520 for (int isd=nsd;isd--;) {
521 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
522 if (fSensMap->IsDisabled(sd)) continue;
523 fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,cycle);
524 int chip = fSeg->GetChipFromChannel(0,col);
525 // if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
526 if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
531 //______________________________________________________________________
532 void AliITSUSimulationPix::AddNoisyPixels()
534 // Adds noisy pixels on each chip (ladder)
535 // This should be called before going from sdigits to digits (i.e. FrompListToDigits)
536 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
537 if (!calObj) { AliDebug(10,Form(" No Calib Object for Noise!!! ")); return;}
538 for (Int_t i=calObj->GetNrBad(); i--;)
540 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 )
541 UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),10*fSimuParam->GetPixThreshold(fChip->GetIndex()));
543 UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),kNoisyPixOCDB );
548 //______________________________________________________________________
549 void AliITSUSimulationPix::FrompListToDigits()
551 // add noise and electronics, perform the zero suppression and add the digits to the list
553 // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
555 int nsd = fSensMap->GetEntriesUnsorted(); // sdigits added from the signal
557 // add different kinds of noise.
558 Bool_t addNoisy = fSimuParam->GetPixAddNoisyFlag() && (nsd>0 || fSimuParam->GetPixNoiseInAllMod()); // do we generate noise?
560 AddNoisyPixels(); // constantly noisy channels
561 AddRandomNoisePixels(0.0); // random noise: at the moment generate noise only for instance 0
562 nsd = fSensMap->GetEntriesUnsorted();
565 if (nsd && fSimuParam->GetPixRemoveDeadFlag()) {
567 // note that here we shall use GetEntries instead of GetEntriesUnsorted since the
568 // later operates on the array where the elements are not removed by flagged
569 nsd = fSensMap->GetEntries();
571 if (!nsd) return; // nothing to digitize
574 Int_t iCycle,modId = fChip->GetIndex();
576 const Int_t knmaxtrk=AliITSdigit::GetNTracks();
577 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
578 static AliITSUDigitPix dig;
580 for (int i=0;i<nsd;i++) {
581 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
582 if (fSensMap->IsDisabled(sd)) continue;
584 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 &&
585 (sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(modId)) continue; //Threshold only applies in analogue simulation
587 if (Abs(sig)>2147483647.0) { //RS?
588 //PH 2147483647 is the max. integer
589 //PH This apparently is a problem which needs investigation
590 AliWarning(Form("Too big or too small signal value %f",sig));
592 fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,iCycle);
595 dig.SetROCycle(iCycle);
596 dig.SetSignal((Int_t)sig);
597 dig.SetSignalPix((Int_t)sig);
598 int ntr = sd->GetNTracks();
599 for (int j=0;j<ntr;j++) {
600 dig.SetTrack(j,sd->GetTrack(j));
601 dig.SetHit(j,sd->GetHit(j));
603 for (int j=ntr;j<knmaxtrk;j++) {
607 aliITS->AddSimDigit(AliITSUGeomTGeo::kChipTypePix, &dig);
612 //______________________________________________________________________
613 Int_t AliITSUSimulationPix::AddRandomNoisePixels(Double_t tof)
615 // create random noisy sdigits above threshold
617 int modId = fChip->GetIndex();
618 int npix = fSeg->GetNPads();
619 int ncand = gRandom->Poisson( npix*fSimuParam->GetPixFakeRate() );
620 if (ncand<1) return 0;
622 double probNoisy,noiseSig,noiseMean,thresh;
626 static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
627 ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
628 int* ordV = ordSample.GetArray();
629 int* ordI = ordSampleInd.GetArray();
631 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 ) {
632 thresh = fSimuParam->GetPixThreshold(modId);
633 fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
634 probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
636 for (int j=0;j<ncand;j++) {
637 fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],col,row,iCycle); // create noisy digit
638 iCycle = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(row,col,tof);
639 UpdateMapNoise(col,row,AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,noiseMean,noiseSig), iCycle);
645 for (int j=0;j<ncand;j++) {
646 fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],col,row,iCycle); // create noisy digit
647 iCycle = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(row,col,tof);
648 UpdateMapNoise(col,row,kNoisyPixRnd, iCycle);
656 //______________________________________________________________________
657 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old)
659 // Take into account the coupling between adiacent pixels.
660 // The parameters probcol and probrow are the probability of the
661 // signal in one pixel shared in the two adjacent pixels along
662 // the column and row direction, respectively.
663 // Note pList is goten via GetMap() and chip is not need any more.
664 // Otherwise it is identical to that coded by Tiziano Virgili (BSN).
667 Double_t pulse1,pulse2;
668 Double_t couplR=0.0,couplC=0.0;
671 fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
673 fSimuParam->GetPixCouplingParam(couplC,couplR);
674 if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,couplC=%e couplR=%e",
675 col,row,couplC,couplR));
676 pulse2 = pulse1 = old->GetSignal();
677 if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
678 for (Int_t isign=-1;isign<=1;isign+=2) {
680 // loop in col direction
681 int j1 = int(col) + isign;
682 xr = gRandom->Rndm();
683 if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
685 // loop in row direction
686 int j2 = int(row) + isign;
687 xr = gRandom->Rndm();
688 if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
693 //______________________________________________________________________
694 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old)
696 // Take into account the coupling between adiacent pixels.
697 // The parameters probcol and probrow are the fractions of the
698 // signal in one pixel shared in the two adjacent pixels along
699 // the column and row direction, respectively.
701 // old existing AliITSUSDigit
702 // ntrack track incex number
703 // idhit hit index number
708 Int_t modId = fChip->GetIndex();
709 Double_t pulse1,pulse2;
710 Double_t couplR=0.0,couplC=0.0;
712 fSensMap->GetMapIndex(old->GetUniqueID(),col,row,cycle);
713 fSimuParam->GetPixCouplingParam(couplC,couplR);
714 if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,roCycle=%d) couplC=%e couplR=%e",col,row,cycle,couplC,couplR));
716 if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
717 for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
718 pulse2 = pulse1 = old->GetSignal();
720 int j1 = int(col)+isign;
722 if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
723 else UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
725 // loop in row direction
726 int j2 = int(row) + isign;
728 if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
729 else UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
733 //______________________________________________________________________
734 void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
736 // attach response parameterisation data
737 fResponseParam = resp;
739 int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
740 const char* hname = 0;
745 case kSpreadFunHisto:
746 fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
747 hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
748 if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname)))
749 AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
752 case kSpreadFunDoubleGauss2D:
753 fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
756 case kSpreadFunGauss2D:
757 fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
760 default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
763 int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
764 switch (readoutType) {
766 fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
768 case kReadOutRollingShutter:
769 fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
771 default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
774 //___ Set the Rolling Shutter read-out window
775 fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
776 //___ Pixel discrimination threshold, and the S/N cut
777 fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all chips
778 //___ Minimum number of electrons to add
779 fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
780 //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
781 fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all chips
782 //___ Pixel fake hit rate
783 fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
784 //___ To apply the noise or not
785 if ( fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01) fSimuParam->SetPixAddNoisyFlag(kTRUE);
786 else fSimuParam->SetPixAddNoisyFlag(kFALSE);
788 if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
789 else fSimuParam->SetPixNoiseInAllMod(kFALSE);
791 // Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
792 fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
794 AliDebug(10,Form("=============== Setting the response start ============================"));
795 AliDebug(10,Form("=============== Digital (1) / Analogue (0) simu: %f",fResponseParam->GetParameter(kDigitalSim)));
796 AliDebug(10,Form("=============== RO type: %d",readoutType));
797 AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
798 AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
799 AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
800 AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
801 AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
802 AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
803 AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
804 AliDebug(10,Form("=============== Setting the response done ============================"));
808 //______________________________________________________________________
809 Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
812 // Get the read-out cycle of the hit in the given column/row of the sensor.
813 // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
814 // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
815 // GetRollingShutterWindow give the with of the rolling shutter read out window
817 double tmin = fReadOutCycleOffset + fReadOutCycleLength*(double(row)/fSeg->Npx()-1.);
818 int cycle = Nint( (hitTime-tmin)/fReadOutCycleLength - 0.5 );
819 AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : tmax %e -> cycle:%d",row,col,hitTime,tmin,
820 tmin+fReadOutCycleLength,cycle));
825 //______________________________________________________________________
826 Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
829 // Check whether the hit is in the read out window of the given column/row of the sensor
831 AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
832 row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
833 hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
834 return (hitTime<0 || hitTime>fReadOutCycleLength) ? kMaxROCycleAccept+1 : 0;
838 //_______________________________________________________________________
839 void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xrow, Int_t zcol, Float_t &x, Float_t &z)
842 // Calculates the shift of the diode wrt the geometric center of the pixel.
843 // It is needed for staggerred pixel layout or double diode pixels with assymetric center
844 // The shift can depend on the column or line or both...
845 // The x and z are passed in cm
847 ((AliITSUSegmentationPix*)fSeg)->GetDiodShift(xrow,zcol,x,z);
851 //______________________________________________________________________
852 void AliITSUSimulationPix::Hits2SDigitsFastDigital()
854 // Does the digital chip response simulation
855 // AliITSUChip *mod Pointer to this chip
859 TObjArray *hits = fChip->GetHits();
860 Int_t nhits = hits->GetEntriesFast();
861 if (nhits<=0) return;
865 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
866 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
867 Double_t st,el,de=0.0;
870 for (h=0;h<nhits;h++) {
872 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
874 st = Sqrt(x1*x1+y1*y1+z1*z1);
876 //___ place hit to the middle of the path segment - CHANGE LATER !
881 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
883 // Check if the hit is inside readout window
884 int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
885 if (Abs(cycleRO)>kMaxROCycleAccept) continue;
888 el = de / fSimuParam->GetGeVToCharge();
890 PlaceDigitalPixels(x,z,el,tof,idtrack,h);
892 } // Loop over all hits h
895 //______________________________________________________________________
896 void AliITSUSimulationPix::PlaceDigitalPixels(Double_t x0hit,Double_t z0hit, Double_t el, Double_t tof, Int_t tID, Int_t hID)
898 // Place the digital pixel positions on the sensor
900 // Double_t x0hit x position of point where charge is liberated (local) - hit
901 // Double_t z0hit z position of point where charge is liberated (local) - hit
902 // Double_t el number of electrons liberated in this step
903 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
904 // Double_t sigz Sigma difusion along z for this step (z0 dependent)
905 // Int_t tID track number
906 // Int_t hID hit "hit" index number
911 Float_t x,z; // keep coordinates float (required by AliSegmentation)
912 Float_t distX = 0, distZ = 0;
914 //___ TEMPORARY - place a fixed pattern cluster COG to a distance d(x,z) away from hit - TEMPORARY
915 // Pattern used (not realistic ebye but averaging over events yes):
920 //___ This list should come from a look up table based on CluTypeID as well as COG coord
924 distX = rnd.Uniform(-5.0*1e-4, 5.0*1e-4); //in um
925 distZ = rnd.Uniform(-5.0*1e-4, 5.0*1e-4); //in um
930 if(!fSeg->LocalToDet(x,z,ix,iz)) return; // if clu CoG is outside of the chip skipp the cluster -> refine later
932 const Int_t nCluPixels = 5;
933 Int_t aPixListX[nCluPixels] = { 0, -1, 0, 1, 0};
934 Int_t aPixListZ[nCluPixels] = { 1, 0, 0, 0, -1};
936 Double_t s = el / 1.0 / nCluPixels;
940 for (Int_t ipix = 0 ; ipix < nCluPixels; ipix++)
942 nx = ix + aPixListX[ipix];
943 nz = iz + aPixListZ[ipix];
944 cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
945 if ( nx >= 0 && nx <= fSeg -> Npx() && nz >= 0 && nz <= fSeg -> Npz() ) UpdateMapSignal(nz,nx,tID,hID,s,cycleRO); //if the pixel is in the detector