Fix coverity warnings (F. Prino)
[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   trigDelay  = simpar->GetSDDTrigDelay(); // compensation for MC time zero
403   Float_t   timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule);
404
405   // Piergiorgio's part (apart for few variables which I made float
406   // when i thought that can be done
407   // Fill detector maps with GEANT hits
408   // loop over hits in the module
409   
410   const Float_t kconv = 1.0e+6;  // GeV->KeV
411   Int_t     itrack      = 0;
412   Int_t     iWing;       // which detector wing/side.
413   Int_t     ii,kk,ka,kt; // loop indexs
414   Int_t     ia,it,index; // sub-pixel integration indexies
415   Int_t     iAnode;      // anode number.
416   Int_t     timeSample;  // time buckett.
417   Int_t     anodeWindow; // anode direction charge integration width
418   Int_t     timeWindow;  // time direction charge integration width
419   Int_t     jamin,jamax; // anode charge integration window
420   Int_t     jtmin,jtmax; // time charge integration window
421   Int_t     nsplitAn;    // the number of splits in anode and time windows
422   Int_t     nsplitTb;    // the number of splits in anode and time windows
423   Int_t     nOfSplits;   // number of times track length is split into
424   Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
425   Float_t   kkF;         // Floating point version of loop index kk.
426   Double_t  pathInSDD; // Track length in SDD.
427   Double_t  drPath; // average position of track in detector. in microns
428   Double_t  drTime; // Drift time
429   Double_t  avDrft;  // x position of path length segment in cm.
430   Double_t  avAnode; // Anode for path length segment in Anode number (float)
431   Double_t  zAnode;  // Floating point anode number.
432   Double_t  driftPath; // avDrft in microns.
433   Double_t  width;     // width of signal at anodes.
434   Double_t  depEnergy; // Energy deposited in this GEANT step.
435   Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
436   Double_t  sigA; // sigma of signal at anode.
437   Double_t  sigT; // sigma in time/drift direction for track segment
438   Double_t  aStep,aConst; // sub-pixel size and offset anode
439   Double_t  tStep,tConst; // sub-pixel size and offset time
440   Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
441   Double_t  chargeloss; // charge loss for track segment.
442   Double_t  anodeAmplitude; // signal amplitude in anode direction
443   Double_t  aExpo;          // exponent of Gaussian anode direction
444   Double_t  timeAmplitude;  // signal amplitude in time direction
445   Double_t  tExpo;          // exponent of Gaussian time direction
446   Double_t  tof;            // Time of flight in ns of this step.    
447   
448   for(ii=0; ii<nhits; ii++) {
449     if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
450                           depEnergy,itrack)) continue;
451     Float_t xloc=xL[0];
452     if(xloc>0) iWing=0; // left side, carlos channel 0
453     else iWing=1; // right side
454     
455     Float_t zloc=xL[2]+0.5*dxL[2];
456     zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
457     driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
458     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
459       AliWarning("Time Interval > Allowed Time Interval");
460     }
461     depEnergy  *= kconv;
462     if (!depEnergy) {
463       AliDebug(1,
464                Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
465                     itrack,ii,mod->GetIndex()));
466       continue;
467       // continue if the particle did not lose energy
468       // passing through detector
469     } // end if !depEnergy
470      
471     tof=0.;
472     AliITShit* h=(AliITShit*)hits->At(ii);
473     if(h){ 
474       tof=h->GetTOF()*1E9; 
475       AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof));
476     }
477
478     xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
479     pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
480     
481     if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
482     drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
483     drPath = sddLength-drPath;
484     if(drPath < 0) {
485       AliDebug(1, // this should be fixed at geometry level
486                Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
487                     drPath,sddLength,dxL[0],xL[0]));
488       continue;
489     } // end if drPath < 0
490     
491     // Compute number of segments to brake step path into
492     drTime = drPath/driftSpeed;  //   Drift Time
493     sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
494     // calcuate the number of time the path length should be split into.
495     nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
496     if(fFlag) nOfSplits = 1;
497     
498     // loop over path segments, init. some variables.
499     depEnergy /= nOfSplits;
500     nOfSplitsF = (Float_t) nOfSplits;
501     Float_t theAverage=0.,theSteps=0.;
502     for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
503       kkF       = (Float_t) kk + 0.5;
504       avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
505       avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
506       theSteps+=1.;
507       theAverage+=avAnode;
508       zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
509       driftSpeed = res->GetDriftSpeedAtAnode(zAnode);   
510       driftPath = TMath::Abs(10000.*avDrft);
511       driftPath = sddLength-driftPath;
512       if(driftPath < 0) {
513         AliDebug(1, // this should be fixed at geometry level
514                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
515                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
516         continue;
517       } // end if driftPath < 0
518       drTime     = driftPath/driftSpeed; // drift time for segment.
519       // Sigma along the anodes for track segment.
520       sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
521       sigT       = sigA/driftSpeed;
522
523       drTime+=tof; // take into account Time Of Flight from production point
524       drTime-=trigDelay;
525       drTime+=timeZero;
526       timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
527       if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
528       iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
529       
530         // Peak amplitude in nanoAmpere
531       amplitude  = fScaleSize*160.*depEnergy/
532         (timeStep*eVpairs*2.*acos(-1.));
533       chargeloss = 1.-cHloss*driftPath/1000.;
534       amplitude *= chargeloss;
535       width  = 2.*nsigma/(nlookups-1);
536       // Spread the charge 
537       nsplitAn = 4; 
538       nsplitTb=4;
539       aStep  = anodePitch/(nsplitAn*sigA);
540       aConst = zAnode*anodePitch/sigA;
541       tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
542       tConst = drTime/sigT;
543       // Define SDD window corresponding to the hit
544       anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
545       timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
546       jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
547       if(jamin <= 0) jamin = 1;
548       if(jamin > nofAnodes*nsplitAn){ 
549         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
550         continue;
551       }
552       jamax = (iAnode + anodeWindow + 2)*nsplitAn;
553       if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn;
554       if(jamax <=0){ 
555         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax));
556         continue;
557       }
558       jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
559       if(jtmin <= 0) jtmin = 1;
560       if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){ 
561         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d  tof=%f",jtmin,tof));
562         continue; 
563       }
564       jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
565       if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
566       if(jtmax <= 0){
567         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d  tof=%f",jtmax,tof));
568         continue; 
569       }
570
571       // Spread the charge in the anode-time window
572       for(ka=jamin; ka <=jamax; ka++) {   
573         ia = (ka-1)/nsplitAn + 1;
574         if(ia <= 0) ia=1; 
575         if(ia > nofAnodes) ia = nofAnodes;
576         aExpo     = (aStep*(ka-0.5)-aConst);
577         if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
578         else {
579           Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
580           anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
581         }
582         // index starts from 0
583         index = iWing*nofAnodes+ia-1;
584         if(anodeAmplitude){
585           for(kt=jtmin; kt<=jtmax; kt++) {
586             it = (kt-1)/nsplitTb+1;  // it starts from 1
587             if(it<=0) it=1;
588             if(it>fScaleSize*fMaxNofSamples)
589               it = fScaleSize*fMaxNofSamples;
590             tExpo    = (tStep*(kt-0.5)-tConst);
591             if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
592             else {
593               Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
594               timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
595             }
596             timeAmplitude *= nanoampToADC;
597             //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
598             Double_t charge = timeAmplitude;
599             charge += fHitMap2->GetSignal(index,it-1);
600             fHitMap2->SetHit(index, it-1, charge);
601             fpList->AddSignal(index,it-1,itrack,ii-1,
602                               mod->GetIndex(),timeAmplitude);
603             fAnodeFire[index] = kTRUE;
604           }  // end loop over time in window               
605         } // end if anodeAmplitude 
606       } // loop over anodes in window
607     } // end loop over "sub-hits"
608   } // end loop over hits
609 }
610
611 //____________________________________________
612 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
613   // Adds a Digit.
614   Int_t size = AliITSdigit::GetNTracks();
615
616   Int_t digits[3];
617   Int_t * tracks = new Int_t[size];
618   Int_t * hits = new Int_t[size];
619   Float_t phys;
620   Float_t * charges = new Float_t[size];
621
622   digits[0] = i;
623   digits[1] = j;
624   digits[2] = signalc;
625
626   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
627   if( pItem == 0 ) {
628     phys = 0.0;
629     for( Int_t l=0; l<size; l++ ) {
630       tracks[l]  = 0;
631       hits[l]    = 0;
632       charges[l] = 0.0;
633     }
634   } else {
635     Int_t idtrack =  pItem->GetTrack( 0 );
636     if( idtrack >= 0 ) phys = pItem->GetSignal();  
637     else phys = 0.0;
638
639     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
640       tracks[l]  = pItem->GetTrack( l );
641       hits[l]    = pItem->GetHit( l );
642       charges[l] = pItem->GetSignal( l );
643     }else{
644       tracks[l]  = -3;
645       hits[l]    = -1;
646       charges[l] = 0.0;
647     }// end for if
648   }
649
650   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale ); 
651   delete [] tracks;
652   delete [] hits;
653   delete [] charges;
654 }
655 //______________________________________________________________________
656 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
657   // add baseline, noise, gain, electronics and ADC saturation effects
658   // apply dead channels
659
660   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
661   Double_t baseline=0; 
662   Double_t noise=0; 
663   Double_t gain=0; 
664   Float_t contrib=0;
665   Int_t i,k,kk;
666   AliITSSimuParam* simpar = fDetType->GetSimuParam();
667   Float_t maxadc = simpar->GetSDDMaxAdc();    
668   Int_t nGroup=fScaleSize;
669   if(res->IsAMAt20MHz()){
670     nGroup=fScaleSize/2;
671   }
672
673   for (i=0;i<fNofMaps;i++) {
674     if( !fAnodeFire[i] ) continue;
675     baseline = res->GetBaseline(i);
676     noise = res->GetNoise(i);
677     gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD();
678     if(res->IsBad()) gain=0.;
679     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
680     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
681       fInZR[k]  = fHitMap2->GetSignal(i,k);
682       if(bAddGain) fInZR[k]*=gain;
683       if( bAddNoise ) {
684         contrib   = (baseline + noise*gRandom->Gaus());
685         fInZR[k] += contrib;
686       }
687       fInZI[k]  = 0.;
688     } // end for k
689     if(!fDoFFT) {      
690       for(k=0; k<fMaxNofSamples; k++) {
691         Double_t newcont = 0.;
692         Double_t maxcont = 0.;
693         for(kk=0;kk<fScaleSize;kk++) {
694           newcont = fInZR[fScaleSize*k+kk];
695           if(newcont > maxcont) maxcont = newcont;
696         } // end for kk
697         newcont = maxcont;
698         if (newcont >= maxadc) newcont = maxadc -1;
699         if(newcont >= baseline){
700           Warning("","newcont=%f>=baseline=%f",newcont,baseline);
701         } // end if
702           // back to analog: ?
703         fHitMap2->SetHit(i,k,newcont);
704       }  // end for k
705     }else{
706       FastFourierTransform(&fInZR[0],&fInZI[0],1);
707       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
708         Double_t rw = fElectronics->GetTraFunReal(k);
709         Double_t iw = fElectronics->GetTraFunImag(k);
710         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
711         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
712       } // end for k
713       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
714       for(k=0; k<fMaxNofSamples; k++) {
715         Double_t newcont1 = 0.;
716         Double_t maxcont1 = 0.;
717         for(kk=0;kk<nGroup;kk++) {
718           newcont1 = fOutZR[fScaleSize*k+kk];
719           if(newcont1 > maxcont1) maxcont1 = newcont1;
720         } // end for kk
721         newcont1 = maxcont1;
722         if (newcont1 >= maxadc) newcont1 = maxadc -1;
723         fHitMap2->SetHit(i,k,newcont1);
724       } // end for k
725     }
726   } // end for i loop over anodes
727   return;
728 }
729
730 //______________________________________________________________________
731 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
732     // function add the crosstalk effect to signal
733     // temporal function, should be checked...!!!
734   
735     // create and inizialice crosstalk map
736     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
737     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
738     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
739     for( Int_t z=0; z<fNofMaps; z++ ) {
740       Double_t baseline = calibr->GetBaseline(z);
741         Bool_t on = kFALSE;
742         Int_t tstart = 0;
743         Int_t tstop = 0;
744         Int_t nTsteps = 0;
745         
746         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
747             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
748             if( fadc > baseline ) {
749                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
750                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
751                     if( fadc1 < fadc ) continue;
752                     on = kTRUE;
753                     nTsteps = 0;
754                     tstart = l;
755                 }
756                 nTsteps++;
757             }
758             else { // end fadc > baseline
759                 if( on == kTRUE ) {        
760                     if( nTsteps > 2 ) {
761                         tstop = l;
762                         // make smooth derivative
763                         Float_t* dev = new Float_t[fMaxNofSamples+1];
764                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
765                         for( Int_t i=tstart; i<tstop; i++ ) {   
766                             if( i > 2 && i < fMaxNofSamples-2 )
767                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 ) 
768                                     -0.1*fHitMap2->GetSignal( z,i-1 ) 
769                                     +0.1*fHitMap2->GetSignal( z,i+1 ) 
770                                     +0.2*fHitMap2->GetSignal( z,i+2 );
771                         }
772                         
773                         // add crosstalk contribution to neibourg anodes  
774                         for( Int_t i=tstart; i<tstop; i++ ) {
775                             Int_t anode = z - 1;
776                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); // 
777                             Float_t ctktmp =  -dev[i1] * 0.25;
778                             if( anode > 0 ) {
779                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
780                             }
781                             anode = z + 1;
782                             if( anode < fNofMaps ) {
783                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
784                             }
785                         }
786                         delete [] dev;
787                         
788                     } // if( nTsteps > 2 )
789                     on = kFALSE;
790                 }  // if( on == kTRUE )
791             }  // else
792         }
793     }
794     
795     for( Int_t a=0; a<fNofMaps; a++ )
796         for( Int_t t=0; t<fMaxNofSamples; t++ ) {     
797             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
798             fHitMap2->SetHit( a, t, signal );
799         }
800
801     delete [] ctk;
802 }
803
804 //______________________________________________________________________
805 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
806     // To the 10 to 8 bit lossive compression.
807     // code from Davide C. and Albert W.
808
809     if (signal < 128)  return signal;
810     if (signal < 256)  return (128+((signal-128)>>1));
811     if (signal < 512)  return (192+((signal-256)>>3));
812     if (signal < 1024) return (224+((signal-512)>>4));
813     return 0;
814 }
815 //______________________________________________________________________
816 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
817   // Decompression from 8 to 10 bit
818
819   if (signal < 0 || signal > 255) {
820     AliWarning(Form("Signal value %d out of range",signal));
821     return 0;
822   } // end if signal <0 || signal >255
823
824   if (signal < 128) return signal;
825   if (signal < 192) {
826     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
827     else  return (128+((signal-128)<<1)+1);
828   } // end if signal < 192
829   if (signal < 224) {
830     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
831     else  return (256+((signal-192)<<3)+4);
832   } // end if signal < 224
833   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
834   return (512+((signal-224)<<4)+8);
835 }
836 //______________________________________________________________________
837 void AliITSsimulationSDD::Compress2D(){
838   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
839   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);  
840   for (Int_t iWing=0; iWing<2; iWing++) {
841     Int_t tL=res->GetZSLowThreshold(iWing);
842     Int_t tH=res->GetZSHighThreshold(iWing);
843     for (Int_t i=0; i<fNofMaps/2; i++) {  
844       Int_t ian=i+iWing*fNofMaps/2;
845       if( !fAnodeFire[ian] ) continue;
846       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
847         Int_t nLow=0, nHigh=0;      
848         Float_t cC=fHitMap2->GetSignal(ian,itb);
849         if(cC<=tL) continue;
850         nLow++; // cC is greater than tL
851         if(cC>tH) nHigh++;
852         //                     N
853         // Get "quintuple":   WCE
854         //                     S
855         Float_t wW=0.;
856         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
857         if(wW>tL) nLow++;
858         if(wW>tH) nHigh++;
859         Float_t eE=0.;
860         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
861         if(eE>tL) nLow++;
862         if(eE>tH) nHigh++;
863         Float_t nN=0.;
864         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
865         if(nN>tL) nLow++;
866         if(nN>tH) nHigh++;
867         Float_t sS=0.;
868         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
869         if(sS>tL) nLow++;
870         if(sS>tH) nHigh++;
871         
872         if(nLow>=2 && nHigh>=1){
873           Int_t signal=(Int_t)cC;
874           Int_t signalc = Convert10to8(signal);
875           Int_t signale = Convert8to10(signalc);
876           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
877           if(signalc>=4) AddDigit(ian,itb,signalc,signale);  // store C 
878         }
879       }
880     }
881   }
882 }
883
884
885 //______________________________________________________________________
886 void AliITSsimulationSDD::StoreAllDigits(){
887   // store digits for non-zero-suppressed data
888   for (Int_t ian=0; ian<fNofMaps; ian++) {
889     for (Int_t itb=0; itb<fMaxNofSamples; itb++){
890       Int_t signal=(Int_t)(fHitMap2->GetSignal(ian,itb));
891       Int_t signalc = Convert10to8(signal);
892       Int_t signale = Convert8to10(signalc);
893       AddDigit(ian,itb,signalc,signale);  
894     } 
895   }
896
897 //______________________________________________________________________
898 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
899   // Creates histograms of maps for debugging
900   Int_t i;
901   
902   fHis=new TObjArray(fNofMaps);
903   for (i=0;i<fNofMaps;i++) {
904     TString sddName;
905     sddName.Form("sdd_%d",i+1);
906     fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
907                          0.,(Float_t) scale*fMaxNofSamples), i);
908   } // end for i
909 }
910 //______________________________________________________________________
911 void AliITSsimulationSDD::FillHistograms(){
912     // fill 1D histograms from map
913
914     if (!fHis) return;
915
916     for( Int_t i=0; i<fNofMaps; i++) {
917         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
918         Int_t nsamples = hist->GetNbinsX();
919         for( Int_t j=0; j<nsamples; j++) {
920             Double_t signal=fHitMap2->GetSignal(i,j);
921             hist->Fill((Float_t)j,signal);
922         } // end for j
923     } // end for i
924 }
925 //______________________________________________________________________
926 void AliITSsimulationSDD::ResetHistograms(){
927     // Reset histograms for this detector
928     Int_t i;
929
930     for (i=0;i<fNofMaps;i++ ) {
931         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
932     } // end for i
933 }
934 //______________________________________________________________________
935 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) { 
936     // Fills a histogram from a give anode.  
937
938     if (!fHis) return 0;
939
940     if(wing <=0 || wing > 2) {
941         Warning("GetAnode","Wrong wing number: %d",wing);
942         return NULL;
943     } // end if wing <=0 || wing >2
944     if(anode <=0 || anode > fNofMaps/2) {
945         Warning("GetAnode","Wrong anode number: %d",anode);
946         return NULL;
947     } // end if ampde <=0 || andoe > fNofMaps/2
948
949     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
950     return (TH1F*)(fHis->At(index));
951 }
952 //______________________________________________________________________
953 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
954     // Writes the histograms to a file
955
956     if (!fHis) return;
957
958     hfile->cd();
959     Int_t i;
960     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
961     return;
962 }
963 //______________________________________________________________________
964 void AliITSsimulationSDD::WriteSDigits(){
965     // Fills the Summable digits Tree
966     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
967
968     for( Int_t i=0; i<fNofMaps; i++ ) {
969         if( !fAnodeFire[i] ) continue;
970         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
971             Double_t sig = fHitMap2->GetSignal( i, j );
972             if( sig > 0.2 ) {
973                 Int_t jdx = j*fScaleSize;
974                 Int_t index = fpList->GetHitIndex( i, j );
975                 AliITSpListItem pItemTmp2( fModule, index, 0. );
976                 // put the fScaleSize analog digits in only one
977                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
978                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
979                     if( pItemTmp == 0 ) continue;
980                     pItemTmp2.Add( pItemTmp );
981                 }
982                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
983                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
984                 aliITS->AddSumDigit( pItemTmp2 );
985             } // end if (sig > 0.2)
986         }
987     }
988     return;
989 }
990 //______________________________________________________________________
991 void AliITSsimulationSDD::PrintStatus() const {
992     // Print SDD simulation Parameters
993
994     cout << "**************************************************" << endl;
995     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
996     cout << "**************************************************" << endl;
997     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
998     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
999     cout << "Number of Anodes used: " << fNofMaps << endl;
1000     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1001     cout << "Scale size factor: " << fScaleSize << endl;
1002     cout << "**************************************************" << endl;
1003 }