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