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