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