]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ITS/UPGRADE/AliITSUSimulationPix.cxx
Fix for the case of non-existent calibration files
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimulationPix.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
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  **************************************************************************/
16
17 #include <TGeoGlobalMagField.h>
18 #include <TH1.h>
19 #include <TH2.h>
20 #include <TString.h>
21 #include "AliITSU.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"
29 #include "AliLog.h"
30 #include "AliRun.h"
31 #include "AliMagF.h"
32 #include "AliMathBase.h"
33 #include "AliITSUSimuParam.h"
34 #include "AliITSUSDigit.h"
35 #include "AliITSUParamList.h"
36
37 using std::cout;
38 using std::endl;
39 using namespace TMath;
40
41 ClassImp(AliITSUSimulationPix)
42 ////////////////////////////////////////////////////////////////////////
43 //  Version: 1
44 //  Modified by D. Elia, G.E. Bruno, H. Tydesjo
45 //  Fast diffusion code by Bjorn S. Nilsen
46 //  March-April 2006
47 //  October     2007: GetCalibrationObjects() removed
48 //
49 //  Version: 0
50 //  Written by Boris Batyunya
51 //  December 20 1999
52 //
53 //  Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
54 //
55 //  AliITSUSimulationPix is to do the simulation of pixels
56 //
57 //  2013 Feb: Added MonoPix response and nois calculation al la MIMOSA32 (levente.molnar@cern.ch)
58 //
59 ////////////////////////////////////////////////////////////////////////
60
61 //______________________________________________________________________
62 AliITSUSimulationPix::AliITSUSimulationPix()
63 :  fTanLorAng(0)
64 ,fGlobalChargeScale(1.0)
65 ,fSpread2DHisto(0)
66 ,fSpreadFun(0)
67 ,fROTimeFun(0)
68 {
69     // Default constructor.
70     SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
71 }
72
73 //______________________________________________________________________
74 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
75 :AliITSUSimulation(sim,map)
76 ,fTanLorAng(0)
77 ,fGlobalChargeScale(1.0)
78 ,fSpread2DHisto(0)
79 ,fSpreadFun(0)
80 ,fROTimeFun(0)
81 {
82     // standard constructor
83     SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
84     Init();
85 }
86
87 //______________________________________________________________________
88 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
89 :AliITSUSimulation(s)
90 ,fTanLorAng(s.fTanLorAng)
91 ,fGlobalChargeScale(s.fGlobalChargeScale)
92 ,fSpread2DHisto(s.fSpread2DHisto)
93 ,fSpreadFun(s.fSpreadFun)
94 ,fROTimeFun(s.fROTimeFun)
95 {
96     //     Copy Constructor
97 }
98
99
100 //______________________________________________________________________
101 AliITSUSimulationPix::~AliITSUSimulationPix()
102 {
103     // destructor
104     // only the sens map is owned and it is deleted by ~AliITSUSimulation
105 }
106
107 //______________________________________________________________________
108 AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
109 {
110     //    Assignment operator
111     if (&s == this) return *this;
112     AliITSUSimulation::operator=(s);
113     fSpread2DHisto = s.fSpread2DHisto;
114     //
115     fGlobalChargeScale = s.fGlobalChargeScale;
116     fSpreadFun    = s.fSpreadFun;
117     fROTimeFun    = s.fROTimeFun;
118     //
119     return *this;
120 }
121
122 //______________________________________________________________________
123 void AliITSUSimulationPix::Init()
124 {
125     // Initilization
126     if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
127 }
128
129 //______________________________________________________________________
130 Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
131 {
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
136     //
137     if (weightHole<0) {
138         weightHole=0.;
139         AliWarning("You have asked for negative Hole weight");
140         AliWarning("I'm going to use only electrons");
141     }
142     if (weightHole>1) {
143         weightHole=1.;
144         AliWarning("You have asked for weight > 1");
145         AliWarning("I'm going to use only holes");
146     }
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));
153     return kTRUE;
154 }
155
156 //_____________________________________________________________________
157 void AliITSUSimulationPix::SDigitiseChip()
158 {
159     //  This function begins the work of creating S-Digits.
160     
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
165     }
166     if (!fSensMap->GetEntries()) return;
167     WriteSDigits();
168     ClearMap();
169     //
170 }
171
172 //______________________________________________________________________
173 void AliITSUSimulationPix::WriteSDigits()
174 {
175     //  This function adds each S-Digit to pList
176     static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
177     int nsd = fSensMap->GetEntries();
178     
179     
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);
184     }
185     return;
186 }
187
188 //______________________________________________________________________
189 void AliITSUSimulationPix::FinishSDigitiseChip()
190 {
191     //  This function calls SDigitsToDigits which creates Digits from SDigits
192     FrompListToDigits();
193     ClearMap();
194     return;
195 }
196
197 //______________________________________________________________________
198 void AliITSUSimulationPix::DigitiseChip()
199 {
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
203     //
204     if(fResponseParam->GetParameter(kDigitalSim) == 0 ) Hits2SDigitsFast(); // analogue chip response simulation
205     else Hits2SDigitsFastDigital();                                         // digital chip response
206     FinishSDigitiseChip();
207 }
208
209 //______________________________________________________________________
210 void AliITSUSimulationPix::Hits2SDigits()
211 {
212     // Does the charge distributions using Gaussian diffusion charge charing.
213     Int_t nhits = fChip->GetNHits();
214     if (!nhits) return;
215     //
216     Int_t h,ix,iz,i;
217     Int_t idtrack;
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
222     
223     //
224     for (h=0;h<nhits;h++) {
225         //
226         if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
227         st = Sqrt(x1*x1+y1*y1+z1*z1);
228         if (st>0.0) {
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;
233             y = -0.5*dy;
234             for (t=0.0;t<1.0;t+=dt) { // Integrate over t
235                 tp  = t+0.5*dt;
236                 x   = x0+x1*tp;
237                 y  += dy;
238                 z   = z0+z1*tp;
239                 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
240                 el  = fGlobalChargeScale * dt * de / fSimuParam->GetGeVToCharge();
241                 //
242                 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
243                 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
244             } // end for t
245         } else { // st == 0.0 deposit it at this point
246             x   = x0;
247             y   = y0 + 0.5*thick;
248             z   = z0;
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);
253         } // end if st>0.0
254     } // Loop over all hits h
255     //
256     // Coupling
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 :
261             break;
262         case AliITSUSimuParam::kNewCouplingPix :
263             for (i=nd;i--;) {
264                 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
265                 if (fSensMap->IsDisabled(dg)) continue;
266                 SetCoupling(dg);
267             }
268             break;
269         case AliITSUSimuParam::kOldCouplingPix:
270             for (i=nd;i--;) {
271                 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
272                 if (fSensMap->IsDisabled(dg)) continue;
273                 SetCouplingOld(dg);
274             }
275             break;
276         default:
277             break;
278             
279     } // end switch
280     if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
281 }
282
283 //______________________________________________________________________
284 void AliITSUSimulationPix::Hits2SDigitsFast()
285 {
286     // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
287     //    AliITSUChip *mod  Pointer to this chip
288     //
289     TObjArray *hits = fChip->GetHits();
290     Int_t nhits = hits->GetEntriesFast();
291     if (nhits<=0) return;
292     //
293     Int_t h,ix,iz,i;
294     Int_t idtrack;
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();
300     //
301     for (h=0;h<nhits;h++) {
302         //
303         if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
304         //
305         st = Sqrt(x1*x1+y1*y1+z1*z1);
306         if (st>0.0) {
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;
312             y = -0.5*dy;
313             step = -0.5*dstep;
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];
317                 x   = x0+x1*step;
318                 y  += dy;
319                 z   = z0+z1*step;
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
325         }
326         else { // st == 0.0 deposit it at this point
327             x   = x0;
328             y   = y0+0.5*thick;
329             z   = z0;
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);
334         } // end if st>0.0
335         
336     } // Loop over all hits h
337     
338     // Coupling
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 :
343             break;
344         case AliITSUSimuParam::kNewCouplingPix :
345             for (i=nd;i--;) {
346                 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
347                 if (fSensMap->IsDisabled(dg)) continue;
348                 SetCoupling(dg);
349             }
350         case AliITSUSimuParam::kOldCouplingPix:
351             for (i=nd;i--;) {
352                 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
353                 if (fSensMap->IsDisabled(dg)) continue;
354                 SetCouplingOld(dg);
355             }
356             break;
357         default:
358             break;
359     } // end switch
360     if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
361 }
362
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)
366 {
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.
371     // Inputs:
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
382     //
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
387     //
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));
389     //
390     Double_t &x1 = dtIn[kCellX1];
391     Double_t &x2 = dtIn[kCellX2];
392     Double_t &z1 = dtIn[kCellZ1];
393     Double_t &z2 = dtIn[kCellZ2];
394     //
395     int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
396     int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
397     //
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++) {
405             //
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;
409             //
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
415             xdioshift *= dxi;
416             zdioshift *= dzi;
417             dxi *= 0.5;
418             dzi *= 0.5;
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
423             x1 -= dxi;      // Lower
424             z2  = z1 + dzi; // Upper
425             z1 -= dzi;      // Lower
426             s   = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
427             if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s,cycleRO);
428         } // end for ix, iz
429     //
430 }
431
432 //______________________________________________________________________
433 Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
434 {
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
439     //
440     // 1st gaussian
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
447     // 2nd gaussian
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);
456     //
457 }
458
459
460 //______________________________________________________________________
461 Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
462 {
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.
469     //
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;
473     //
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
475     //
476     return qpixfrac;
477 }
478
479 //______________________________________________________________________
480 Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
481 {
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
492     //
493 }
494
495 //______________________________________________________________________
496 void AliITSUSimulationPix::RemoveDeadPixels()
497 {
498     // Removes dead pixels on each chip (ladder)
499     // This should be called before going from sdigits to digits (i.e. from FrompListToDigits)
500     
501     AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
502     if (!calObj) return;
503     //
504     if (calObj->IsBad()) {ClearMap(); return;} // whole chip is masked
505     //
506     // prepare the list of r/o cycles seen
507     Char_t cyclesSeen[2*kMaxROCycleAccept+1];
508     int ncycles = 0;
509     for (int i=(2*kMaxROCycleAccept+1);i--;) if (fCyclesID[i]) cyclesSeen[ncycles++]=i-kMaxROCycleAccept;
510     
511     // remove single bad pixels one by one
512     int nsingle = calObj->GetNrBadSingle();
513     UInt_t col,row;
514     Int_t cycle;
515     for (int i=nsingle;i--;) {
516         calObj->GetBadPixelSingle(i,row,col);
517         for (int icl=ncycles;icl--;) fSensMap->DeleteItem(col,row,cyclesSeen[icl]);
518     }
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
527     }
528     //
529 }
530
531 //______________________________________________________________________
532 void AliITSUSimulationPix::AddNoisyPixels()
533 {
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--;)
539     {
540         if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 )
541             UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),10*fSimuParam->GetPixThreshold(fChip->GetIndex()));
542         else
543             UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),kNoisyPixOCDB );
544     }
545     //
546 }
547
548 //______________________________________________________________________
549 void AliITSUSimulationPix::FrompListToDigits()
550 {
551     // add noise and electronics, perform the zero suppression and add the digits to the list
552     //
553     // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
554     //
555     int nsd = fSensMap->GetEntriesUnsorted();  // sdigits added from the signal
556     //
557     // add different kinds of noise.
558     Bool_t addNoisy = fSimuParam->GetPixAddNoisyFlag() && (nsd>0 || fSimuParam->GetPixNoiseInAllMod()); // do we generate noise?
559     if (addNoisy) {
560         AddNoisyPixels();       // constantly noisy channels
561         AddRandomNoisePixels(0.0); // random noise: at the moment generate noise only for instance 0
562         nsd = fSensMap->GetEntriesUnsorted();
563     }
564     //
565     if (nsd && fSimuParam->GetPixRemoveDeadFlag()) {
566         RemoveDeadPixels();
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();
570     }
571     if (!nsd) return; // nothing to digitize
572     //
573     UInt_t row,col;
574     Int_t iCycle,modId = fChip->GetIndex();
575     Double_t sig;
576     const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
577     static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
578     static AliITSUDigitPix dig;
579     //
580     for (int i=0;i<nsd;i++) {
581         AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
582         if (fSensMap->IsDisabled(sd)) continue;
583         //
584         if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 &&
585             (sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(modId)) continue;   //Threshold only applies in analogue simulation
586         //
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));
591         }
592         fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,iCycle);
593         dig.SetCoord1(col);
594         dig.SetCoord2(row);
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));
602         }
603         for (int j=ntr;j<knmaxtrk;j++) {
604             dig.SetTrack(j,-3);
605             dig.SetHit(j,-1);
606         }
607         aliITS->AddSimDigit(AliITSUGeomTGeo::kChipTypePix, &dig);
608     }
609     //
610 }
611
612 //______________________________________________________________________
613 Int_t AliITSUSimulationPix::AddRandomNoisePixels(Double_t tof)
614 {
615     // create random noisy sdigits above threshold
616     //
617     int modId = fChip->GetIndex();
618     int npix = fSeg->GetNPads();
619     int ncand = gRandom->Poisson( npix*fSimuParam->GetPixFakeRate() );
620     if (ncand<1) return 0;
621     //
622     double probNoisy,noiseSig,noiseMean,thresh;
623     //
624     UInt_t row,col;
625     Int_t iCycle;
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();
630     //
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
635        //
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);
640         }
641         return ncand;
642     }
643     else
644     {
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);
649         }
650         
651         
652     }
653 }
654
655
656 //______________________________________________________________________
657 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old)
658 {
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).
665     UInt_t col,row;
666     Int_t iCycle;
667     Double_t pulse1,pulse2;
668     Double_t couplR=0.0,couplC=0.0;
669     Double_t xr=0.;
670     //
671     fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
672     int cycle = 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) {
679         //
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);
684         //
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);
689     }
690     //
691 }
692
693 //______________________________________________________________________
694 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old)
695 {
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.
700     // Inputs:
701     // old            existing AliITSUSDigit
702     // ntrack         track incex number
703     // idhit          hit index number
704     // chip         chip number
705     //
706     UInt_t col,row;
707     Int_t cycle;
708     Int_t modId = fChip->GetIndex();
709     Double_t pulse1,pulse2;
710     Double_t couplR=0.0,couplC=0.0;
711     //
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));
715     //
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();
719         //
720         int j1 = int(col)+isign;
721         pulse1 *= couplC;
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);
724         
725         // loop in row direction
726         int j2 = int(row) + isign;
727         pulse2 *= couplR;
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);
730     } // for isign
731 }
732
733 //______________________________________________________________________
734 void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
735 {
736     // attach response parameterisation data
737     fResponseParam = resp;
738     //
739     int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
740     const char* hname = 0;
741     fSpread2DHisto = 0;
742     //
743     switch (spreadID) {
744             //
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));
750             break;
751             //
752         case kSpreadFunDoubleGauss2D:
753             fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
754             break;
755             //
756         case kSpreadFunGauss2D:
757             fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
758             break;
759             //
760         default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
761     }
762     //
763     int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
764     switch (readoutType) {
765         case kReadOutStrobe:
766             fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
767             break;
768         case kReadOutRollingShutter:
769             fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
770             break;
771         default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
772     }
773     
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);
787     //
788     if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
789     else fSimuParam->SetPixNoiseInAllMod(kFALSE);
790     //
791     //  Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
792     fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
793     
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  ============================"));
805     
806 }
807
808 //______________________________________________________________________
809 Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
810 {
811     //
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
816     //
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));
821     return cycle;
822     //
823 }
824
825 //______________________________________________________________________
826 Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
827 {
828     //
829     // Check whether the hit is in the read out window of the given column/row of the sensor
830     //
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;
835     //
836 }
837
838 //_______________________________________________________________________
839 void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xrow, Int_t zcol, Float_t &x, Float_t &z)
840 {
841     //
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
846     //
847     ((AliITSUSegmentationPix*)fSeg)->GetDiodShift(xrow,zcol,x,z);
848     //
849 }
850
851 //______________________________________________________________________
852 void AliITSUSimulationPix::Hits2SDigitsFastDigital()
853 {
854     // Does the digital chip response simulation
855     //    AliITSUChip *mod  Pointer to this chip
856     //
857     
858     
859     TObjArray *hits = fChip->GetHits();
860     Int_t nhits = hits->GetEntriesFast();
861     if (nhits<=0) return;
862     //
863     Int_t h,ix,iz;
864     Int_t idtrack;
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;
868     
869     //
870     for (h=0;h<nhits;h++) {
871         //
872         if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
873         //
874         st = Sqrt(x1*x1+y1*y1+z1*z1);
875         
876         //___ place hit to the middle of the path segment - CHANGE LATER !
877         x   = (x0+x1)/2.0;
878         y   = (y0+y1)/2.0;
879         z   = (z0+z1)/2.0;
880         //
881         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
882         //
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;
886         //
887         //
888         el  = de / fSimuParam->GetGeVToCharge();
889         //
890         PlaceDigitalPixels(x,z,el,tof,idtrack,h);
891         //
892     } // Loop over all hits h
893     
894 }
895 //______________________________________________________________________
896 void AliITSUSimulationPix::PlaceDigitalPixels(Double_t x0hit,Double_t z0hit, Double_t el, Double_t tof, Int_t tID, Int_t hID)
897 {
898     // Place the digital pixel positions on the sensor
899     // Inputs:
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
907     //
908     
909     
910     Int_t ix,iz,nx,nz;
911     Float_t x,z;   // keep coordinates float (required by AliSegmentation)
912     Float_t distX = 0, distZ = 0;
913     
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):
916     //  -+-
917     //  +++
918     //  -+-
919     //
920     //___ This list should come from a look up table based on CluTypeID as well as COG coord
921     
922     //
923     TRandom3 rnd;
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
926     //
927     x = x0hit + distX;
928     z = z0hit + distZ;
929     //
930     if(!fSeg->LocalToDet(x,z,ix,iz)) return; // if clu CoG is outside of the chip skipp the cluster -> refine later
931     //
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};
935     //
936     Double_t s = el / 1.0 / nCluPixels;
937     //
938     int cycleRO;
939     //
940     for (Int_t ipix = 0 ; ipix < nCluPixels; ipix++)
941     {
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
946     }
947    
948
949     
950 }
951
952
953