Fixes for Coverity warnings (M. van Leeuwen)
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSDD.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 /* $Id$ */
17
18 #include <Riostream.h>
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include <cstring>
22
23 #include <TCanvas.h>
24 #include <TF1.h>
25 #include <TH1.h>
26 #include <TFile.h>
27 #include <TRandom.h>
28 #include <TROOT.h>
29 #include "AliITS.h"
30 #include "AliITSMapA2.h"
31 #include "AliITSRawData.h"
32 #include "AliITSdigitSPD.h"
33 #include "AliITSetfSDD.h"
34 #include "AliITSmodule.h"
35 #include "AliITShit.h"
36 #include "AliITSpList.h"
37 #include "AliITSCalibrationSDD.h"
38 #include "AliITSresponseSDD.h"
39 #include "AliITSsimulationSDD.h"
40 #include "AliLog.h"
41 #include "AliRun.h"
42
43 ClassImp(AliITSsimulationSDD)
44 ////////////////////////////////////////////////////////////////////////
45 // Version: 0                                                         //
46 // Written by Piergiorgio Cerello                                     //
47 // November 23 1999                                                   //
48 //                                                                    //
49 // AliITSsimulationSDD is the simulation of SDDs.                     //
50 ////////////////////////////////////////////////////////////////////////
51
52 //______________________________________________________________________
53 AliITSsimulationSDD::AliITSsimulationSDD():
54 AliITSsimulation(),
55 fITS(0),
56 fHitMap2(0),
57 fHitSigMap2(0),
58 fHitNoiMap2(0),
59 fElectronics(0),
60 fInZR(0),
61 fInZI(0),
62 fOutZR(0),
63 fOutZI(0),
64 fAnodeFire(0),
65 fHis(0),
66 fFlag(kFALSE),
67 fCrosstalkFlag(kFALSE),
68 fDoFFT(1),
69 fNofMaps(0),
70 fMaxNofSamples(0),
71 fScaleSize(0){
72     // Default constructor
73     SetPerpendTracksFlag();
74     SetCrosstalkFlag();
75     SetDoFFT();
76 }
77 //______________________________________________________________________
78 AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
79     AliITSsimulation(source),
80 fITS(source.fITS),
81 fHitMap2(source.fHitMap2),
82 fHitSigMap2(source.fHitSigMap2),
83 fHitNoiMap2(source.fHitNoiMap2),
84 fElectronics(source.fElectronics),
85 fInZR(source.fInZR),
86 fInZI(source.fInZI),
87 fOutZR(source.fOutZR),
88 fOutZI(source.fOutZI),
89 fAnodeFire(source.fAnodeFire),
90 fHis(source.fHis),
91 fFlag(source.fFlag),
92 fCrosstalkFlag(source.fCrosstalkFlag),
93 fDoFFT(source.fDoFFT),
94 fNofMaps(source.fNofMaps),
95 fMaxNofSamples(source.fMaxNofSamples),
96 fScaleSize(source.fScaleSize){
97     // Copy constructor to satify Coding roules only.
98
99 }
100 //______________________________________________________________________
101 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
102     // Assignment operator to satify Coding roules only.
103
104     if(this==&src) return *this;
105     Error("AliITSsimulationSDD","Not allowed to make a = with "
106           "AliITSsimulationSDD Using default creater instead");
107     return *this ;
108 }
109 /*
110 //______________________________________________________________________
111 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
112     // Assignment operator to satify Coding roules only.
113
114     if(this==&src) return *this;
115     Error("AliITSsimulationSSD","Not allowed to make a = with "
116           "AliITSsimulationSDD Using default creater instead");
117     return *this ;
118 }
119 */
120 //______________________________________________________________________
121 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
122 AliITSsimulation(dettyp),
123 fITS(0),
124 fHitMap2(0),
125 fHitSigMap2(0),
126 fHitNoiMap2(0),
127 fElectronics(0),
128 fInZR(0),
129 fInZI(0),
130 fOutZR(0),
131 fOutZI(0),
132 fAnodeFire(0),
133 fHis(0),
134 fFlag(kFALSE),
135 fCrosstalkFlag(kFALSE),
136 fDoFFT(1),
137 fNofMaps(0),
138 fMaxNofSamples(0),
139 fScaleSize(0){
140     // Default Constructor
141   Init();
142 }
143 //______________________________________________________________________
144 void AliITSsimulationSDD::Init(){
145     // Standard Constructor
146
147     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
148     fScaleSize = ScaleFourier(seg);
149     SetPerpendTracksFlag();
150     SetCrosstalkFlag();
151     SetDoFFT();
152
153     AliITSSimuParam* simpar = fDetType->GetSimuParam();
154     fpList = new AliITSpList( seg->Npz(),
155                               fScaleSize*seg->Npx() );
156     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
157     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
158     fHitMap2 = fHitSigMap2;
159
160     fNofMaps = seg->Npz();
161     fMaxNofSamples = seg->Npx();
162     fAnodeFire = new Bool_t [fNofMaps];
163     
164     Float_t sddWidth  = seg->Dz();
165     Float_t anodePitch = seg->Dpz(0);
166     Double_t timeStep  = (Double_t)seg->Dpx(0);
167
168     if(anodePitch*(fNofMaps/2) > sddWidth) {
169       AliWarning(Form("Too many anodes %d or too big pitch %f ",
170                 fNofMaps/2,anodePitch));
171     } // end if
172
173
174     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
175                                     simpar->GetSDDElectronics());
176
177
178     fITS       = (AliITS*)gAlice->GetModule("ITS");
179  
180     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
181     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
182     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
183     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];  
184 }
185 //______________________________________________________________________
186 AliITSsimulationSDD::~AliITSsimulationSDD() { 
187     // destructor
188
189     //    delete fpList;
190     delete fHitSigMap2;
191     delete fHitNoiMap2;
192     delete fElectronics;
193
194     fITS = 0;
195
196     if (fHis) {
197         fHis->Delete(); 
198         delete fHis;     
199     } // end if fHis
200     if(fInZR)  delete [] fInZR;
201     if(fInZI)  delete [] fInZI;        
202     if(fOutZR) delete [] fOutZR;
203     if(fOutZI) delete [] fOutZI;
204     if(fAnodeFire) delete [] fAnodeFire;
205 }
206 //______________________________________________________________________
207 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
208     // create maps to build the lists of tracks for each summable digit
209     fModule = module;
210     fEvent  = event;
211     ClearMaps();
212     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);    
213 }
214 //______________________________________________________________________
215 void AliITSsimulationSDD::ClearMaps() {
216     // clear maps
217     fpList->ClearMap();
218     fHitSigMap2->ClearMap();
219     fHitNoiMap2->ClearMap();
220 }
221 //______________________________________________________________________
222 void AliITSsimulationSDD::FastFourierTransform(Double_t *real,
223                           Double_t *imag,Int_t direction) {
224     // Do a Fast Fourier Transform
225
226     Int_t samples = fElectronics->GetSamples();
227     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
228     Int_t m1 = samples;
229     Int_t m  = samples/2;
230     Int_t m2 = samples/m1;
231     Int_t i,j,k;
232     for(i=1; i<=l; i++) {
233         for(j=0; j<samples; j += m1) {
234             Int_t p = 0;
235             for(k=j; k<= j+m-1; k++) {
236                 Double_t wsr = fElectronics->GetWeightReal(p);
237                 Double_t wsi = fElectronics->GetWeightImag(p);
238                 if(direction == -1) wsi = -wsi;
239                 Double_t xr = *(real+k+m);
240                 Double_t xi = *(imag+k+m);
241                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
242                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
243                 *(real+k) += xr;
244                 *(imag+k) += xi;
245                 p += m2;
246             } // end for k
247         } // end for j
248         m1 = m;
249         m /= 2;
250         m2 += m2;
251     } // end for i
252     for(j=0; j<samples; j++) {
253         Int_t j1 = j;
254         Int_t p = 0;
255         Int_t i1;
256         for(i1=1; i1<=l; i1++) {
257             Int_t j2 = j1;
258             j1 /= 2;
259             p = p + p + j2 - j1 - j1;
260         } // end for i1
261         if(p >= j) {
262             Double_t xr = *(real+j);
263             Double_t xi = *(imag+j);
264             *(real+j) = *(real+p);
265             *(imag+j) = *(imag+p);
266             *(real+p) = xr;
267             *(imag+p) = xi;
268         } // end if p>=j
269     } // end for j
270     if(direction == -1) {
271         for(i=0; i<samples; i++) {
272             *(real+i) /= samples;
273             *(imag+i) /= samples;
274         } // end for i
275     } // end if direction == -1
276     return;
277 }
278
279 //______________________________________________________________________
280 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
281     // digitize module using the "slow" detector simulator creating
282     // summable digits.
283
284     TObjArray *fHits = mod->GetHits();
285     Int_t nhits      = fHits->GetEntriesFast();
286     if( !nhits ) return;
287
288     InitSimulationModule( md, ev );
289     HitsToAnalogDigits( mod );  // fills fHitMap2 which is = fHitSigmap2
290     ChargeToSignal( fModule,kFALSE,kTRUE ); // - Process signal adding gain without adding noise
291     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
292     ChargeToSignal( fModule,kTRUE,kFALSE );  // - Process only noise
293     fHitMap2 = fHitSigMap2;   // - Return to signal map
294     WriteSDigits();
295     ClearMaps();
296 }
297 //______________________________________________________________________
298 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
299                                                Int_t mask ) {
300     // Add Summable digits to module maps.
301     AliITSSimuParam* simpar = fDetType->GetSimuParam();
302     Int_t    nItems = pItemArray->GetEntries();
303     Double_t maxadc = simpar->GetSDDMaxAdc();
304     Bool_t sig = kFALSE;
305     
306     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
307     for( Int_t i=0; i<nItems; i++ ) {
308         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
309         if( pItem->GetModule() != fModule ) {
310             Error( "AliITSsimulationSDD","Error reading, SDigits module "
311                    "%d != current module %d: exit",
312                    pItem->GetModule(), fModule );
313             return sig;
314         } // end if
315
316         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
317         
318         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
319         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
320         Double_t sigAE = pItem2->GetSignalAfterElect();
321         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
322         Int_t ia;
323         Int_t it;
324         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
325         fHitMap2->SetHit( ia, it, sigAE );
326         fAnodeFire[ia] = kTRUE;
327     }
328     return sig;
329 }
330 //______________________________________________________________________
331 void AliITSsimulationSDD::FinishSDigitiseModule() {
332     // digitize module using the "slow" detector simulator from
333     // the sum of summable digits.
334     FinishDigits() ;
335     ClearMaps();
336 }
337 //______________________________________________________________________
338 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
339     // create maps to build the lists of tracks for each digit
340
341     TObjArray *fHits = mod->GetHits();
342     Int_t nhits      = fHits->GetEntriesFast();
343
344     InitSimulationModule( md, ev );
345     if( !nhits ) return;
346         
347     HitsToAnalogDigits( mod );
348     ChargeToSignal( fModule,kTRUE,kTRUE );  // process signal + noise
349
350     for( Int_t i=0; i<fNofMaps; i++ ) {
351         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
352             Int_t jdx = j*fScaleSize;
353             Int_t index = fpList->GetHitIndex( i, j );
354             AliITSpListItem pItemTmp2( fModule, index, 0. );
355             // put the fScaleSize analog digits in only one
356             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
357                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
358                 if( pItemTmp == 0 ) continue;
359                 pItemTmp2.Add( pItemTmp );
360             }
361             fpList->DeleteHit( i, j );
362             fpList->AddItemTo( 0, &pItemTmp2 );
363         }
364     }
365     FinishDigits();
366     ClearMaps();
367 }
368 //______________________________________________________________________
369 void AliITSsimulationSDD::FinishDigits() {
370     // introduce the electronics effects and do zero-suppression if required
371
372     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
373
374     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
375     Bool_t isZeroSupp = res->GetZeroSupp();
376     if (isZeroSupp) Compress2D();
377     else StoreAllDigits();
378 }
379 //______________________________________________________________________
380 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
381     // create maps to build the lists of tracks for each digit
382   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
383   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
384   AliITSSimuParam* simpar = fDetType->GetSimuParam();
385   TObjArray *hits     = mod->GetHits();
386   Int_t      nhits    = hits->GetEntriesFast();
387
388   //    Int_t      arg[6]   = {0,0,0,0,0,0};
389   Int_t     nofAnodes  = fNofMaps/2;
390   Double_t  sddLength  = seg->Dx();
391   Double_t  anodePitch = seg->Dpz(0);
392   Double_t  timeStep   = seg->Dpx(0);
393   Double_t  driftSpeed ;  // drift velocity (anode dependent)
394   Double_t  nanoampToADC       = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); //   maxadc/topValue;
395   Double_t  cHloss     = simpar->GetSDDChargeLoss();
396   Float_t   dfCoeff, s1; 
397   simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape
398   Double_t  eVpairs    = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
399   Double_t  nsigma     = simpar->GetNSigmaIntegration(); //
400   Int_t     nlookups   = simpar->GetGausNLookUp();       //
401   Float_t   jitter     = simpar->GetSDDJitterError(); // 
402   Float_t   mapsmear   = simpar->GetSDDCorrMapPrecision(); // 
403   Float_t   trigDelay  = simpar->GetSDDTrigDelay(); // compensation for MC time zero
404   if(res->IsAMAt20MHz()) trigDelay+=12.5; // compensation for discretization step
405
406   Float_t   timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule);
407   Float_t   adcscale   = fDetType->GetResponseSDD()->GetADCtokeV(fModule);
408   adcscale/=simpar->GetSDDkeVtoADC();
409
410   // Piergiorgio's part (apart for few variables which I made float
411   // when i thought that can be done
412   // Fill detector maps with GEANT hits
413   // loop over hits in the module
414   
415   const Float_t kconv = 1.0e+6;  // GeV->KeV
416   Int_t     itrack      = 0;
417   Int_t     iWing;       // which detector wing/side.
418   Int_t     ii,kk,ka,kt; // loop indexs
419   Int_t     ia,it,index; // sub-pixel integration indexies
420   Int_t     iAnode;      // anode number.
421   Int_t     timeSample;  // time buckett.
422   Int_t     anodeWindow; // anode direction charge integration width
423   Int_t     timeWindow;  // time direction charge integration width
424   Int_t     jamin,jamax; // anode charge integration window
425   Int_t     jtmin,jtmax; // time charge integration window
426   Int_t     nsplitAn;    // the number of splits in anode and time windows
427   Int_t     nsplitTb;    // the number of splits in anode and time windows
428   Int_t     nOfSplits;   // number of times track length is split into
429   Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
430   Float_t   kkF;         // Floating point version of loop index kk.
431   Double_t  pathInSDD; // Track length in SDD.
432   Double_t  drPath; // average position of track in detector. in microns
433   Double_t  drTime; // Drift time
434   Double_t  avDrft;  // x position of path length segment in cm.
435   Double_t  avAnode; // Anode for path length segment in Anode number (float)
436   Double_t  zAnode;  // Floating point anode number.
437   Double_t  driftPath; // avDrft in microns.
438   Double_t  width;     // width of signal at anodes.
439   Double_t  depEnergy; // Energy deposited in this GEANT step.
440   Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
441   Double_t  sigA; // sigma of signal at anode.
442   Double_t  sigT; // sigma in time/drift direction for track segment
443   Double_t  aStep,aConst; // sub-pixel size and offset anode
444   Double_t  tStep,tConst; // sub-pixel size and offset time
445   Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
446   Double_t  chargeloss; // charge loss for track segment.
447   Double_t  anodeAmplitude; // signal amplitude in anode direction
448   Double_t  aExpo;          // exponent of Gaussian anode direction
449   Double_t  timeAmplitude;  // signal amplitude in time direction
450   Double_t  tExpo;          // exponent of Gaussian time direction
451   Double_t  tof;            // Time of flight in ns of this step.    
452   
453   for(ii=0; ii<nhits; ii++) {
454     if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
455                           depEnergy,itrack)) continue;
456     Float_t xloc=xL[0];
457     if(xloc>0) iWing=0; // left side, carlos channel 0
458     else iWing=1; // right side
459     
460     Float_t zloc=xL[2]+0.5*dxL[2];
461     zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
462     driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
463     driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
464
465     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
466       AliWarning("Time Interval > Allowed Time Interval");
467     }
468     depEnergy  *= kconv;
469     if (!depEnergy) {
470       AliDebug(1,
471                Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
472                     itrack,ii,mod->GetIndex()));
473       continue;
474       // continue if the particle did not lose energy
475       // passing through detector
476     } // end if !depEnergy
477      
478     tof=0.;
479     AliITShit* h=(AliITShit*)hits->At(ii);
480     if(h){ 
481       tof=h->GetTOF()*1E9; 
482       AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof));
483     }
484
485     Float_t corrx=0, corrz=0;
486     res->GetShiftsForSimulation(xL[2],xL[0],corrz,corrx,seg);
487     xL[2]-=corrz;
488     xL[0]-=corrx;
489     xL[0] += 0.0001*gRandom->Gaus( 0, mapsmear); //
490     xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
491
492     pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
493     
494     if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
495     drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
496     drPath = sddLength-drPath;
497     if(drPath < 0) {
498       AliInfo( // this should be fixed at geometry level
499                Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
500                     drPath,sddLength,dxL[0],xL[0]));
501       continue;
502     } // end if drPath < 0
503     
504     // Compute number of segments to brake step path into
505     drTime = drPath/driftSpeed;  //   Drift Time
506     sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
507     // calcuate the number of time the path length should be split into.
508     nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
509     if(fFlag) nOfSplits = 1;
510     
511     // loop over path segments, init. some variables.
512     depEnergy /= nOfSplits;
513     nOfSplitsF = (Float_t) nOfSplits;
514     Float_t theAverage=0.,theSteps=0.;
515     for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
516       kkF       = (Float_t) kk + 0.5;
517       avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
518       avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
519       theSteps+=1.;
520       theAverage+=avAnode;
521       zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
522       driftSpeed = res->GetDriftSpeedAtAnode(zAnode);   
523       driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
524       driftPath = TMath::Abs(10000.*avDrft);
525       driftPath = sddLength-driftPath;
526       if(driftPath < 0) {
527         AliDebug(1, // this should be fixed at geometry level
528                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
529                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
530         continue;
531       } // end if driftPath < 0
532       drTime     = driftPath/driftSpeed; // drift time for segment.
533       // Sigma along the anodes for track segment.
534       sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
535       sigT       = sigA/driftSpeed;
536
537       drTime+=tof; // take into account Time Of Flight from production point
538       drTime-=trigDelay;
539       drTime+=timeZero;
540       timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
541       if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
542       iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
543       
544         // Peak amplitude in nanoAmpere
545       amplitude  = fScaleSize*160.*depEnergy/
546         (timeStep*eVpairs*2.*acos(-1.));
547       chargeloss = 1.-cHloss*driftPath/1000.;
548       amplitude *= chargeloss;
549       amplitude *= adcscale;
550       width  = 2.*nsigma/(nlookups-1);
551       // Spread the charge 
552       nsplitAn = 4; 
553       nsplitTb=4;
554       aStep  = anodePitch/(nsplitAn*sigA);
555       aConst = zAnode*anodePitch/sigA;
556       tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
557       tConst = drTime/sigT;
558       // Define SDD window corresponding to the hit
559       anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
560       timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
561       jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
562       if(jamin <= 0) jamin = 1;
563       if(jamin > nofAnodes*nsplitAn){ 
564         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
565         continue;
566       }
567       jamax = (iAnode + anodeWindow + 2)*nsplitAn;
568       if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn;
569       if(jamax <=0){ 
570         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax));
571         continue;
572       }
573       jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
574       if(jtmin <= 0) jtmin = 1;
575       if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){ 
576         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d  tof=%f",jtmin,tof));
577         continue; 
578       }
579       jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
580       if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
581       if(jtmax <= 0){
582         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d  tof=%f",jtmax,tof));
583         continue; 
584       }
585
586       // Spread the charge in the anode-time window
587       for(ka=jamin; ka <=jamax; ka++) {   
588         ia = (ka-1)/nsplitAn + 1;
589         if(ia <= 0) ia=1; 
590         if(ia > nofAnodes) ia = nofAnodes;
591         aExpo     = (aStep*(ka-0.5)-aConst);
592         if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
593         else {
594           Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
595           anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
596         }
597         // index starts from 0
598         index = iWing*nofAnodes+ia-1;
599         if(anodeAmplitude){
600           for(kt=jtmin; kt<=jtmax; kt++) {
601             it = (kt-1)/nsplitTb+1;  // it starts from 1
602             if(it<=0) it=1;
603             if(it>fScaleSize*fMaxNofSamples)
604               it = fScaleSize*fMaxNofSamples;
605             tExpo    = (tStep*(kt-0.5)-tConst);
606             if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
607             else {
608               Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
609               timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
610             }
611             timeAmplitude *= nanoampToADC;
612             //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
613             Double_t charge = timeAmplitude;
614             charge += fHitMap2->GetSignal(index,it-1);
615             fHitMap2->SetHit(index, it-1, charge);
616             fpList->AddSignal(index,it-1,itrack,ii-1,
617                               mod->GetIndex(),timeAmplitude);
618             fAnodeFire[index] = kTRUE;
619           }  // end loop over time in window               
620         } // end if anodeAmplitude 
621       } // loop over anodes in window
622     } // end loop over "sub-hits"
623   } // end loop over hits
624 }
625
626 //____________________________________________
627 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
628   // Adds a Digit.
629   Int_t size = AliITSdigit::GetNTracks();
630
631   Int_t digits[3];
632   Int_t * tracks = new Int_t[size];
633   Int_t * hits = new Int_t[size];
634   Float_t phys;
635   Float_t * charges = new Float_t[size];
636
637   digits[0] = i;
638   digits[1] = j;
639   digits[2] = signalc;
640
641   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
642   if( pItem == 0 ) {
643     phys = 0.0;
644     for( Int_t l=0; l<size; l++ ) {
645       tracks[l]  = 0;
646       hits[l]    = 0;
647       charges[l] = 0.0;
648     }
649   } else {
650     Int_t idtrack =  pItem->GetTrack( 0 );
651     if( idtrack >= 0 ) phys = pItem->GetSignal();  
652     else phys = 0.0;
653
654     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
655       tracks[l]  = pItem->GetTrack( l );
656       hits[l]    = pItem->GetHit( l );
657       charges[l] = pItem->GetSignal( l );
658     }else{
659       tracks[l]  = -3;
660       hits[l]    = -1;
661       charges[l] = 0.0;
662     }// end for if
663   }
664
665   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale ); 
666   delete [] tracks;
667   delete [] hits;
668   delete [] charges;
669 }
670 //______________________________________________________________________
671 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
672   // add baseline, noise, gain, electronics and ADC saturation effects
673   // apply dead channels
674
675   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
676   Double_t baseline=0; 
677   Double_t noise=0; 
678   Double_t gain=0; 
679   Float_t contrib=0;
680   Int_t i,k,kk;
681   AliITSSimuParam* simpar = fDetType->GetSimuParam();
682   Float_t maxadc = simpar->GetSDDMaxAdc();    
683   Int_t nGroup=fScaleSize;
684   if(res->IsAMAt20MHz()){
685     nGroup=fScaleSize/2;
686   }
687
688   for (i=0;i<fNofMaps;i++) {
689     if( !fAnodeFire[i] ) continue;
690     baseline = res->GetBaseline(i);
691     noise = res->GetNoise(i);
692     gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD();
693     if(res->IsBad()) gain=0.;
694     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
695     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
696       fInZR[k]  = fHitMap2->GetSignal(i,k);
697       if(bAddGain) fInZR[k]*=gain;
698       if( bAddNoise ) {
699         contrib   = (baseline + noise*gRandom->Gaus());
700         fInZR[k] += contrib;
701       }
702       fInZI[k]  = 0.;
703     } // end for k
704     if(!fDoFFT) {      
705       for(k=0; k<fMaxNofSamples; k++) {
706         Double_t newcont = 0.;
707         Double_t maxcont = 0.;
708         for(kk=0;kk<fScaleSize;kk++) {
709           newcont = fInZR[fScaleSize*k+kk];
710           if(newcont > maxcont) maxcont = newcont;
711         } // end for kk
712         newcont = maxcont;
713         if (newcont >= maxadc) newcont = maxadc -1;
714         if(newcont >= baseline){
715           Warning("","newcont=%f>=baseline=%f",newcont,baseline);
716         } // end if
717           // back to analog: ?
718         fHitMap2->SetHit(i,k,newcont);
719       }  // end for k
720     }else{
721       FastFourierTransform(&fInZR[0],&fInZI[0],1);
722       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
723         Double_t rw = fElectronics->GetTraFunReal(k);
724         Double_t iw = fElectronics->GetTraFunImag(k);
725         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
726         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
727       } // end for k
728       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
729       for(k=0; k<fMaxNofSamples; k++) {
730         Double_t newcont1 = 0.;
731         Double_t maxcont1 = 0.;
732         for(kk=0;kk<nGroup;kk++) {
733           newcont1 = fOutZR[fScaleSize*k+kk];
734           if(newcont1 > maxcont1) maxcont1 = newcont1;
735         } // end for kk
736         newcont1 = maxcont1;
737         if (newcont1 >= maxadc) newcont1 = maxadc -1;
738         fHitMap2->SetHit(i,k,newcont1);
739       } // end for k
740     }
741   } // end for i loop over anodes
742   return;
743 }
744
745 //______________________________________________________________________
746 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
747     // function add the crosstalk effect to signal
748     // temporal function, should be checked...!!!
749   
750     // create and inizialice crosstalk map
751     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
752     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
753     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
754     for( Int_t z=0; z<fNofMaps; z++ ) {
755       Double_t baseline = calibr->GetBaseline(z);
756         Bool_t on = kFALSE;
757         Int_t tstart = 0;
758         Int_t tstop = 0;
759         Int_t nTsteps = 0;
760         
761         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
762             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
763             if( fadc > baseline ) {
764                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
765                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
766                     if( fadc1 < fadc ) continue;
767                     on = kTRUE;
768                     nTsteps = 0;
769                     tstart = l;
770                 }
771                 nTsteps++;
772             }
773             else { // end fadc > baseline
774                 if( on == kTRUE ) {        
775                     if( nTsteps > 2 ) {
776                         tstop = l;
777                         // make smooth derivative
778                         Float_t* dev = new Float_t[fMaxNofSamples+1];
779                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
780                         for( Int_t i=tstart; i<tstop; i++ ) {   
781                             if( i > 2 && i < fMaxNofSamples-2 )
782                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 ) 
783                                     -0.1*fHitMap2->GetSignal( z,i-1 ) 
784                                     +0.1*fHitMap2->GetSignal( z,i+1 ) 
785                                     +0.2*fHitMap2->GetSignal( z,i+2 );
786                         }
787                         
788                         // add crosstalk contribution to neibourg anodes  
789                         for( Int_t i=tstart; i<tstop; i++ ) {
790                             Int_t anode = z - 1;
791                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); // 
792                             Float_t ctktmp =  -dev[i1] * 0.25;
793                             if( anode > 0 ) {
794                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
795                             }
796                             anode = z + 1;
797                             if( anode < fNofMaps ) {
798                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
799                             }
800                         }
801                         delete [] dev;
802                         
803                     } // if( nTsteps > 2 )
804                     on = kFALSE;
805                 }  // if( on == kTRUE )
806             }  // else
807         }
808     }
809     
810     for( Int_t a=0; a<fNofMaps; a++ )
811         for( Int_t t=0; t<fMaxNofSamples; t++ ) {     
812             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
813             fHitMap2->SetHit( a, t, signal );
814         }
815
816     delete [] ctk;
817 }
818
819 //______________________________________________________________________
820 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
821     // To the 10 to 8 bit lossive compression.
822     // code from Davide C. and Albert W.
823
824     if (signal < 128)  return signal;
825     if (signal < 256)  return (128+((signal-128)>>1));
826     if (signal < 512)  return (192+((signal-256)>>3));
827     if (signal < 1024) return (224+((signal-512)>>4));
828     return 0;
829 }
830 //______________________________________________________________________
831 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
832   // Decompression from 8 to 10 bit
833
834   if (signal < 0 || signal > 255) {
835     AliWarning(Form("Signal value %d out of range",signal));
836     return 0;
837   } // end if signal <0 || signal >255
838
839   if (signal < 128) return signal;
840   if (signal < 192) {
841     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
842     else  return (128+((signal-128)<<1)+1);
843   } // end if signal < 192
844   if (signal < 224) {
845     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
846     else  return (256+((signal-192)<<3)+4);
847   } // end if signal < 224
848   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
849   return (512+((signal-224)<<4)+8);
850 }
851 //______________________________________________________________________
852 void AliITSsimulationSDD::Compress2D(){
853   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
854   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);  
855   for (Int_t iWing=0; iWing<2; iWing++) {
856     Int_t tL=res->GetZSLowThreshold(iWing);
857     Int_t tH=res->GetZSHighThreshold(iWing);
858     for (Int_t i=0; i<fNofMaps/2; i++) {  
859       Int_t ian=i+iWing*fNofMaps/2;
860       if( !fAnodeFire[ian] ) continue;
861       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
862         Int_t nLow=0, nHigh=0;      
863         Float_t cC=fHitMap2->GetSignal(ian,itb);
864         if(cC<=tL) continue;
865         nLow++; // cC is greater than tL
866         if(cC>tH) nHigh++;
867         //                     N
868         // Get "quintuple":   WCE
869         //                     S
870         Float_t wW=0.;
871         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
872         if(wW>tL) nLow++;
873         if(wW>tH) nHigh++;
874         Float_t eE=0.;
875         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
876         if(eE>tL) nLow++;
877         if(eE>tH) nHigh++;
878         Float_t nN=0.;
879         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
880         if(nN>tL) nLow++;
881         if(nN>tH) nHigh++;
882         Float_t sS=0.;
883         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
884         if(sS>tL) nLow++;
885         if(sS>tH) nHigh++;
886         
887         if(nLow>=2 && nHigh>=1){
888           Int_t signal=(Int_t)cC;
889           Int_t signalc = Convert10to8(signal);
890           Int_t signale = Convert8to10(signalc);
891           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
892           if(signalc>=4) AddDigit(ian,itb,signalc,signale);  // store C 
893         }
894       }
895     }
896   }
897 }
898
899
900 //______________________________________________________________________
901 void AliITSsimulationSDD::StoreAllDigits(){
902   // store digits for non-zero-suppressed data
903   for (Int_t ian=0; ian<fNofMaps; ian++) {
904     for (Int_t itb=0; itb<fMaxNofSamples; itb++){
905       Int_t signal=(Int_t)(fHitMap2->GetSignal(ian,itb));
906       Int_t signalc = Convert10to8(signal);
907       Int_t signale = Convert8to10(signalc);
908       AddDigit(ian,itb,signalc,signale);  
909     } 
910   }
911
912 //______________________________________________________________________
913 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
914   // Creates histograms of maps for debugging
915   Int_t i;
916   
917   fHis=new TObjArray(fNofMaps);
918   for (i=0;i<fNofMaps;i++) {
919     TString sddName;
920     sddName.Form("sdd_%d",i+1);
921     fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
922                          0.,(Float_t) scale*fMaxNofSamples), i);
923   } // end for i
924 }
925 //______________________________________________________________________
926 void AliITSsimulationSDD::FillHistograms(){
927     // fill 1D histograms from map
928
929     if (!fHis) return;
930
931     for( Int_t i=0; i<fNofMaps; i++) {
932         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
933         Int_t nsamples = hist->GetNbinsX();
934         for( Int_t j=0; j<nsamples; j++) {
935             Double_t signal=fHitMap2->GetSignal(i,j);
936             hist->Fill((Float_t)j,signal);
937         } // end for j
938     } // end for i
939 }
940 //______________________________________________________________________
941 void AliITSsimulationSDD::ResetHistograms(){
942     // Reset histograms for this detector
943     Int_t i;
944
945     for (i=0;i<fNofMaps;i++ ) {
946         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
947     } // end for i
948 }
949 //______________________________________________________________________
950 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) { 
951     // Fills a histogram from a give anode.  
952
953     if (!fHis) return 0;
954
955     if(wing <=0 || wing > 2) {
956         Warning("GetAnode","Wrong wing number: %d",wing);
957         return NULL;
958     } // end if wing <=0 || wing >2
959     if(anode <=0 || anode > fNofMaps/2) {
960         Warning("GetAnode","Wrong anode number: %d",anode);
961         return NULL;
962     } // end if ampde <=0 || andoe > fNofMaps/2
963
964     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
965     return (TH1F*)(fHis->At(index));
966 }
967 //______________________________________________________________________
968 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
969     // Writes the histograms to a file
970
971     if (!fHis) return;
972
973     hfile->cd();
974     Int_t i;
975     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
976     return;
977 }
978 //______________________________________________________________________
979 void AliITSsimulationSDD::WriteSDigits(){
980     // Fills the Summable digits Tree
981     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
982
983     for( Int_t i=0; i<fNofMaps; i++ ) {
984         if( !fAnodeFire[i] ) continue;
985         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
986             Double_t sig = fHitMap2->GetSignal( i, j );
987             if( sig > 0.2 ) {
988                 Int_t jdx = j*fScaleSize;
989                 Int_t index = fpList->GetHitIndex( i, j );
990                 AliITSpListItem pItemTmp2( fModule, index, 0. );
991                 // put the fScaleSize analog digits in only one
992                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
993                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
994                     if( pItemTmp == 0 ) continue;
995                     pItemTmp2.Add( pItemTmp );
996                 }
997                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
998                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
999                 aliITS->AddSumDigit( pItemTmp2 );
1000             } // end if (sig > 0.2)
1001         }
1002     }
1003     return;
1004 }
1005 //______________________________________________________________________
1006 void AliITSsimulationSDD::PrintStatus() const {
1007     // Print SDD simulation Parameters
1008
1009     cout << "**************************************************" << endl;
1010     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1011     cout << "**************************************************" << endl;
1012     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1013     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1014     cout << "Number of Anodes used: " << fNofMaps << endl;
1015     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1016     cout << "Scale size factor: " << fScaleSize << endl;
1017     cout << "**************************************************" << endl;
1018 }