0d2b99f8e15e779c05a00eab7e63ee457666510d
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimulationPix.cxx
1 /*
2  questions to experts: why RemoveDeadPixels should be called before FrompListToDigits ? 
3
4  
5 */
6
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 *                                                                        *
10 * Author: The ALICE Off-line Project.                                    *
11 * Contributors are mentioned in the code where appropriate.              *
12 *                                                                        *
13 * Permission to use, copy, modify and distribute this software and its   *
14 * documentation strictly for non-commercial purposes is hereby granted   *
15 * without fee, provided that the above copyright notice appears in all   *
16 * copies and that both the copyright notice and this permission notice   *
17 * appear in the supporting documentation. The authors make no claims     *
18 * about the suitability of this software for any purpose. It is          *
19 * provided "as is" without express or implied warranty.                  *
20 **************************************************************************/
21
22 #include <Riostream.h>
23 #include <TGeoGlobalMagField.h>
24 #include <TH1.h>
25 #include <TString.h>
26 #include "AliITSU.h"
27 #include "AliITSUDigitPix.h"
28 #include "AliITSUHit.h"
29 #include "AliITSUModule.h"
30 #include "AliITSUSensMap.h"
31 #include "AliITSUCalibrationPix.h"
32 #include "AliITSUSegmentationPix.h"
33 #include "AliITSUSimulationPix.h"
34 #include "AliLog.h"
35 #include "AliRun.h"
36 #include "AliMagF.h"
37 #include "AliMathBase.h"
38 #include "AliITSUSimuParam.h"
39 #include "AliITSUSDigit.h"
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 ////////////////////////////////////////////////////////////////////////
58
59 //______________________________________________________________________
60 AliITSUSimulationPix::AliITSUSimulationPix()
61 :  fTanLorAng(0)
62   ,fStrobe(kTRUE)
63   ,fStrobeLenght(4)
64   ,fStrobePhase(-12.5e-9)
65 {
66    // Default constructor.
67 }
68
69 //______________________________________________________________________
70 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
71   :AliITSUSimulation(sim,map)
72   ,fTanLorAng(0)
73   ,fStrobe(kTRUE)
74   ,fStrobeLenght(4)
75   ,fStrobePhase(-12.5e-9)
76 {
77   // standard constructor
78   Init();
79 }
80
81 //______________________________________________________________________
82 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s) 
83   :AliITSUSimulation(s)
84   ,fTanLorAng(s.fTanLorAng)
85   ,fStrobe(s.fStrobe)
86   ,fStrobeLenght(s.fStrobeLenght)
87   ,fStrobePhase(s.fStrobePhase)
88 {
89   //     Copy Constructor
90 }
91
92
93 //______________________________________________________________________
94 AliITSUSimulationPix::~AliITSUSimulationPix()
95 {
96   // destructor
97   // only the sens map is owned and it is deleted by ~AliITSUSimulation
98 }
99
100 //______________________________________________________________________
101 AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
102 {
103   //    Assignment operator
104   if (&s == this) return *this;
105   AliITSUSimulation::operator=(s);
106   fStrobe       = s.fStrobe;
107   fStrobeLenght = s.fStrobeLenght;
108   fStrobePhase  = s.fStrobePhase;
109   return *this;
110 }
111
112 //______________________________________________________________________
113 void AliITSUSimulationPix::Init()
114 {
115   // Initilization
116   if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
117 }
118
119 //______________________________________________________________________
120 Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole) 
121 {
122   // This function set the Tangent of the Lorentz angle. 
123   // A weighted average is used for electrons and holes 
124   // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
125   // output: Bool_t : kTRUE in case of success
126   //
127   if (weightHole<0) {
128     weightHole=0.;
129     AliWarning("You have asked for negative Hole weight");
130     AliWarning("I'm going to use only electrons");
131   }
132   if (weightHole>1) {
133     weightHole=1.;
134     AliWarning("You have asked for weight > 1");
135     AliWarning("I'm going to use only holes");
136   }
137   Double_t weightEle=1.-weightHole;
138   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
139   if (!fld) AliFatal("The field is not initialized");
140   Double_t bz = fld->SolenoidField();
141   fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
142                           weightEle*fSimuParam->LorentzAngleElectron(bz));
143   return kTRUE;
144 }
145
146 //_____________________________________________________________________
147 void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg)
148 {
149   //  This function begins the work of creating S-Digits.
150   if (!(mod->GetNHits())) {
151     AliDebug(1,Form("In event %d module %d there are %d hits returning.",
152                     event, mod->GetIndex(),mod->GetNHits()));
153     return;
154   } 
155   //
156   if (fStrobe) if (event != fEvent) GenerateStrobePhase(); 
157   InitSimulationModule(mod->GetIndex(), event, seg );
158   // Hits2SDigits(mod);
159   Hits2SDigitsFast(mod);
160   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
161   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();  
162   WriteSDigits();
163   ClearMap();
164 }
165
166 //______________________________________________________________________
167 void AliITSUSimulationPix::WriteSDigits()
168 {
169   //  This function adds each S-Digit to pList
170   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
171   int nsd = fSensMap->GetEntries();
172   for (int i=0;i<nsd;i++) {
173     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
174     if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
175     aliITS->AddSumDigit(*sd);
176   }
177   return; 
178 }
179
180 //______________________________________________________________________
181 void AliITSUSimulationPix::FinishSDigitiseModule()
182 {
183    //  This function calls SDigitsToDigits which creates Digits from SDigits
184   FrompListToDigits();
185   ClearMap();
186   return;
187 }
188
189 //______________________________________________________________________
190 void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg)
191 {
192   //  This function creates Digits straight from the hits and then adds
193   //  electronic noise to the digits before adding them to pList
194   //  Each of the input variables is passed along to Hits2SDigits
195   //
196   if (fStrobe) if (event != fEvent) GenerateStrobePhase();
197   InitSimulationModule( mod->GetIndex(), event, seg );
198   // Hits2SDigits(mod);
199   Hits2SDigitsFast(mod);
200   //
201   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
202   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
203   FrompListToDigits();
204   ClearMap();
205 }
206
207 //______________________________________________________________________
208 void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
209 {
210   // Does the charge distributions using Gaussian diffusion charge charing.
211   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
212   const Double_t kcmtomic = 1.e4;
213   const Double_t kBunchLenght = 25e-9; // LHC clock
214   Int_t nhits = mod->GetNHits();
215   if (!nhits) return;
216   //
217   Int_t h,ix,iz,i;
218   Int_t idtrack;
219   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
220   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;
221   Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
222   Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half Thickness
223   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
224   //
225   for (h=0;h<nhits;h++) {
226     if (fStrobe && 
227         ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
228          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
229         ) continue;
230     //
231     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
232     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
233     if (st>0.0) {
234       st = (Double_t)((Int_t)(st*kcmtomic)); // number of microns
235       if (st<=1.0) st = 1.0;
236       dt = 1.0/st;
237       for (t=0.0;t<1.0;t+=dt) { // Integrate over t
238         tp  = t+0.5*dt;
239         x   = x0+x1*tp;
240         y   = y0+y1*tp;
241         z   = z0+z1*tp;
242         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
243         el  = dt * de / fSimuParam->GetGeVToCharge();
244         //
245         if (GetDebug(1)) if (el<=0.0) cout<<"el="<<el<<" dt="<<dt<<" de="<<de<<endl; // end if GetDebug
246         //
247         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); 
248         sigx=sig;
249         sigz=sig*fda;
250         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
251         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
252       } // end for t
253     } else { // st == 0.0 deposit it at this point
254       x   = x0;
255       y   = y0;
256       z   = z0;
257       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
258       el  = de / fSimuParam->GetGeVToCharge();
259       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
260       sigx = sig;
261       sigz = sig*fda;
262       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
263       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
264     } // end if st>0.0    
265   } // Loop over all hits h
266   //
267   // Coupling
268   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
269   AliITSUSDigit* dg = 0;
270   switch (fSimuParam->GetPixCouplingOption()) {
271   case AliITSUSimuParam::kNewCouplingPix :
272     for (i=nd;i--;) {
273       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
274       if (fSensMap->IsDisabled(dg)) continue;
275       SetCoupling(dg,idtrack,h);
276     } 
277     break;
278   case AliITSUSimuParam::kOldCouplingPix:
279     for (i=nd;i--;) {
280       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
281       if (fSensMap->IsDisabled(dg)) continue;
282       SetCouplingOld(dg,idtrack,h);
283     } 
284     break;
285   default:
286     break;
287     
288   } // end switch
289   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
290 }
291
292 //______________________________________________________________________
293 void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
294 {
295   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
296   //    AliITSUModule *mod  Pointer to this module
297   // Output:
298   //    none.
299   // Return:
300   //    none.
301   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
302   const Int_t kn10=10;
303   const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
304                             4.325316833e-1,4.869532643e-1,5.130467358e-1,
305                             5.674683167e-1,6.602952159e-1,7.833023029e-1,
306                             9.255628306e-1};
307   const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
308                             7.472567455e-2,3.333567215e-2,3.333567215e-2,
309                             7.472567455e-2,1.095431813e-1,1.346333597e-1,
310                             1.477621124e-1};
311   const Double_t kBunchLenght = 25e-9; // LHC clock
312   TObjArray *hits = mod->GetHits();
313   Int_t nhits = hits->GetEntriesFast();
314   if (nhits<=0) return;
315   //
316   Int_t h,ix,iz,i;
317   Int_t idtrack;
318   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
319   Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; 
320   Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
321   Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half thickness
322   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
323   //
324   for (h=0;h<nhits;h++) {
325     // Check if the hit is inside readout window
326     if (fStrobe && 
327         ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
328          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
329         ) continue;
330     //
331     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
332     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
333     if (st>0.0) 
334       for (i=0;i<kn10;i++) { // Integrate over t
335         t   = kti[i];
336         x   = x0+x1*t;
337         y   = y0+y1*t;
338         z   = z0+z1*t;
339         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
340         el  = kwi[i]*de/fSimuParam->GetGeVToCharge();
341         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
342         sigx=sig;
343         sigz=sig*fda;
344         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
345         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
346         //                cout << "sigx sigz " << sigx << " " << sigz << endl; // dom
347       } // end for i // End Integrate over t
348     else { // st == 0.0 deposit it at this point
349       x   = x0;
350       y   = y0;
351       z   = z0;
352       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
353       el  = de / fSimuParam->GetGeVToCharge();
354       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
355       sigx=sig;
356       sigz=sig*fda;
357       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
358       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
359     } // end if st>0.0
360     
361   } // Loop over all hits h
362   
363   // Coupling
364   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
365   AliITSUSDigit* dg = 0;
366   switch (fSimuParam->GetPixCouplingOption()) {
367   case AliITSUSimuParam::kNewCouplingPix :
368     for (i=nd;i--;) {
369       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
370       if (fSensMap->IsDisabled(dg)) continue;
371       SetCoupling(dg,idtrack,h);
372     } 
373   case AliITSUSimuParam::kOldCouplingPix:
374     for (i=nd;i--;) {
375       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
376       if (fSensMap->IsDisabled(dg)) continue;
377       SetCouplingOld(dg,idtrack,h);
378     } 
379     break;
380   default:
381     break;
382   } // end switch
383   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
384 }
385
386 //______________________________________________________________________
387 void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
388                                           Int_t ix0,Int_t iz0,
389                                           Double_t el,Double_t sig,Double_t ld,
390                                           Int_t t,Int_t hi)
391 {
392    // Spreads the charge over neighboring cells. Assume charge is distributed
393    // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
394    // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
395    // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
396    // Defined this way, the integral over all x and z is el.
397    // Inputs:
398    //    Double_t x0   x position of point where charge is liberated
399    //    Double_t z0   z position of point where charge is liberated
400    //    Int_t    ix0  row of cell corresponding to point x0
401    //    Int_t    iz0  columb of cell corresponding to point z0
402    //    Double_t el   number of electrons liberated in this step
403    //    Double_t sig  Sigma difusion for this step (y0 dependent)
404    //    Double_t ld   lorentz drift in x for this step (y0 dependent)
405    //    Int_t    t    track number
406    //    Int_t    ti   hit track index number
407    //    Int_t    hi   hit "hit" index number
408    // Outputs:
409    //     none.
410    // Return:
411    //     none.
412    const Int_t knx = 3,knz = 2;
413    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
414    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
415    Int_t ix,iz,ixs,ixe,izs,ize;
416    Float_t x,z;  // keep coordinates float (required by AliSegmentation)
417    Double_t s,sp,x1,x2,z1,z2; 
418    //
419    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));
420    if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
421      UpdateMapSignal(iz0,ix0,t,hi,el);
422      return;
423    } // end if
424    sp = 1.0/(sig*kRoot2);
425    ixs = TMath::Max(-knx+ix0,0);
426    ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
427    izs = TMath::Max(-knz+iz0,0);
428    ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
429    for (ix=ixs;ix<=ixe;ix++) 
430      for (iz=izs;iz<=ize;iz++) {
431        fSeg->DetToLocal(ix,iz,x,z); // pixel center
432        double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
433        double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
434        x1  = x;
435        z1  = z;
436        x2  = x1 + dxi; // Upper
437        x1 -= dxi;  // Lower
438        z2  = z1 + dzi; // Upper
439        z1 -= dzi;  // Lower
440        x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
441        x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
442        z1 -= z0; // Distance from where track traveled
443        z2 -= z0; // Distance from where track traveled
444        s   = el*0.25; // Correction based on definision of Erfc
445        s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
446        if (GetDebug(3)) {
447          cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
448            " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<< 
449            " sp*x1="<<sp*x1<<" sp*x2="<<sp*x2<<" s="<<s;
450        } // end if GetDebug
451        s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
452        if (GetDebug(3)) cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl; // end if GetDebug
453        if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
454      } // end for ix, iz
455    //
456 }
457
458 //______________________________________________________________________
459 void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
460                                             Int_t ix0,Int_t iz0,
461                                             Double_t el,Double_t sigx,Double_t sigz,
462                                             Double_t ld,Int_t t,Int_t hi)
463 {
464   // Spreads the charge over neighboring cells. Assume charge is distributed
465   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
466   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
467   // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
468   // Defined this way, the integral over all x and z is el.
469   // Inputs:
470   //    Double_t x0   x position of point where charge is liberated
471   //    Double_t z0   z position of point where charge is liberated
472   //    Int_t    ix0  row of cell corresponding to point x0
473   //    Int_t    iz0  columb of cell corresponding to point z0
474   //    Double_t el   number of electrons liberated in this step
475   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
476   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
477   //    Double_t ld   lorentz drift in x for this stip (y0 dependent)
478   //    Int_t    t    track number
479   //    Int_t    ti   hit track index number
480   //    Int_t    hi   hit "hit" index number
481   // Outputs:
482   //     none.
483   // Return:
484   //     none.
485   const Int_t knx = 3,knz = 3; // RS: TO TUNE
486   const Double_t kRoot2 = 1.414213562; // Sqrt(2).
487   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
488   Int_t ix,iz,ixs,ixe,izs,ize;
489   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
490   Double_t s,spx,spz,x1,x2,z1,z2; 
491   //
492   if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d,ld=%e)",
493                                 x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
494   if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
495     UpdateMapSignal(iz0,ix0,t,hi,el);
496     if (GetDebug(2)) {
497       cout << "sigx<=0.0=" << sigx << endl;
498       cout << "sigz<=0.0=" << sigz << endl;
499     } // end if GetDebug
500     return;
501   } // end if
502   spx = 1.0/(sigx*kRoot2);     
503   spz = 1.0/(sigz*kRoot2);
504   if (GetDebug(2)) {
505     cout << "sigx=" << sigx << " spx=" << spx << endl;
506     cout << "sigz=" << sigz << " spz=" << spz << endl;
507   } // end if GetDebug
508   ixs = TMath::Max(-knx+ix0,0);
509   ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
510   izs = TMath::Max(-knz+iz0,0);
511   ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
512   for (ix=ixs;ix<=ixe;ix++) 
513     for (iz=izs;iz<=ize;iz++) {
514       fSeg->DetToLocal(ix,iz,x,z); // pixel center
515       double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
516       double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
517       x1  = x;
518       z1  = z;
519       x2  = x1 + dxi; // Upper
520       x1 -= dxi;      // Lower
521       z2  = z1 + dzi; // Upper
522       z1 -= dzi;      // Lower
523       x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
524       x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
525       z1 -= z0; // Distance from where track traveled
526       z2 -= z0; // Distance from where track traveled
527       s   = el*0.25; // Correction based on definision of Erfc
528       s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
529       if (GetDebug(3)) cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<" iz0="<<iz0<<" iz="<<iz<<" z0="<<z<<" spx*x1="<<spx*x1<<" spx*x2="<<spx*x2<<" s="<<s; // end if GetDebug
530       s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
531       if (GetDebug(3)) cout<<" spz*z1="<<spz*z1<<" spz*z2="<<spz*z2<<" s="<<s<< endl; // end if GetDebug
532       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
533     } // end for ix, iz
534   //
535 }
536
537 //______________________________________________________________________
538 void AliITSUSimulationPix::RemoveDeadPixels() 
539 {
540   // Removes dead pixels on each module (ladder)
541   // This should be called before going from sdigits to digits (FrompListToDigits)
542   
543   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
544   if (!calObj) return;
545   //
546   if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
547   //
548   // remove single bad pixels one by one
549   int nsingle = calObj->GetNrBadSingle();
550   UInt_t col,row;
551   for (int i=nsingle;i--;) {
552     calObj->GetBadPixelSingle(i,row,col);
553     fSensMap->DeleteItem(col,row);
554   }
555   int nsd = fSensMap->GetEntriesUnsorted();
556   for (int isd=nsd;isd--;) {
557     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
558     if (fSensMap->IsDisabled(sd)) continue;
559     fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
560     int chip = fSeg->GetChipFromChannel(0,col);
561     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
562     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
563   }
564   //
565 }
566
567 //______________________________________________________________________
568 void AliITSUSimulationPix::AddNoisyPixels() 
569 {
570   // Adds noisy pixels on each module (ladder)
571   // This should be called before going from sdigits to digits (FrompListToDigits)
572   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
573   if (!calObj) return;
574   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), 
575                                                         10*fSimuParam->GetPixThreshold(fModule));
576   //
577 }
578
579 //______________________________________________________________________
580 void AliITSUSimulationPix::FrompListToDigits() 
581 {
582   // add noise and electronics, perform the zero suppression and add the
583   // digit to the list
584   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
585   UInt_t ix,iz;
586   Double_t sig;
587   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
588   static AliITSUDigitPix dig;
589   // RS: in principle:
590   // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold. 
591   // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
592   // With many channels this will be too time consuming, hence I do the following
593   // 1) Precalculate the probability that the nois alone will exceed the threshold. 
594   // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold.
595   // 3) For pixels having a hits apply the usual noise and compare with threshold
596   //
597   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
598   //
599   int maxInd = fSensMap->GetMaxIndex();
600   double minProb = 0.1/maxInd;
601   //
602   int nsd = fSensMap->GetEntries();
603   Int_t prevID=0,curID=0;
604   TArrayI ordSampleInd(100),ordSample(100);
605   //
606   double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
607   fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
608   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
609   //
610   for (int i=0;i<nsd;i++) {
611     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
612     if (fSensMap->IsDisabled(sd)) continue;
613     curID = (int)sd->GetUniqueID();
614     //
615     if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current
616       CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
617       prevID = curID+1;
618     }
619     //
620     if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
621     if (TMath::Abs(sig)>2147483647.0) { //RS?
622       //PH 2147483647 is the max. integer
623       //PH This apparently is a problem which needs investigation
624       AliWarning(Form("Too big or too small signal value %f",sig));
625     }
626     fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix);
627     dig.SetCoord1(iz);
628     dig.SetCoord2(ix);
629     dig.SetSignal(1);
630     dig.SetSignalPix((Int_t)sig);
631     int ntr = sd->GetNTracks();
632     for (int j=0;j<ntr;j++) {
633       dig.SetTrack(j,sd->GetTrack(j));
634       dig.SetHit(j,sd->GetHit(j));
635     }
636     for (int j=ntr;j<knmaxtrk;j++) {
637       dig.SetTrack(j,-3);
638       dig.SetHit(j,-1);
639     }
640     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
641   }
642   // if needed, add noisy pixels with id from last real hit to maxID
643   if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
644   // 
645 }
646
647 //______________________________________________________________________
648 Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
649 {
650   // create random noisy digits above threshold within id range [minID,maxID[
651   // see FrompListToDigits for details
652   //
653   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
654   UInt_t ix,iz;
655   static AliITSUDigitPix dig;
656   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
657   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
658   //
659   Int_t ncand = 0;
660   int npix = maxID-minID;
661   if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added
662   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd); 
663   int* ordV = ordSample.GetArray();
664   int* ordI = ordSampleInd.GetArray();
665   for (int j=0;j<ncand;j++) {
666     fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix);   // create noisy digit
667     dig.SetCoord1(iz);
668     dig.SetCoord2(ix);
669     dig.SetSignal(1);
670     dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
671     for (int k=knmaxtrk;k--;) {
672       dig.SetTrack(k,-3);
673       dig.SetHit(k,-1);
674     }
675     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
676     if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
677   }
678   return ncand;
679 }
680
681 //______________________________________________________________________
682 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit) 
683 {
684   //  Take into account the coupling between adiacent pixels.
685   //  The parameters probcol and probrow are the probability of the
686   //  signal in one pixel shared in the two adjacent pixels along
687   //  the column and row direction, respectively.
688   //  Note pList is goten via GetMap() and module is not need any more.
689   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
690   //Begin_Html
691   /*
692     <img src="picts/ITS/barimodel_3.gif">
693      </pre>
694      <br clear=left>
695      <font size=+2 color=red>
696      <a href="mailto:tiziano.virgili@cern.ch"></a>.
697      </font>
698      <pre>
699    */
700    //End_Html
701    // Inputs:
702   // old                  existing AliITSUSDigit
703   // Int_t ntrack         track incex number
704   // Int_t idhit          hit index number
705   UInt_t col,row;
706   Double_t pulse1,pulse2;
707   Double_t couplR=0.0,couplC=0.0;
708   Double_t xr=0.;
709   //
710   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
711   fSimuParam->GetPixCouplingParam(couplC,couplR);
712   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
713                                 col,row,ntrack,idhit,couplC,couplR));
714   pulse2 = pulse1 = old->GetSignal();
715   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
716   for (Int_t isign=-1;isign<=1;isign+=2) {
717     //
718     // loop in col direction
719     int j1 = int(col) + isign;
720     xr = gRandom->Rndm();
721     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
722     //
723     // loop in row direction
724     int j2 = int(row) + isign;
725     xr = gRandom->Rndm();
726     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
727   } 
728   //
729 }
730
731 //______________________________________________________________________
732 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit) 
733 {
734   //  Take into account the coupling between adiacent pixels.
735   //  The parameters probcol and probrow are the fractions of the
736   //  signal in one pixel shared in the two adjacent pixels along
737   //  the column and row direction, respectively.
738   //Begin_Html
739   /*
740     <img src="picts/ITS/barimodel_3.gif">
741     </pre>
742     <br clear=left>
743     <font size=+2 color=red>
744     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
745     </font>
746     <pre>
747   */
748   //End_Html
749   // Inputs:
750   // old            existing AliITSUSDigit
751   // ntrack         track incex number
752   // idhit          hit index number
753   // module         module number
754   //
755   UInt_t col,row;
756   Double_t pulse1,pulse2;
757   Double_t couplR=0.0,couplC=0.0;
758   //
759   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
760   fSimuParam->GetPixCouplingParam(couplC,couplR);
761   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
762                                 col,row,ntrack,idhit,couplC,couplR));
763  //
764  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
765  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
766    pulse2 = pulse1 = old->GetSignal();
767    //
768    int j1 = int(col)+isign;
769    pulse1 *= couplC;    
770    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
771    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
772    
773    // loop in row direction
774    int j2 = int(row) + isign;
775    pulse2 *= couplR;
776    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
777    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
778  } // for isign
779 }
780
781 //______________________________________________________________________
782 void AliITSUSimulationPix::GenerateStrobePhase()
783 {
784   // Generate randomly the strobe
785   // phase w.r.t to the LHC clock
786   // Done once per event
787   const Double_t kBunchLenght = 25e-9; // LHC clock
788   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
789     (Double_t)fStrobeLenght*kBunchLenght+
790     kBunchLenght/2;
791 }
792