Introducing Riostream.h
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPDdubna.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 $Log$
18 Revision 1.8  2002/10/14 14:57:08  hristov
19 Merging the VirtualMC branch to the main development branch (HEAD)
20
21 Revision 1.3.8.2  2002/10/14 13:14:08  hristov
22 Updating VirtualMC to v3-09-02
23
24 Revision 1.7  2002/09/09 17:23:28  nilsen
25 Minor changes in support of changes to AliITSdigitS?D class'.
26
27 Revision 1.6  2002/08/21 22:09:58  nilsen
28 Updated SPD simulation with difusion effects. ReWritten Hit to SDigits
29 code.
30
31 */
32 #include <Riostream.h>
33 #include <TRandom.h>
34 #include <TH1.h>
35 #include <TMath.h>
36 #include <TString.h>
37 #include <TParticle.h>
38
39 #include "AliRun.h"
40 #include "AliITS.h"
41 #include "AliITShit.h"
42 #include "AliITSdigit.h"
43 #include "AliITSmodule.h"
44 #include "AliITSMapA2.h" 
45 #include "AliITSpList.h"
46 #include "AliITSsimulationSPDdubna.h"
47 #include "AliITSsegmentationSPD.h"
48 #include "AliITSresponseSPDdubna.h"
49
50 //#define DEBUG
51
52 ClassImp(AliITSsimulationSPDdubna)
53 ////////////////////////////////////////////////////////////////////////
54 // Version: 0
55 // Written by Boris Batyunya
56 // December 20 1999
57 //
58 // AliITSsimulationSPDdubna is the simulation of SPDs
59 //______________________________________________________________________
60
61
62 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(){
63     // constructor
64
65     fResponse = 0;
66     fSegmentation = 0;
67     fMapA2 = 0;
68     fpList = 0;
69     fModule = 0;
70     fEvent = 0;
71     fHis = 0;
72     fNoise = 0.;
73     fBaseline = 0.;
74     fNPixelsZ = 0;
75     fNPixelsX = 0;
76 }
77 //______________________________________________________________________
78 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg,
79                                                    AliITSresponse *resp){
80     // standard constructor
81     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
82
83     fHis = 0;
84     fResponse = resp;
85     fSegmentation = seg;
86     fModule = 0;
87     fEvent = 0;
88
89     fNPixelsZ=fSegmentation->Npz();
90     fNPixelsX=fSegmentation->Npx();
91
92     fResponse->GetNoiseParam(fNoise,fBaseline);
93     fResponse->SetDistanceOverVoltage(kmictocm*fSegmentation->Dy(),50.0);
94
95 //    fMapA2 = new AliITSMapA2(fSegmentation);
96     fMapA2 = 0;
97
98     fpList = new AliITSpList(fNPixelsZ+1,fNPixelsX+1);
99
100 }
101 //______________________________________________________________________
102 AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna(){
103     // destructor
104
105     if(fMapA2) delete fMapA2;
106
107     if (fHis) {
108         fHis->Delete(); 
109         delete fHis;     
110     } // end if fHis
111 }
112 //______________________________________________________________________
113 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const 
114                                                    AliITSsimulationSPDdubna 
115                                                    &source){
116     //     Copy Constructor 
117     if(&source == this) return;
118     this->fMapA2 = source.fMapA2;
119     this->fNoise = source.fNoise;
120     this->fBaseline = source.fBaseline;
121     this->fNPixelsX = source.fNPixelsX;
122     this->fNPixelsZ = source.fNPixelsZ;
123     this->fHis = source.fHis;
124     return;
125 }
126 //______________________________________________________________________
127 AliITSsimulationSPDdubna&  AliITSsimulationSPDdubna::operator=(const 
128                                            AliITSsimulationSPDdubna &source){
129     //    Assignment operator
130     if(&source == this) return *this;
131     this->fMapA2 = source.fMapA2;
132     this->fNoise = source.fNoise;
133     this->fBaseline = source.fBaseline;
134     this->fNPixelsX = source.fNPixelsX;
135     this->fNPixelsZ = source.fNPixelsZ;
136     this->fHis = source.fHis;
137     return *this;
138 }
139 //______________________________________________________________________
140 void AliITSsimulationSPDdubna::InitSimulationModule(Int_t module, Int_t event){
141     //  This function creates maps to build the list of tracks for each
142     //  summable digit.
143     //
144     //  Inputs:
145     //    Int_t module   // Module number to be simulated
146     //    Int_t event    // Event number to be simulated
147     //
148     //  Outputs:
149     //    none
150     //
151     //  Returns:
152     //    none
153
154     fModule = module;
155     fEvent  = event;
156 //    fMapA2->ClearMap();
157     fpList->ClearMap();
158 }
159 //_____________________________________________________________________
160 void AliITSsimulationSPDdubna::SDigitiseModule(AliITSmodule *mod, Int_t mask,
161                                                Int_t event){
162     //  This function begins the work of creating S-Digits
163     //
164     //  Inputs:
165     //    AliITSmodule *mod  //  module
166     //    Int_t mask         //  mask to be applied to the module
167     //
168     //  Outputs:
169     //    none
170     //
171     //  Return:
172     //    test              //  test returns kTRUE if the module contained hits
173     //                      //  test returns kFALSE if it did not contain hits
174
175     Int_t module = 0;
176
177     if(!(mod->GetNhits())) return;// if module has no hits don't create Sdigits
178     fModule = mod->GetIndex();
179     HitToSDigit(mod, module, mask, fpList);
180     WriteSDigits(fpList);
181 //    fMapA2->ClearMap();
182     fpList->ClearMap();
183 }
184 //______________________________________________________________________
185 void AliITSsimulationSPDdubna::WriteSDigits(AliITSpList *pList){
186     //  This function adds each S-Digit to pList
187     //
188     //  Inputs:
189     //    AliITSpList *pList
190     //
191     //  Outputs:
192     //    none
193     //
194     //  Return:
195     //    none
196     Int_t ix, nix, iz, niz;
197     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
198
199     pList->GetMaxMapIndex(niz, nix);
200     for(iz=0; iz<niz; iz++)for(ix=0; ix<nix; ix++){
201         if(pList->GetSignalOnly(iz,ix)>0.0){
202             aliITS->AddSumDigit(*(pList->GetpListItem(iz,ix)));
203 #ifdef DEBUG
204             cout <<"SDigits " << iz << "," << ix << "," << 
205                 *(pList->GetpListItem(iz,ix)) << endl;
206 #endif
207         } // end if pList
208     } // end for iz,ix
209     return; 
210 }
211 //______________________________________________________________________
212 void AliITSsimulationSPDdubna::FinishSDigitiseModule(){
213     //  This function calls SDigitsToDigits which creates Digits from SDigits
214     //
215     //  Inputs:
216     //    none
217     //
218     //  Outputs:
219     //    none
220     //  Return
221     //    none
222
223     SDigitsToDigits(fModule, fpList);
224     return;
225 }
226 //______________________________________________________________________
227 void AliITSsimulationSPDdubna::SDigitsToDigits(Int_t module,
228                                                AliITSpList *pList){
229     //  This function adds electronic noise to the S-Digits and then adds them
230     // to  a new pList
231     //
232     //  Inputs:
233     //    Int_t       module  // module number
234     //    AliITSpList *pList  // pList
235     //
236     //  Outputs:
237     //    pList is passed along to the functions ChargeToSignal and GetList
238     //
239     //  Return:
240     //    none
241
242     fModule = module;
243     ChargeToSignal(pList); // Charge To Signal both adds noise and
244 //    fMapA2->ClearMap();
245     pList->ClearMap();
246 }
247 //______________________________________________________________________
248 void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module,
249                                               Int_t dummy){
250     //  This function creates Digits straight from the hits and then adds
251     //  electronic noise to the digits before adding them to pList
252     //
253     //  Inputs:
254     //    AliITSmodule *mod    // module
255     //    Int_t        module  // module number  Dummy.
256     //    Int_t        dummy
257     //
258     //  Outputs:
259     //    Each of the input variables is passed along to HitToSDigit
260     //
261     //  Return:
262     //    none
263
264     fModule = mod->GetIndex();  //This calls the module for HitToSDigit
265     HitToSDigit(mod,fModule, dummy, fpList);
266     ChargeToSignal(fpList);
267 //    fMapA2->ClearMap();
268     fpList->ClearMap();
269 }
270 //______________________________________________________________________
271 void AliITSsimulationSPDdubna::UpdateMapSignal(Int_t iz, Int_t ix, Int_t trk,
272                                                Int_t ht, Int_t module,
273                                                Double_t signal,
274                                                AliITSpList *pList){
275     //  This function adds a signal to the pList from the pList class
276     //
277     //  Inputs:
278     //    Int_t       iz     // row number
279     //    Int_t       ix     // column number
280     //    Int_t       trk    // track number
281     //    Int_t       ht     // hit number
282     //    Double_t    signal // signal strength
283     //    AliITSpList *pList // pList
284     //
285     //  Outputs:
286     //    All of the inputs are passed to AliITSpList::AddSignal
287     //    Int_t    ix  // row number
288     //    Int_t    iz  // column number
289     //    Double_t sig // signal strength
290     //          // These three variables are defined to preserve the
291     //          // assignments used in the function AliITSMapA2::AddSignal
292     //
293     //  Return:
294     //    none
295
296 //    fMapA2->AddSignal(iz, ix, signal);
297     pList->AddSignal(iz,ix, trk, ht, fModule, signal);
298 }
299 //______________________________________________________________________
300 void AliITSsimulationSPDdubna::UpdateMapNoise(Int_t iz,
301                                               Int_t ix, Int_t fModule,
302                                               Double_t sig, Float_t noise,
303                                               AliITSpList *pList){
304     //  This function adds noise to data in the MapA2 as well as the pList
305     //
306     //  Inputs:
307     //    Int_t       iz       // row number
308     //    Int_t       ix       // column number
309     //    Int_t       mod     // module number
310     //    Double_t    sig     // signal strength
311     //    Double_t    noise   // electronic noise generated by ChargeToSignal
312     //    AliITSpList *pList  // pList
313     //
314     //  Outputs:
315     //    All of the inputs are passed to AliITSMapA2::AddSignal or
316     //    AliITSpList::AddNoise
317     //
318     //  Return:
319     //    none
320
321 //    fMapA2->AddSignal(iz, ix, noise);
322     pList->AddNoise(iz,ix, fModule, noise);
323 }
324 //______________________________________________________________________
325 void AliITSsimulationSPDdubna::HitToDigit(AliITSmodule *mod, Int_t module,
326                                           Int_t dummy){
327     DigitiseModule(mod, module, dummy);
328 }
329 //______________________________________________________________________
330 void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod, Int_t module,
331                                             Int_t dummy,AliITSpList *pList){
332     // Does the charge distributions using Gaussian diffusion charge charing.
333     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
334     TObjArray *hits = mod->GetHits();
335     Int_t nhits = hits->GetEntriesFast();
336     Int_t h,ix,iz;
337     Int_t idtrack;
338     Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
339     Double_t x,y,z,t,tp,st,dt=0.2,el,sig;
340     Double_t thick = kmictocm*GetSeg()->Dy();
341
342     if(nhits<=0) return;
343     for(h=0;h<nhits;h++){
344 #ifdef DEBUG
345         cout << "Hits=" << h << "," << *(mod->GetHit(h)) << endl;
346 #endif
347         if(mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)){
348         st =TMath::Sqrt(x1*x1+y1*y1+z1*z1);
349         if(st>0.0) for(t=0;t<1.0;t+=dt){ // Integrate over t
350             tp = t+0.5*dt;
351             el = GetResp()->GeVToCharge((Float_t)(dt*de));
352 #ifdef DEBUG
353             if(el<=0.0) cout << "el="<<el<<" dt="<<dt<<" de="<<de<<endl;
354 #endif
355             x = x0+x1*tp;
356             y = y0+y1*tp;
357             z = z0+z1*tp;
358             GetSeg()->LocalToDet(x,z,ix,iz);
359             sig = GetResp()->SigmaDiffusion1D(thick + y);
360             SpreadCharge(x,y,z,ix,iz,el,sig,
361                          idtrack,mod->GetHitTrackIndex(h),h,mod->GetIndex());
362         } else { // st == 0.0 deposit it at this point
363             el = GetResp()->GeVToCharge((Float_t)de);
364             x = x0;
365             y = y0;
366             z = z0;
367             GetSeg()->LocalToDet(x,z,ix,iz);
368             sig = GetResp()->SigmaDiffusion1D(thick + y);
369             SpreadCharge(x,y,z,ix,iz,el,sig,
370                          idtrack,mod->GetHitTrackIndex(h),h,mod->GetIndex());
371         } // end if st>0.0
372     }} // Loop over all hits h
373 }
374 //______________________________________________________________________
375 void AliITSsimulationSPDdubna::SpreadCharge(Double_t x0,Double_t y0,
376                                             Double_t z0,Int_t ix0,Int_t iz0,
377                                             Double_t el,Double_t sig,Int_t t,
378                                             Int_t ti,Int_t hi,Int_t mod){
379     // Spreads the charge over neighboring cells. Assume charge is distributed
380     // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
381     // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
382     // Defined this way, the integral over all x and z is el.
383     const Int_t knx = 3,knz = 2;
384     const Double_t kRoot2 = 1.414213562; // Sqrt(2).
385     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
386     Int_t ix,iz,ixs,ixe,izs,ize;
387     Float_t x,z;
388     Double_t x1,x2,z1,z2,s,sp;
389
390     if(sig<=0.0) {
391         fpList->AddSignal(iz0,ix0,t,hi,mod,el);
392         return;
393     } // end if
394     sp = 1.0/(sig*kRoot2);
395 #ifdef DEBUG
396     cout << "sig=" << sig << " sp=" << sp << endl;
397 #endif
398     ixs = TMath::Max(-knx+ix0,0);
399     ixe = TMath::Min(knx+ix0,GetSeg()->Npx()-1);
400     izs = TMath::Max(-knz+iz0,0);
401     ize = TMath::Min(knz+iz0,GetSeg()->Npz()-1);
402     for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
403         GetSeg()->DetToLocal(ix,iz,x,z); // pixel center
404         x1 = x;
405         z1 = z;
406         x2  = x1 + 0.5*kmictocm*GetSeg()->Dpx(ix); // Upper
407         x1 -= 0.5*kmictocm*GetSeg()->Dpx(ix);  // Lower
408         z2  = z1 + 0.5*kmictocm*GetSeg()->Dpz(iz); // Upper
409         z1 -= 0.5*kmictocm*GetSeg()->Dpz(iz);  // Lower
410         x1 -= x0; // Distance from where track traveled
411         x2 -= x0; // Distance from where track traveled
412         z1 -= z0; // Distance from where track traveled
413         z2 -= z0; // Distance from where track traveled
414         s = 0.25; // Correction based on definision of Erfc
415         s *= TMath::Erfc(sp*x1) - TMath::Erfc(sp*x2);
416 #ifdef DEBUG
417         cout << "el=" << el << " ix0=" << ix0 << " ix=" << ix << " x0="<< x <<
418             " iz0=" << iz0 << " iz=" << iz << " z0=" << z  << 
419             " sp*x1=" << sp*x1 <<" sp*x2=" << sp*x2 << " s=" << s;
420 #endif
421         s *= TMath::Erfc(sp*z1) - TMath::Erfc(sp*z2);
422 #ifdef DEBUG
423         cout << " sp*z1=" << sp*z1 <<" sp*z2=" << sp*z2 << " s=" << s << endl;
424 #endif
425         fpList->AddSignal(iz,ix,t,hi,mod,s*el);
426     } // end for ix, iz
427 }
428 //______________________________________________________________________
429 void AliITSsimulationSPDdubna::HitToSDigitOld(AliITSmodule *mod, Int_t module,
430                                            Int_t dummy, AliITSpList *pList){
431     // digitize module 
432     const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons 
433                                       // for 3.6 eV/pair 
434     const Float_t kconv = 10000.;     // cm -> microns
435
436     Float_t spdLength = fSegmentation->Dz();
437     Float_t spdWidth = fSegmentation->Dx();
438     Float_t spdThickness = fSegmentation->Dy();
439     Float_t difCoef, dum;       
440     fResponse->DiffCoeff(difCoef,dum); 
441     if(spdThickness > 290) difCoef = 0.00613;  
442
443     Float_t zPix0 = 1e+6;
444     Float_t xPix0 = 1e+6;
445     Float_t yPrev = 1e+6;   
446
447     Float_t zPitch = fSegmentation->Dpz(0);
448     Float_t xPitch = fSegmentation->Dpx(0);
449   
450     TObjArray *fHits = mod->GetHits();
451     module = mod->GetIndex();
452     Int_t nhits = fHits->GetEntriesFast();
453     if (!nhits) return;
454 #ifdef DEBUG
455     cout<<"len,wid,thickness,nx,nz,pitchx,pitchz,difcoef ="<<spdLength<<","
456         <<spdWidth<<","<<spdThickness<<","<<fNPixelsX<<","<<fNPixelsZ<<","
457         <<xPitch<<","<<zPitch<<","<<difCoef<<endl;
458 #endif
459     //  Array of pointers to the label-signal list
460     Int_t indexRange[4] = {0,0,0,0};
461
462     // Fill detector maps with GEANT hits
463     // loop over hits in the module
464     static Bool_t first;
465     Int_t lasttrack=-2;
466     Int_t hit, iZi, jz, jx;
467     Int_t idhit=-1; //!
468 #ifdef DEBUG
469     cout<<"SPDdubna: module,nhits ="<<module<<","<<nhits<<endl;
470 #endif
471     for (hit=0;hit<nhits;hit++) {
472         AliITShit *iHit = (AliITShit*) fHits->At(hit);
473 #ifdef DEBUG
474         cout << "Hits=" << hit << "," << *iHit << endl;
475 #endif
476         //Int_t layer = iHit->GetLayer();
477         Float_t yPix0 = -spdThickness/2; 
478
479         // work with the idtrack=entry number in the TreeH
480         //Int_t idhit,idtrack; //!
481         //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);  //!    
482         //Int_t idtrack=mod->GetHitTrackIndex(hit);  
483         // or store straight away the particle position in the array
484         // of particles : 
485         if(iHit->StatusEntering()) idhit=hit;
486         Int_t itrack = iHit->GetTrack();
487         Int_t dray = 0;
488    
489         if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE; 
490
491         //Int_t parent = iHit->GetParticle()->GetFirstMother();
492         Int_t partcode = iHit->GetParticle()->GetPdgCode();
493
494         //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0,
495         // 211 - pi+,  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
496
497         Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
498                                                    // vertex
499         pmod *= 1000;
500
501         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
502                                                  // at p < 6 MeV/c
503
504         //  Get hit z and x(r*phi) cordinates for each module (detector)
505         //  in local system.
506
507         Float_t zPix = kconv*iHit->GetZL();
508         Float_t xPix = kconv*iHit->GetXL();
509         Float_t yPix = kconv*iHit->GetYL();
510
511         // Get track status
512         Int_t status = iHit->GetTrackStatus();      
513
514         // Check boundaries
515         if(zPix  > spdLength/2) {
516 #ifdef DEBUG
517             cout<<"!!! SPD: z outside ="<<zPix<<endl;
518 #endif
519             zPix = spdLength/2 - 10;
520         }
521         if(zPix  < 0 && zPix < -spdLength/2) {
522 #ifdef DEBUG
523             cout<<"!!! SPD: z outside ="<<zPix<<endl;
524 #endif
525             zPix = -spdLength/2 + 10;
526         }
527         if(xPix  > spdWidth/2) {
528 #ifdef DEBUG
529             cout<<"!!! SPD: x outside ="<<xPix<<endl;
530 #endif
531             xPix = spdWidth/2 - 10;
532         }
533         if(xPix  < 0 && xPix < -spdWidth/2) {
534 #ifdef DEBUG
535             cout<<"!!! SPD: x outside ="<<xPix<<endl;
536 #endif
537             xPix = -spdWidth/2 + 10;
538         }
539         Int_t trdown = 0;
540
541         // enter Si or after event in Si
542         if (status == 66 ) {  
543             zPix0 = zPix;
544             xPix0 = xPix;
545             yPrev = yPix; 
546         } // end if status == 66
547
548         Float_t depEnergy = iHit->GetIonization();
549         // skip if the input point to Si       
550
551         if(depEnergy <= 0.) continue;        
552
553         // if track returns to the opposite direction:
554         if (yPix < yPrev) {
555             trdown = 1;
556         } // end if yPix < yPrev
557
558         // take into account the holes diffusion inside the Silicon
559         // the straight line between the entrance and exit points in Si is
560         // divided into the several steps; the diffusion is considered 
561         // for each end point of step and charge
562         // is distributed between the pixels through the diffusion.
563
564         //  ---------- the diffusion in Z (beam) direction -------
565         Float_t charge = depEnergy*kEnToEl;         // charge in e-
566         Float_t drPath = 0.;   
567         Float_t tang = 0.;
568         Float_t sigmaDif = 0.; 
569         Float_t zdif = zPix - zPix0;
570         Float_t xdif = xPix - xPix0;
571         Float_t ydif = TMath::Abs(yPix - yPrev);
572         Float_t ydif0 = TMath::Abs(yPrev - yPix0);
573
574         if(ydif < 1) continue; // ydif is not zero
575
576         Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
577
578         Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
579         Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1; 
580
581         // number of the steps along the track:
582         Int_t nsteps = ndZ;
583         if(ndX > ndZ) nsteps = ndX;
584         if(nsteps < 20) nsteps = 20;  // minimum number of the steps 
585
586         if (projDif < 5 ) {
587             drPath = (yPix-yPix0)*1.e-4;  
588             drPath = TMath::Abs(drPath);        // drift path in cm
589             sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm
590             sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
591             nsteps = 1;
592         }  // end if projDif < 5
593
594         if(projDif > 5) tang = ydif/projDif;
595         Float_t dCharge = charge/nsteps;       // charge in e- for one step
596         Float_t dZ = zdif/nsteps;
597         Float_t dX = xdif/nsteps;
598
599         for (iZi = 1; iZi <= nsteps;iZi++) {
600             Float_t dZn = iZi*dZ;
601             Float_t dXn = iZi*dX;
602             Float_t zPixn = zPix0 + dZn;
603             Float_t xPixn = xPix0 + dXn;
604
605             if(projDif >= 5) {
606                 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
607                 drPath = dProjn*tang*1.e-4; // drift path for iZi+1 step in cm 
608                 if(trdown == 0) {
609                     drPath = TMath::Abs(drPath) + ydif0*1.e-4;
610                 }// end if trdow ==0
611                 if(trdown == 1) {
612                     drPath = ydif0*1.e-4 - TMath::Abs(drPath);
613                     drPath = TMath::Abs(drPath);
614                 } // end if trdown == 1
615                 sigmaDif = difCoef*sqrt(drPath);    
616                 sigmaDif = sigmaDif*kconv;       // sigma diffusion in microns
617             } // end if projdif >= 5
618
619             zPixn = (zPixn + spdLength/2.);  
620             xPixn = (xPixn + spdWidth/2.);  
621             Int_t nZpix, nXpix;
622             fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
623             zPitch = fSegmentation->Dpz(nZpix);
624             fSegmentation->GetPadTxz(xPixn,zPixn);
625             // set the window for the integration
626             Int_t jzmin = 1;  
627             Int_t jzmax = 3; 
628             if(nZpix == 1) jzmin =2;
629             if(nZpix == fNPixelsZ) jzmax = 2; 
630
631             Int_t jxmin = 1;  
632             Int_t jxmax = 3; 
633             if(nXpix == 1) jxmin =2;
634             if(nXpix == fNPixelsX) jxmax = 2; 
635
636             Float_t zpix = nZpix; 
637             Float_t dZright = zPitch*(zpix - zPixn);
638             Float_t dZleft = zPitch - dZright;
639
640             Float_t xpix = nXpix; 
641             Float_t dXright = xPitch*(xpix - xPixn);
642             Float_t dXleft = xPitch - dXright;
643
644             Float_t dZprev = 0.;
645             Float_t dZnext = 0.;
646             Float_t dXprev = 0.;
647             Float_t dXnext = 0.;
648
649             for(jz=jzmin; jz <=jzmax; jz++) {
650                 if(jz == 1) {
651                     dZprev = -zPitch - dZleft;
652                     dZnext = -dZleft;
653                 } else if(jz == 2) {
654                     dZprev = -dZleft;
655                     dZnext = dZright;
656                 } else if(jz == 3) {
657                     dZprev = dZright;
658                     dZnext = dZright + zPitch;
659                 } // end if jz
660                 // kz changes from 1 to the fNofPixels(270)  
661                 Int_t kz = nZpix + jz -2; 
662
663                 Float_t zArg1 = dZprev/sigmaDif;
664                 Float_t zArg2 = dZnext/sigmaDif;
665                 Float_t zProb1 = TMath::Erfc(zArg1);
666                 Float_t zProb2 = TMath::Erfc(zArg2);
667                 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge; 
668
669
670                 // ----------- holes diffusion in X(r*phi) direction  --------
671
672                 if(dZCharge > 1.) { 
673                     for(jx=jxmin; jx <=jxmax; jx++) {
674                         if(jx == 1) {
675                             dXprev = -xPitch - dXleft;
676                             dXnext = -dXleft;
677                         } else if(jx == 2) {
678                             dXprev = -dXleft;
679                             dXnext = dXright;
680                         } else if(jx == 3) {
681                             dXprev = dXright;
682                             dXnext = dXright + xPitch;
683                         }  // end if jx
684                         Int_t kx = nXpix + jx -2;
685                         Float_t xArg1 = dXprev/sigmaDif;
686                         Float_t xArg2 = dXnext/sigmaDif;
687                         Float_t xProb1 = TMath::Erfc(xArg1);
688                         Float_t xProb2 = TMath::Erfc(xArg2);
689                         Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
690
691                         if(dXCharge > 1.) {
692                             if (first) {
693                                 indexRange[0]=indexRange[1]=kz-1;
694                                 indexRange[2]=indexRange[3]=kx-1;
695                                 first=kFALSE;
696                             } // end if first
697                             indexRange[0]=TMath::Min(indexRange[0],kz-1);
698                             indexRange[1]=TMath::Max(indexRange[1],kz-1);
699                             indexRange[2]=TMath::Min(indexRange[2],kx-1);
700                             indexRange[3]=TMath::Max(indexRange[3],kx-1);
701 /*
702                             // build the list of digits for this module 
703                             Double_t signal = fMapA2->GetSignal(kz-1,kx-1);
704                             signal+=dXCharge;
705                             fMapA2->SetHit(kz-1,kx-1,(double)signal);
706 */
707                             // The calling sequence for UpdateMapSignal was 
708                             // moved into the (dx > 1 e-) loop because it 
709                             // needs to call signal which is defined inside 
710                             // this loop
711                             fModule   = module;//Defined because functions 
712                                                // called by UpdateMapSignal 
713                                                // expect module to be an 
714                                                // integer
715                             UpdateMapSignal(kz-1,kx-1,
716 //                                          mod->GetHitTrackIndex(hit),
717                              ((AliITShit*)(mod->GetHit(hit)))->GetTrack(),
718                                             hit,fModule,dXCharge,pList);
719                         }      // dXCharge > 1 e-
720                     }       // jx loop
721                 }       // dZCharge > 1 e-
722             }        // jz loop
723         }         // iZi loop
724         if (status == 65) {   // the step is inside of Si
725             zPix0 = zPix;
726             xPix0 = xPix;
727         } // end if status == 65
728         yPrev = yPix;
729     }   // hit loop inside the module
730 }
731 //______________________________________________________________________
732 void AliITSsimulationSPDdubna::ChargeToSignal(AliITSpList *pList){
733     // add noise and electronics, perform the zero suppression and add the
734     // digit to the list
735     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
736     Float_t threshold = (float)fResponse->MinVal();
737     Int_t j;
738 //    Int_t    digits[3], tracks[3], hits[3];
739 //    Float_t  charges[3];
740     Float_t  electronics;
741 //    Float_t  phys; 
742     Double_t sig;
743     const Int_t    nmaxtrk=AliITSdigitSPD::GetNTracks();
744     static AliITSdigitSPD dig;
745
746     for(Int_t iz=0; iz<fNPixelsZ; iz++){
747         for(Int_t ix=0; ix<fNPixelsX; ix++){
748             electronics = fBaseline + fNoise*gRandom->Gaus();
749             sig = pList->GetSignalOnly(iz,ix);
750             UpdateMapNoise(iz,ix,fModule,sig,electronics,pList);
751 #ifdef DEBUG
752 //          cout << sig << "+" << electronics <<">threshold=" << threshold 
753 //               << endl;
754 #endif
755             if (sig+electronics > threshold) {
756                 dig.fCoord1 = iz;
757                 dig.fCoord2 = ix;
758                 dig.fSignal = 1;
759                 dig.fSignalSPD = (Int_t) pList->GetSignal(iz,ix);
760                 /*
761                 digits[0] = iz;
762                 digits[1] = ix;
763                 digits[2] = 1; */
764                 for(j=0;j<nmaxtrk;j++){
765 //                  charges[j] = 0.0;
766                     if (j<pList->GetNEnteries()) {
767                         dig.fTracks[j] = pList->GetTrack(iz,ix,j);
768                         dig.fHits[j]   = pList->GetHit(iz,ix,j);
769                         /*
770                         tracks[j] = pList->GetTrack(iz,ix,j);
771                         hits[j]   = pList->GetHit(iz,ix,j);
772                         */
773                     }else { // Default values
774                         dig.fTracks[j] = -3;
775                         dig.fHits[j]   = -1;
776 /*                      tracks[j] = -2; //noise
777                         hits[j]   = -1;  */
778                     } // end if pList
779                 } // end for j
780 //              charges[0] = (Float_t) pList->GetSumSignal(iz,ix);
781 /*
782                 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
783                     tracks[1] = -3;
784                     hits[1] = -1;
785                     tracks[2] = -3;
786                     hits[2] = -1;
787                 } else if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
788                     tracks[1] = -3;
789                     hits[1] = -1;   
790                 } else if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
791                     tracks[2] = -3;
792                     hits[2] = -1;   
793                 } else if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
794                     tracks[2] = -3;
795                     hits[2] = -1;   
796                 } // end if
797 */
798 //              phys = 0.0;
799 #ifdef DEBUG
800                 cout << iz << "," << ix << "," << 
801                     *(pList->GetpListItem(iz,ix)) << endl;
802 #endif
803 //              aliITS->AddSimDigit(0, phys, digits, tracks, hits, charges);
804                 aliITS->AddSimDigit(0,&dig);
805             } // 
806         } // 
807     } //
808 }
809 //______________________________________________________________________
810 void AliITSsimulationSPDdubna::CreateHistograms(){
811     // create 1D histograms for tests
812
813     printf("SPD - create histograms\n");
814
815     fHis=new TObjArray(fNPixelsZ);
816     TString spdName("spd_");
817     for (Int_t i=0;i<fNPixelsZ;i++) {
818         Char_t pixelz[4];
819         sprintf(pixelz,"%d",i+1);
820         spdName.Append(pixelz);
821         //PH       (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
822         //PH                       fNPixelsX,0.,(Float_t) fNPixelsX);
823         fHis->AddAt(new TH1F(spdName.Data(),"SPD maps",
824                              fNPixelsX,0.,(Float_t) fNPixelsX), i);
825     } // end for i
826 }
827 //______________________________________________________________________
828 void AliITSsimulationSPDdubna::ResetHistograms(){
829     //
830     // Reset histograms for this detector
831     //
832
833     for ( int i=0;i<fNPixelsZ;i++ ) {
834         //PH    if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
835         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
836     } // end for i
837 }