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