Decoupled ITS/UPGRADE cmake stuff from ITS
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSsimulationPixUpg.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 #include <Riostream.h>
17 #include <TGeoGlobalMagField.h>
18 #include <TH1.h>
19 #include <TString.h>
20 #include "AliITS.h"
21 #include "AliITSdigitPixUpg.h"
22 #include "AliITShit.h"
23 #include "AliITSmodule.h"
24 #include "AliITSpList.h"
25 #include "AliITSCalibrationPixUpg.h"
26 #include "AliITSsegmentationPixUpg.h"
27 #include "AliITSsimulationPixUpg.h"
28 #include "AliLog.h"
29 #include "AliRun.h"
30 #include "AliMagF.h"
31 #include "AliMathBase.h"
32
33 //#define DEBUG
34
35 using std::endl;
36 using std::cout;
37 ClassImp(AliITSsimulationPixUpg)
38 ////////////////////////////////////////////////////////////////////////
39 //  Version: 1
40 //  Modified by D. Elia, G.E. Bruno, H. Tydesjo 
41 //  Fast diffusion code by Bjorn S. Nilsen
42 //  March-April 2006
43 //  October     2007: GetCalibrationObjects() removed
44 //
45 //  Version: 0
46 //  Written by Boris Batyunya
47 //  December 20 1999
48 //
49 //  Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
50 //
51 // AliITSsimulationPixUpg is to do the simulation of pixels
52 //
53 ////////////////////////////////////////////////////////////////////////
54
55 //______________________________________________________________________
56 AliITSsimulationPixUpg::AliITSsimulationPixUpg():
57 AliITSsimulation(),
58 fHistos(0),
59 fHistName(),
60 fCoupling(),
61 fLorentz(kFALSE),
62 fTanLorAng(0),
63 fStrobe(kTRUE),
64 fStrobeLenght(4),
65 fStrobePhase(-12.5e-9)
66 {
67    // Default constructor.
68    // Inputs:
69    //    none.
70    // Outputs:
71    //    none.
72    // Return:
73    //    A default constructed AliITSsimulationPixUpg class.
74
75    AliDebug(1,Form("Calling default constructor"));
76 //    Init();
77 }
78
79 //______________________________________________________________________
80 AliITSsimulationPixUpg::AliITSsimulationPixUpg(AliITSDetTypeSim *dettyp):
81 AliITSsimulation(dettyp),
82 fHistos(0),
83 fHistName(),
84 fCoupling(),
85 fLorentz(kFALSE),
86 fTanLorAng(0),
87 fStrobe(kTRUE),
88 fStrobeLenght(4),
89 fStrobePhase(-12.5e-9)
90 {
91    // standard constructor
92    // Inputs:
93    //    AliITSsegmentation *seg  A pointer to the segmentation class
94    //                             to be used for this simulation
95    //    AliITSCalibration     *resp A pointer to the responce class to
96    //                             be used for this simulation
97    // Outputs:
98    //    none.
99    // Return:
100    //    A default constructed AliITSsimulationPixUpg class.
101
102    AliDebug(1,Form("Calling standard constructor "));
103    Init();
104 }
105
106 //______________________________________________________________________
107 void AliITSsimulationPixUpg::Init()
108 {
109   // Initilization
110   // Inputs:
111   //    none.
112   // Outputs:
113   //    none.
114   // Return:
115   //    none.
116   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
117   
118   SetModuleNumber(0);
119   SetEventNumber(0);
120   SetMap(new AliITSpList(GetNPixelsZ(),GetNPixelsX()));
121   AliITSSimuParam* simpar = fDetType->GetSimuParam();
122   AliITSsegmentationPixUpg* seg = (AliITSsegmentationPixUpg*)GetSegmentationModel(-1);
123   Double_t bias = simpar->GetPixBiasVoltage();
124 //    cout << "Bias Voltage --> " << bias << endl; // dom    
125    simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),bias);
126 // set kind of coupling ("old" or "new")
127    char opt[20];
128    simpar->GetPixCouplingOption(opt);
129    char *old = strstr(opt,"old");
130    if (old) {
131      fCoupling=2;
132    } else {
133      fCoupling=1;
134    } // end if
135    SetLorentzDrift(simpar->GetSPDLorentzDrift());
136    if (fLorentz) SetTanLorAngle(simpar->GetSPDLorentzHoleWeight());
137    //SetStrobeGeneration(kFALSE);
138    if (fStrobe) GenerateStrobePhase();
139 }
140
141 //______________________________________________________________________
142 Bool_t AliITSsimulationPixUpg::SetTanLorAngle(Double_t weightHole) 
143 {
144     // This function set the Tangent of the Lorentz angle. 
145     // A weighted average is used for electrons and holes 
146     // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
147     // output: Bool_t : kTRUE in case of success
148     //
149     if (!fDetType) {
150       AliError("AliITSsimulationPixUpg::SetTanLorAngle: AliITSDetTypeSim* fDetType not set ");
151       return kFALSE;}
152     if (weightHole<0) {
153        weightHole=0.;
154        AliWarning("AliITSsimulationPixUpg::SetTanLorAngle: You have asked for negative Hole weight");
155        AliWarning("AliITSsimulationPixUpg::SetTanLorAngle: I'm going to use only electrons");
156     }
157     if (weightHole>1) {
158        weightHole=1.;
159        AliWarning("AliITSsimulationPixUpg::SetTanLorAngle: You have asked for weight > 1");
160        AliWarning("AliITSsimulationPixUpg::SetTanLorAngle: I'm going to use only holes");
161     }
162     Double_t weightEle=1.-weightHole;
163     AliITSSimuParam* simpar = fDetType->GetSimuParam();
164     AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
165     if (!fld) AliFatal("The field is not initialized");
166     Double_t bz = fld->SolenoidField();
167     fTanLorAng = TMath::Tan(weightHole*simpar->LorentzAngleHole(bz) +
168                              weightEle*simpar->LorentzAngleElectron(bz));
169     return kTRUE;
170 }
171
172 //______________________________________________________________________
173 AliITSsimulationPixUpg::~AliITSsimulationPixUpg()
174 {
175    // destructor
176    // Inputs:
177    //    none.
178    // Outputs:
179    //    none.
180    // Return:
181    //     none.
182
183    if (fHistos) {
184        fHistos->Delete(); 
185        delete fHistos;     
186    } // end if fHistos
187 }
188
189 //______________________________________________________________________
190 AliITSsimulationPixUpg::AliITSsimulationPixUpg(const AliITSsimulationPixUpg &s) : 
191   AliITSsimulation(s),
192   fHistos(s.fHistos),
193   fHistName(s.fHistName),
194   fCoupling(s.fCoupling),
195   fLorentz(s.fLorentz),
196   fTanLorAng(s.fTanLorAng),
197   fStrobe(s.fStrobe),
198   fStrobeLenght(s.fStrobeLenght),
199   fStrobePhase(s.fStrobePhase)
200 {
201   //     Copy Constructor
202    // Inputs:
203    //    AliITSsimulationPixUpg &s The original class for which
204    //                                this class is a copy of
205    // Outputs:
206    //    none.
207    // Return:
208
209 }
210
211 //______________________________________________________________________
212 AliITSsimulationPixUpg&  AliITSsimulationPixUpg::operator=(const AliITSsimulationPixUpg &s)
213 {
214    //    Assignment operator
215    // Inputs:
216    //    AliITSsimulationPixUpg &s The original class for which
217    //                                this class is a copy of
218    // Outputs:
219    //    none.
220    // Return:
221
222    if (&s == this) return *this;
223    this->fHistos = s.fHistos;
224    fCoupling  = s.fCoupling;
225    fHistName   = s.fHistName;
226    fLorentz   = s.fLorentz;
227    fTanLorAng = s.fTanLorAng;
228    fStrobe       = s.fStrobe;
229    fStrobeLenght = s.fStrobeLenght;
230    fStrobePhase  = s.fStrobePhase;
231    return *this;
232 }
233
234 //______________________________________________________________________
235 void AliITSsimulationPixUpg::InitSimulationModule(Int_t module, Int_t event)
236 {
237    //  This function creates maps to build the list of tracks for each
238    //  summable digit. Inputs defined by base class.
239    //  Inputs:
240    //    Int_t module   // Module number to be simulated
241    //    Int_t event    // Event number to be simulated
242    //  Outputs:
243    //    none
244    //  Returns:
245    //    none
246
247    AliDebug(1,Form("(module=%d,event=%d)",module,event));
248    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
249    AliITSSimuParam* simpar = fDetType->GetSimuParam();
250    AliITSsegmentationPixUpg* seg = (AliITSsegmentationPixUpg*)GetSegmentationModel(-1);
251    SetModuleNumber(module);
252    SetEventNumber(event);
253    simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),simpar->GetSPDBiasVoltage(module)); 
254    ClearMap();
255 }
256
257 //_____________________________________________________________________
258 void AliITSsimulationPixUpg::SDigitiseModule(AliITSmodule *mod,Int_t, Int_t event)
259 {
260   //  This function begins the work of creating S-Digits.  Inputs defined
261   //  by base class.
262   //  Inputs:
263   //    AliITSmodule *mod  //  module
264   //    Int_t              //  not used
265   //    Int_t event        //  Event number
266   //  Outputs:
267   //    none
268   //  Return:
269   //    test              //  test returns kTRUE if the module contained hits
270   //                      //  test returns kFALSE if it did not contain hits
271   
272   AliDebug(1,Form("(mod=%p, ,event=%d)",mod,event));
273   if (!(mod->GetNhits())) {
274     AliDebug(1,Form("In event %d module %d there are %d hits returning.",
275                     event, mod->GetIndex(),mod->GetNhits()));
276     return;// if module has no hits don't create Sdigits
277   } // end if
278   SetModuleNumber(mod->GetIndex());
279   if (fStrobe) if (event != GetEventNumber()) GenerateStrobePhase(); 
280   SetEventNumber(event);
281   InitSimulationModule( GetModuleNumber() , event );
282   // HitToSDigit(mod);
283   HitToSDigitFast(mod);
284   if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag())   AddNoisyPixels();
285   if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels();
286   
287   //    cout << "After Remove in SDigitiseModule !!!!!" << endl; // dom
288   //    cout << "Module " << mod->GetIndex() << " Event " << event << endl; // dom
289   WriteSDigits();
290   ClearMap();
291 }
292
293 //______________________________________________________________________
294 void AliITSsimulationPixUpg::WriteSDigits()
295 {
296    //  This function adds each S-Digit to pList
297    //  Inputs:
298    //    none.
299    //  Outputs:
300    //    none.
301    //  Return:
302    //    none
303    Int_t ix, nix, iz, niz;
304    static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
305
306    AliDebug(1,Form("Writing SDigits for module %d",GetModuleNumber()));
307 //    cout << "WriteSDigits for module " << GetModuleNumber() << endl; // dom
308    GetMap()->GetMaxMapIndex(niz, nix);
309    for (iz=0; iz<niz; iz++)for (ix=0; ix<nix; ix++) {
310        if (GetMap()->GetSignalOnly(iz,ix)>0.0) {
311          //            cout << " Signal gt 0  iz ix " << iz << ix << " Module " << GetModuleNumber() << endl; // dom
312          aliITS->AddSumDigit(*(GetMap()->GetpListItem(iz,ix)));
313          if (AliDebugLevel()>0) {
314            AliDebug(1,Form("%d, %d",iz,ix));
315            cout << *(GetMap()->GetpListItem(iz,ix)) << endl;
316          } // end if GetDebug
317        } // end if GetMap()->GetSignalOnly(iz,ix)>0.0
318      } // end for iz,ix
319    return; 
320 }
321
322 //______________________________________________________________________
323 void AliITSsimulationPixUpg::FinishSDigitiseModule()
324 {
325    //  This function calls SDigitsToDigits which creates Digits from SDigits
326    //  Inputs:
327    //    none
328    //  Outputs:
329    //    none
330    //  Return
331    //    none
332
333    AliDebug(1,"()");
334 //    cout << "FinishSDigitiseModule for module " << GetModuleNumber() << endl; // dom
335    FrompListToDigits(); // Charge To Signal both adds noise and
336    ClearMap();
337    return;
338 }
339
340 //______________________________________________________________________
341 void AliITSsimulationPixUpg::DigitiseModule(AliITSmodule *mod,Int_t, Int_t event)
342 {
343    //  This function creates Digits straight from the hits and then adds
344    //  electronic noise to the digits before adding them to pList
345    //  Each of the input variables is passed along to HitToSDigit
346    //  Inputs:
347    //    AliITSmodule *mod     module
348    //    Int_t                 Dummy.
349    //    Int_t                 Dummy
350    //  Outputs:
351    //     none.
352    //  Return:
353    //    none.
354
355   if (fStrobe) if (event != GetEventNumber()) GenerateStrobePhase();
356   AliDebug(1,Form("(mod=%p,,0)",mod));
357   // HitToSDigit(mod);
358   InitSimulationModule( mod->GetIndex(), event );
359   HitToSDigitFast(mod);
360   
361   if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag())   AddNoisyPixels();
362   if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels();
363   //    cout << "After Remove in DigitiseModule in module " << mod->GetIndex() << endl; // dom
364   FrompListToDigits();
365   ClearMap();
366 }
367
368 //______________________________________________________________________
369 void AliITSsimulationPixUpg::HitToSDigit(AliITSmodule *mod)
370 {
371    // Does the charge distributions using Gaussian diffusion charge charing.
372    // Inputs:
373    //    AliITSmodule *mod  Pointer to this module
374    // Output:
375    //    none.
376    // Return:
377    //    none.
378    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
379    const Double_t kBunchLenght = 25e-9; // LHC clock
380    TObjArray *hits = mod->GetHits();
381    Int_t nhits = hits->GetEntriesFast();
382    Int_t h,ix,iz,i;
383    Int_t idtrack;
384    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;
385    Double_t x,y,z,t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
386    AliITSsegmentationPixUpg* seg = (AliITSsegmentationPixUpg*)GetSegmentationModel(-1);
387    AliITSSimuParam *simpar = fDetType->GetSimuParam();
388    Double_t thick = 0.5*kmictocm*seg->Dy();  // Half Thickness
389    simpar->GetSPDSigmaDiffusionAsymmetry(fda);  //
390
391    AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling));
392    if (nhits<=0) return;
393    for (h=0;h<nhits;h++) {
394      if (AliDebugLevel()>0) {
395        AliDebug(1,Form("Hits, %d", h));
396        cout << *(mod->GetHit(h)) << endl;
397      } // end if GetDebug
398      // Check if the hit is inside readout window
399      if (fStrobe)
400        if ((mod->GetHit(h)->GetTOF() < fStrobePhase) ||
401            (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue;
402      if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
403      st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
404      if (st>0.0) {
405        st = (Double_t)((Int_t)(st/kmictocm)); // number of microns
406        if (st<=1.0) st = 1.0;
407        dt = 1.0/st;
408        for (t=0.0;t<1.0;t+=dt) { // Integrate over t
409          tp  = t+0.5*dt;
410          x   = x0+x1*tp;
411          y   = y0+y1*tp;
412          z   = z0+z1*tp;
413          if (!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
414          //el  = res->GeVToCharge((Double_t)(dt*de));
415          el  = dt * de / simpar->GetGeVToCharge();
416          //
417          if (GetDebug(1)) if (el<=0.0) cout<<"el="<<el<<" dt="<<dt<<" de="<<de<<endl; // end if GetDebug
418          //
419          sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y)); 
420          sigx=sig;
421          sigz=sig*fda;
422          if (fLorentz) ld=(y+thick)*fTanLorAng;
423          SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
424        } // end for t
425      } else { // st == 0.0 deposit it at this point
426        x   = x0;
427        y   = y0;
428        z   = z0;
429        if (!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
430        //el  = res->GeVToCharge((Double_t)de);
431        el  = de / simpar->GetGeVToCharge();
432        sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
433        sigx=sig;
434        sigz=sig*fda;
435        if (fLorentz) ld=(y+thick)*fTanLorAng;
436        SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
437      } // end if st>0.0
438      
439    } // Loop over all hits h
440    
441    // Coupling
442    switch (fCoupling) {
443    default:
444      break;
445    case 1: //case 3:
446      for (i=0;i<GetMap()->GetEntries();i++) 
447        if (GetMap()->GetpListItem(i)==0) continue;
448        else{
449          GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
450          SetCoupling(iz,ix,idtrack,h);
451        } // end for i
452      break;
453    case 2: // case 4:
454      for (i=0;i<GetMap()->GetEntries();i++) 
455        if (GetMap()->GetpListItem(i)==0) continue;
456        else{
457          GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
458          SetCouplingOld(iz,ix,idtrack,h);
459        } // end for i
460      break;
461    } // end switch
462    if (GetDebug(2)) Info("HitToSDigit","Finished fCoupling=%d",fCoupling);
463 }
464
465 //______________________________________________________________________
466 void AliITSsimulationPixUpg::HitToSDigitFast(AliITSmodule *mod)
467 {
468   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
469   //    AliITSmodule *mod  Pointer to this module
470   // Output:
471   //    none.
472   // Return:
473   //    none.
474   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
475   const Int_t kn10=10;
476   const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
477                             4.325316833e-1,4.869532643e-1,5.130467358e-1,
478                             5.674683167e-1,6.602952159e-1,7.833023029e-1,
479                             9.255628306e-1};
480   const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
481                             7.472567455e-2,3.333567215e-2,3.333567215e-2,
482                             7.472567455e-2,1.095431813e-1,1.346333597e-1,
483                             1.477621124e-1};
484   const Double_t kBunchLenght = 25e-9; // LHC clock
485   TObjArray *hits = mod->GetHits();
486   Int_t nhits = hits->GetEntriesFast();
487   Int_t h,ix,iz,i;
488   Int_t idtrack;
489   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;
490   Double_t x,y,z,t,st,el,sig,sigx,sigz,fda;
491   AliITSsegmentationPixUpg* seg = (AliITSsegmentationPixUpg*)GetSegmentationModel(-1);
492   AliITSSimuParam* simpar = fDetType->GetSimuParam();
493   Double_t thick = 0.5*kmictocm*seg->Dy();  // Half thickness
494   simpar->GetSPDSigmaDiffusionAsymmetry(fda);
495   //    cout << "Half Thickness " << thick << endl;  // dom
496   //    cout << "Diffusion asymm " << fda << endl;  // dom
497   
498   AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling));
499   if (nhits<=0) return;
500   for (h=0;h<nhits;h++) {
501     if (AliDebugLevel()>0) {
502       AliDebug(1,Form("Hits, %d", h));
503       cout << *(mod->GetHit(h)) << endl;
504     } // end if GetDebug
505     // Check if the hit is inside readout window
506     if (fStrobe)
507       if ((mod->GetHit(h)->GetTOF() < fStrobePhase) ||
508           (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue;
509     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
510     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
511     if (st>0.0) for (i=0;i<kn10;i++) { // Integrate over t
512         t   = kti[i];
513         x   = x0+x1*t;
514         y   = y0+y1*t;
515         z   = z0+z1*t;
516         if (!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
517         el  = kwi[i]*de/simpar->GetGeVToCharge();
518         if (GetDebug(1)) {
519           if (el<=0.0) cout<<"el="<<el<<" kwi["<<i<<"]="<<kwi[i]
520                           <<" de="<<de<<endl;
521         } // end if GetDebug
522         sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
523         sigx=sig;
524         sigz=sig*fda;
525         if (fLorentz) ld=(y+thick)*fTanLorAng;
526         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
527         //                cout << "sigx sigz " << sigx << " " << sigz << endl; // dom
528       } // end for i // End Integrate over t
529     else { // st == 0.0 deposit it at this point
530       x   = x0;
531       y   = y0;
532       z   = z0;
533       if (!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
534       el  = de / simpar->GetGeVToCharge();
535       sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
536       sigx=sig;
537       sigz=sig*fda;
538       if (fLorentz) ld=(y+thick)*fTanLorAng;
539       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
540     } // end if st>0.0
541     
542   } // Loop over all hits h
543   
544   // Coupling
545   switch (fCoupling) {
546   default:
547     break;
548   case 1: // case 3:
549     for (i=0;i<GetMap()->GetEntries();i++)
550       if (GetMap()->GetpListItem(i)==0) continue;
551       else{
552         GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
553         SetCoupling(iz,ix,idtrack,h);
554       } // end for i
555     break;
556   case 2: // case 4:
557     for (i=0;i<GetMap()->GetEntries();i++)
558       if (GetMap()->GetpListItem(i)==0) continue;
559       else {
560         GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);  
561         SetCouplingOld(iz,ix,idtrack,h);
562       } // end for i
563     break;
564   } // end switch
565   if (GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling);
566 }
567
568 //______________________________________________________________________
569 void AliITSsimulationPixUpg::SpreadCharge(Double_t x0,Double_t z0,
570                                           Int_t ix0,Int_t iz0,
571                                           Double_t el,Double_t sig,Double_t ld,
572                                           Int_t t,Int_t hi)
573 {
574    // Spreads the charge over neighboring cells. Assume charge is distributed
575    // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
576    // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
577    // if fLorentz=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
578    // Defined this way, the integral over all x and z is el.
579    // Inputs:
580    //    Double_t x0   x position of point where charge is liberated
581    //    Double_t z0   z position of point where charge is liberated
582    //    Int_t    ix0  row of cell corresponding to point x0
583    //    Int_t    iz0  columb of cell corresponding to point z0
584    //    Double_t el   number of electrons liberated in this step
585    //    Double_t sig  Sigma difusion for this step (y0 dependent)
586    //    Double_t ld   lorentz drift in x for this step (y0 dependent)
587    //    Int_t    t    track number
588    //    Int_t    ti   hit track index number
589    //    Int_t    hi   hit "hit" index number
590    // Outputs:
591    //     none.
592    // Return:
593    //     none.
594    const Int_t knx = 3,knz = 2;
595    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
596    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
597    Int_t ix,iz,ixs,ixe,izs,ize;
598    Float_t x,z;
599    Double_t x1,x2,z1,z2,s,sp;
600    AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(-1);
601    //
602    if (GetDebug(4)) Info("SpreadCharge","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi);
603    if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
604      GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
605      if (GetDebug(2)) cout << "sig<=0.0=" << sig << endl; // end if GetDebug
606      return;
607    } // end if
608    sp = 1.0/(sig*kRoot2);
609    if (GetDebug(2)) cout << "sig=" << sig << " sp=" << sp << endl; // end if GetDebug
610    ixs = TMath::Max(-knx+ix0,0);
611    ixe = TMath::Min(knx+ix0,seg->Npx()-1);
612    izs = TMath::Max(-knz+iz0,0);
613    ize = TMath::Min(knz+iz0,seg->Npz()-1);
614    for (ix=ixs;ix<=ixe;ix++) for (iz=izs;iz<=ize;iz++) {
615        seg->DetToLocal(ix,iz,x,z); // pixel center
616        x1  = x;
617        z1  = z;
618        x2  = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper
619        x1 -= 0.5*kmictocm*seg->Dpx(ix);  // Lower
620        z2  = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper
621        z1 -= 0.5*kmictocm*seg->Dpz(iz);  // Lower
622        x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
623        x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
624        z1 -= z0; // Distance from where track traveled
625        z2 -= z0; // Distance from where track traveled
626        s   = 0.25; // Correction based on definision of Erfc
627        s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
628        if (GetDebug(3)) {
629          cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
630            " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<< 
631            " sp*x1="<<sp*x1<<" sp*x2="<<sp*x2<<" s="<<s;
632        } // end if GetDebug
633        s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
634        if (GetDebug(3)) cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl; // end if GetDebug
635        GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
636      } // end for ix, iz
637 }
638
639 //______________________________________________________________________
640 void AliITSsimulationPixUpg::SpreadChargeAsym(Double_t x0,Double_t z0,
641                                            Int_t ix0,Int_t iz0,
642                                            Double_t el,Double_t sigx,Double_t sigz,
643                                            Double_t ld,Int_t t,Int_t hi)
644 {
645    // Spreads the charge over neighboring cells. Assume charge is distributed
646    // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
647    // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
648    // if fLorentz=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
649    // Defined this way, the integral over all x and z is el.
650    // Inputs:
651    //    Double_t x0   x position of point where charge is liberated
652    //    Double_t z0   z position of point where charge is liberated
653    //    Int_t    ix0  row of cell corresponding to point x0
654    //    Int_t    iz0  columb of cell corresponding to point z0
655    //    Double_t el   number of electrons liberated in this step
656    //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
657    //    Double_t sigz Sigma difusion along z for this step (y0 dependent)
658    //    Double_t ld   lorentz drift in x for this stip (y0 dependent)
659    //    Int_t    t    track number
660    //    Int_t    ti   hit track index number
661    //    Int_t    hi   hit "hit" index number
662    // Outputs:
663    //     none.
664    // Return:
665    //     none.
666    const Int_t knx = 3,knz = 2;
667    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
668    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
669    Int_t ix,iz,ixs,ixe,izs,ize;
670    Float_t x,z;
671    Double_t x1,x2,z1,z2,s,spx,spz;
672    AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(-1);
673    //
674    if (GetDebug(4)) Info("SpreadChargeAsym","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d)",x0,z0,ix0,iz0,el,sigx,sigz,t,hi);
675    if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
676      GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
677      if (GetDebug(2)) {
678        cout << "sigx<=0.0=" << sigx << endl;
679        cout << "sigz<=0.0=" << sigz << endl;
680      } // end if GetDebug
681      return;
682    } // end if
683    spx = 1.0/(sigx*kRoot2);     spz = 1.0/(sigz*kRoot2);
684    if (GetDebug(2)) {
685      cout << "sigx=" << sigx << " spx=" << spx << endl;
686      cout << "sigz=" << sigz << " spz=" << spz << endl;
687    } // end if GetDebug
688    ixs = TMath::Max(-knx+ix0,0);
689    ixe = TMath::Min(knx+ix0,seg->Npx()-1);
690    izs = TMath::Max(-knz+iz0,0);
691    ize = TMath::Min(knz+iz0,seg->Npz()-1);
692    for (ix=ixs;ix<=ixe;ix++) for (iz=izs;iz<=ize;iz++) {
693        seg->DetToLocal(ix,iz,x,z); // pixel center
694        x1  = x;
695        z1  = z;
696        x2  = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper
697        x1 -= 0.5*kmictocm*seg->Dpx(ix);  // Lower
698        z2  = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper
699        z1 -= 0.5*kmictocm*seg->Dpz(iz);  // Lower
700        x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
701        x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
702        z1 -= z0; // Distance from where track traveled
703        z2 -= z0; // Distance from where track traveled
704        s   = 0.25; // Correction based on definision of Erfc
705        s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
706        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
707        s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
708        if (GetDebug(3)) cout<<" spz*z1="<<spz*z1<<" spz*z2="<<spz*z2<<" s="<<s<< endl; // end if GetDebug
709        GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
710      } // end for ix, iz
711 }
712
713 //______________________________________________________________________
714 void AliITSsimulationPixUpg::RemoveDeadPixels() 
715 {
716   // Removes dead pixels on each module (ladder)
717   // This should be called before going from sdigits to digits (FrompListToDigits)
718   Int_t mod = GetModuleNumber();
719   AliITSCalibrationPixUpg* calObj = (AliITSCalibrationPixUpg*) fDetType->GetCalibrationModel(mod);
720
721   Int_t nrDead = calObj->GetNrBad();
722   for (Int_t i=0; i<nrDead; i++) GetMap()->DeleteHit(calObj->GetBadColAt(i), calObj->GetBadRowAt(i));
723   //
724 }
725
726 //______________________________________________________________________
727 void AliITSsimulationPixUpg::AddNoisyPixels() 
728 {
729   // Adds noisy pixels on each module (ladder)
730   // This should be called before going from sdigits to digits (FrompListToDigits)
731   Int_t mod = GetModuleNumber();
732   AliITSCalibrationPixUpg* calObj = (AliITSCalibrationPixUpg*) fDetType->GetPixNoisyModel(mod);
733   //
734   Int_t nrNoisy = calObj->GetNrBad();
735   for (Int_t i=0; i<nrNoisy; i++) {
736     // adding 10 times the threshold will for sure make this pixel fire...
737     GetMap()->AddNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), mod, 10*GetThreshold());
738   }
739 }
740
741 //______________________________________________________________________
742 void AliITSsimulationPixUpg::FrompListToDigits() 
743 {
744   // add noise and electronics, perform the zero suppression and add the
745   // digit to the list
746   // Inputs:
747   //    none.
748   // Outputs:
749   //    none.
750   // Return:
751   //    none.
752   static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
753   Int_t j,ix,iz;
754   Double_t  electronics;
755   Double_t sig;
756   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
757   static AliITSdigitSPD dig;
758   AliITSSimuParam *simpar = fDetType->GetSimuParam();
759   if (GetDebug(1)) Info("FrompListToDigits","()");
760   for (iz=0; iz<GetNPixelsZ(); iz++) 
761     for (ix=0; ix<GetNPixelsX(); ix++) {
762       electronics = simpar->ApplySPDBaselineAndNoise();
763       UpdateMapNoise(ix,iz,electronics);
764       //
765       // Apply Threshold and write Digits.
766       sig = GetMap()->GetSignalOnly(iz,ix);
767       FillHistograms(ix,iz,sig+electronics);
768       if (GetDebug(3)) {
769         cout<<sig<<"+"<<electronics<<">threshold("<<ix<<","<<iz
770                <<")="<<GetThreshold() <<endl;
771       } // end if GetDebug
772       // if (sig+electronics <= GetThreshold()) continue;
773       if (GetMap()->GetSignal(iz,ix) <= GetThreshold()) continue;
774       dig.SetCoord1(iz);
775       dig.SetCoord2(ix);
776       dig.SetSignal(1);
777       
778       //        dig.SetSignalSPD((Int_t) GetMap()->GetSignal(iz,ix));
779       Double_t aSignal =  GetMap()->GetSignal(iz,ix);
780       if (TMath::Abs(aSignal)>2147483647.0) {
781         //PH 2147483647 is the max. integer
782         //PH This apparently is a problem which needs investigation
783         AliWarning(Form("Too big or too small signal value %f",aSignal));
784         aSignal = TMath::Sign((Double_t)2147483647,aSignal);
785       }
786       dig.SetSignalSPD((Int_t)aSignal);
787       
788       for (j=0;j<knmaxtrk;j++) {
789         if (j<GetMap()->GetNEntries()) {
790           dig.SetTrack(j,GetMap()->GetTrack(iz,ix,j));
791           dig.SetHit(j,GetMap()->GetHit(iz,ix,j));
792         }else { // Default values
793           dig.SetTrack(j,-3);
794           dig.SetHit(j,-1);
795         } // end if GetMap()
796       } // end for j
797       if (GetDebug(3)) {
798         cout<<iz<<","<<ix<<","<<*(GetMap()->GetpListItem(iz,ix))<<endl;
799       } // end if GetDebug
800       aliITS->AddSimDigit(0,&dig);
801       // simulate fo signal response for this pixel hit:
802       //RS      fDetType->ProcessSPDDigitForFastOr(fModule, dig.GetCoord1(), dig.GetCoord2());
803     } //  for ix/iz
804 }
805
806 //______________________________________________________________________
807 void AliITSsimulationPixUpg::CreateHistograms() 
808 {
809    // create 1D histograms for tests
810    // Inputs:
811    //    none.
812    // Outputs:
813    //    none.
814    // Return:
815    //     none.
816
817    if (GetDebug(1)) Info("CreateHistograms","create histograms");
818
819    fHistos = new TObjArray(GetNPixelsZ());
820    fHistName="pix_";
821    for (Int_t i=0;i<GetNPixelsZ();i++) {
822      Char_t pixelz[4];
823      snprintf(pixelz,3,"%d",i);
824      fHistName.Append(pixelz);
825      fHistos->AddAt(new TH1F(fHistName.Data(),"pixel maps",
826                              GetNPixelsX(),0.,(Double_t)GetNPixelsX()),i);
827    } // end for i
828 }
829
830 //______________________________________________________________________
831 void AliITSsimulationPixUpg::FillHistograms(Int_t ix,Int_t iz,Double_t v) 
832 {
833    // Fill the histogram
834    // Inputs:
835    //    none.
836    // Outputs:
837    //    none.
838    // Return:
839    //     none.
840   
841   if (!GetHistArray()) return; // Only fill if setup.
842   if (GetDebug(2)) Info("FillHistograms","fill histograms");
843   GetHistogram(iz)->Fill(ix,v);
844 }
845
846 //______________________________________________________________________
847 void AliITSsimulationPixUpg::ResetHistograms() 
848 {
849   // Reset histograms for this detector
850   // Inputs:
851   //    none.
852   // Outputs:
853   //    none.
854   // Return:
855   //     none.
856   
857   if (!GetHistArray()) return; // Only fill if setup.
858   if (GetDebug(2)) Info("FillHistograms","fill histograms");
859   for ( int i=0;i<GetNPixelsZ();i++ ) if (fHistos->At(i))    ((TH1F*)fHistos->At(i))->Reset();
860   //
861 }
862
863 //______________________________________________________________________
864 void AliITSsimulationPixUpg::SetCoupling(Int_t col, Int_t row, Int_t ntrack, Int_t idhit) {
865   //  Take into account the coupling between adiacent pixels.
866   //  The parameters probcol and probrow are the probability of the
867   //  signal in one pixel shared in the two adjacent pixels along
868   //  the column and row direction, respectively.
869   //  Note pList is goten via GetMap() and module is not need any more.
870   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
871   //Begin_Html
872   /*
873     <img src="picts/ITS/barimodel_3.gif">
874      </pre>
875      <br clear=left>
876      <font size=+2 color=red>
877      <a href="mailto:tiziano.virgili@cern.ch"></a>.
878      </font>
879      <pre>
880    */
881    //End_Html
882    // Inputs:
883    //    Int_t col            z cell index
884    //    Int_t row            x cell index
885    //    Int_t ntrack         track incex number
886    //    Int_t idhit          hit index number
887    // Outputs:
888    //    none.
889    // Return:
890    //     none.
891    Int_t j1,j2,flag=0;
892    Double_t pulse1,pulse2;
893    Double_t couplR=0.0,couplC=0.0;
894    Double_t xr=0.;
895
896    GetCouplings(couplC,couplR);
897    if (GetDebug(3)) Info("SetCoupling","(col=%d,row=%d,ntrack=%d,idhit=%d) "
898                          "Calling SetCoupling couplC=%e couplR=%e",
899                          col,row,ntrack,idhit,couplC,couplR);
900    j1 = col;
901    j2 = row;
902    pulse1 = GetMap()->GetSignalOnly(col,row);
903    pulse2 = pulse1;
904    for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
905      do {
906        j1 += isign;
907        xr = gRandom->Rndm();
908        if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplC)) {
909          j1 = col;
910          flag = 1;
911        } else {
912          UpdateMapSignal(row,j1,ntrack,idhit,pulse1);
913          //  flag = 0;
914          flag = 1; // only first next!!
915        } // end if
916      } while (flag == 0);
917      // loop in row direction
918      do {
919        j2 += isign;
920        xr = gRandom->Rndm();
921        if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplR)) {
922          j2 = row;
923          flag = 1;
924        } else {
925          UpdateMapSignal(j2,col,ntrack,idhit,pulse2);
926          //  flag = 0;
927          flag = 1; // only first next!!
928        } // end if
929      } while(flag == 0);
930    } // for isign
931 }
932
933 //______________________________________________________________________
934 void AliITSsimulationPixUpg::SetCouplingOld(Int_t col, Int_t row,Int_t ntrack,Int_t idhit) 
935 {
936   //  Take into account the coupling between adiacent pixels.
937   //  The parameters probcol and probrow are the fractions of the
938   //  signal in one pixel shared in the two adjacent pixels along
939   //  the column and row direction, respectively.
940   //Begin_Html
941   /*
942     <img src="picts/ITS/barimodel_3.gif">
943     </pre>
944     <br clear=left>
945     <font size=+2 color=red>
946     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
947     </font>
948     <pre>
949   */
950   //End_Html
951   // Inputs:
952   //    Int_t col            z cell index
953   //    Int_t row            x cell index
954   //    Int_t ntrack         track incex number
955   //    Int_t idhit          hit index number
956   //    Int_t module         module number
957   // Outputs:
958   //    none.
959   // Return:
960   //     none.
961   Int_t j1,j2,flag=0;
962   Double_t pulse1,pulse2;
963   Double_t couplR=0.0,couplC=0.0;
964   
965   GetCouplings(couplC,couplR);
966   
967   //  Debugging ...
968   //    cout << "Threshold --> " << GetThreshold() << endl;  // dom
969   //    cout << "Couplings --> " << couplC << " " << couplR << endl;  // dom
970   
971   if (GetDebug(3)) Info("SetCouplingOld","(col=%d,row=%d,ntrack=%d,idhit=%d) "
972                         "Calling SetCoupling couplC=%e couplR=%e",
973                         col,row,ntrack,idhit,couplC,couplR);
974   for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
975     pulse1 = GetMap()->GetSignalOnly(col,row);
976     pulse2 = pulse1;
977     j1 = col;
978     j2 = row;
979     do{
980       j1 += isign;
981       pulse1 *= couplC;
982       if ((j1<0)||(j1>GetNPixelsZ()-1)||(pulse1<GetThreshold())) {
983         pulse1 = GetMap()->GetSignalOnly(col,row);
984         j1 = col;
985         flag = 1;
986       }else{
987         UpdateMapSignal(row,j1,ntrack,idhit,pulse1);
988         // flag = 0;
989         flag = 1;  // only first next !!
990       } // end if
991     } while(flag == 0);
992     // loop in row direction
993     do{
994       j2 += isign;
995       pulse2 *= couplR;
996       if ((j2<0)||(j2>(GetNPixelsX()-1))||(pulse2<GetThreshold())) {
997         pulse2 = GetMap()->GetSignalOnly(col,row);
998         j2 = row;
999         flag = 1;
1000       }else{
1001         UpdateMapSignal(j2,col,ntrack,idhit,pulse2);
1002         // flag = 0;
1003         flag = 1; // only first next!!
1004       } // end if
1005     } while(flag == 0);
1006   } // for isign
1007 }
1008
1009 //______________________________________________________________________
1010 void AliITSsimulationPixUpg::GenerateStrobePhase()
1011 {
1012   // Generate randomly the strobe
1013   // phase w.r.t to the LHC clock
1014   // Done once per event
1015   const Double_t kBunchLenght = 25e-9; // LHC clock
1016   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
1017     (Double_t)fStrobeLenght*kBunchLenght+
1018     kBunchLenght/2;
1019 }
1020