]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ITS/AliITSsimulationSDD.cxx
One instance of AliITSresponse per sensor - AliITSresponse objects read/stored from...
[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 <string.h>
22
23 #include <TCanvas.h>
24 #include <TF1.h>
25 #include <TH1.h>
26 #include <TFile.h>
27
28 #include "AliITS.h"
29 #include "AliITSMapA2.h"
30 #include "AliITSRawData.h"
31 #include "AliITSdigitSPD.h"
32 #include "AliITSetfSDD.h"
33 #include "AliITSgeom.h"
34 #include "AliITSmodule.h"
35 #include "AliITSpList.h"
36 #include "AliITSresponseSDD.h"
37 #include "AliITSsegmentationSDD.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 Int_t power(Int_t b, Int_t e) {
53     // compute b to the e power, where both b and e are Int_ts.
54     Int_t power = 1,i;
55
56     for(i=0; i<e; i++) power *= b;
57     return power;
58 }
59 //______________________________________________________________________
60 void FastFourierTransform(AliITSetfSDD *alisddetf,Double_t *real,
61                           Double_t *imag,Int_t direction) {
62     // Do a Fast Fourier Transform
63
64     Int_t samples = alisddetf->GetSamples();
65     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
66     Int_t m1 = samples;
67     Int_t m  = samples/2;
68     Int_t m2 = samples/m1;
69     Int_t i,j,k;
70     for(i=1; i<=l; i++) {
71         for(j=0; j<samples; j += m1) {
72             Int_t p = 0;
73             for(k=j; k<= j+m-1; k++) {
74                 Double_t wsr = alisddetf->GetWeightReal(p);
75                 Double_t wsi = alisddetf->GetWeightImag(p);
76                 if(direction == -1) wsi = -wsi;
77                 Double_t xr = *(real+k+m);
78                 Double_t xi = *(imag+k+m);
79                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
80                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
81                 *(real+k) += xr;
82                 *(imag+k) += xi;
83                 p += m2;
84             } // end for k
85         } // end for j
86         m1 = m;
87         m /= 2;
88         m2 += m2;
89     } // end for i
90   
91     for(j=0; j<samples; j++) {
92         Int_t j1 = j;
93         Int_t p = 0;
94         Int_t i1;
95         for(i1=1; i1<=l; i1++) {
96             Int_t j2 = j1;
97             j1 /= 2;
98             p = p + p + j2 - j1 - j1;
99         } // end for i1
100         if(p >= j) {
101             Double_t xr = *(real+j);
102             Double_t xi = *(imag+j);
103             *(real+j) = *(real+p);
104             *(imag+j) = *(imag+p);
105             *(real+p) = xr;
106             *(imag+p) = xi;
107         } // end if p>=j
108     } // end for j
109     if(direction == -1) {
110         for(i=0; i<samples; i++) {
111             *(real+i) /= samples;
112             *(imag+i) /= samples;
113         } // end for i
114     } // end if direction == -1
115     return;
116 }
117 //______________________________________________________________________
118 AliITSsimulationSDD::AliITSsimulationSDD():
119 AliITSsimulation(),
120 fITS(0),
121 fHitMap2(0),
122 fHitSigMap2(0),
123 fHitNoiMap2(0),
124 fStream(0),
125 fElectronics(0),
126 fInZR(0),
127 fInZI(0),
128 fOutZR(0),
129 fOutZI(0),
130 fAnodeFire(0),
131 fHis(0),
132 fD(),
133 fT1(),
134 fTol(),
135 fBaseline(),
136 fNoise(),
137 fTreeB(0),
138 fFileName(),
139 fFlag(kFALSE),
140 fCheckNoise(kFALSE),
141 fCrosstalkFlag(kFALSE),
142 fDoFFT(1),
143 fNofMaps(0),
144 fMaxNofSamples(0),
145 fScaleSize(0){
146     // Default constructor
147     SetScaleFourier();
148     SetPerpendTracksFlag();
149     SetCrosstalkFlag();
150     SetDoFFT();
151     SetCheckNoise();
152 }
153 //______________________________________________________________________
154 AliITSsimulationSDD::AliITSsimulationSDD(AliITSsimulationSDD &source) :
155     AliITSsimulation(source){
156     // Copy constructor to satify Coding roules only.
157
158     if(this==&source) return;
159     Error("AliITSsimulationSDD","Not allowed to make a copy of "
160           "AliITSsimulationSDD Using default creater instead");
161     AliITSsimulationSDD();
162 }
163 //______________________________________________________________________
164 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
165     // Assignment operator to satify Coding roules only.
166
167     if(this==&src) return *this;
168     Error("AliITSsimulationSDD","Not allowed to make a = with "
169           "AliITSsimulationSDD Using default creater instead");
170     return *this ;
171 }
172 //______________________________________________________________________
173 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
174     // Assignment operator to satify Coding roules only.
175
176     if(this==&src) return *this;
177     Error("AliITSsimulationSSD","Not allowed to make a = with "
178           "AliITSsimulationSDD Using default creater instead");
179     return *this ;
180 }
181
182 //______________________________________________________________________
183 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
184 AliITSsimulation(dettyp),
185 fITS(0),
186 fHitMap2(0),
187 fHitSigMap2(0),
188 fHitNoiMap2(0),
189 fStream(0),
190 fElectronics(0),
191 fInZR(0),
192 fInZI(0),
193 fOutZR(0),
194 fOutZI(0),
195 fAnodeFire(0),
196 fHis(0),
197 fD(),
198 fT1(),
199 fTol(),
200 fBaseline(),
201 fNoise(),
202 fTreeB(0),
203 fFileName(),
204 fFlag(kFALSE),
205 fCheckNoise(kFALSE),
206 fCrosstalkFlag(kFALSE),
207 fDoFFT(1),
208 fNofMaps(0),
209 fMaxNofSamples(0),
210 fScaleSize(0){
211     // Standard Constructor
212
213     Init();
214 }
215 //______________________________________________________________________
216 void AliITSsimulationSDD::Init(){
217     // Standard Constructor
218
219     SetScaleFourier();
220     SetPerpendTracksFlag();
221     SetCrosstalkFlag();
222     SetDoFFT();
223     SetCheckNoise();
224
225     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
226     
227     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
228     fpList = new AliITSpList( seg->Npz(),
229                               fScaleSize*seg->Npx() );
230     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
231     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
232     fHitMap2 = fHitSigMap2;
233
234     fNofMaps = seg->Npz();
235     fMaxNofSamples = seg->Npx();
236     fAnodeFire = new Bool_t [fNofMaps];
237     
238     Float_t sddLength = seg->Dx();
239     Float_t sddWidth  = seg->Dz();
240
241     Int_t dummy        = 0;
242     Float_t anodePitch = seg->Dpz(dummy);
243     Double_t timeStep  = (Double_t)seg->Dpx(dummy);
244     Float_t driftSpeed = res->DriftSpeed();
245
246     if(anodePitch*(fNofMaps/2) > sddWidth) {
247         Warning("AliITSsimulationSDD",
248                 "Too many anodes %d or too big pitch %f \n",
249                 fNofMaps/2,anodePitch);
250     } // end if
251
252     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
253         Error("AliITSsimulationSDD",
254               "Time Interval > Allowed Time Interval: exit\n");
255         return;
256     } // end if
257
258     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
259                                     res->Electronics());
260
261     char opt1[20], opt2[20];
262     res->ParamOptions(opt1,opt2);
263     fParam = opt2;
264     char *same = strstr(opt1,"same");
265     if (same) {
266         fNoise.Set(0);
267         fBaseline.Set(0);
268     } else {
269         fNoise.Set(fNofMaps);
270         fBaseline.Set(fNofMaps);
271     } // end if
272
273     const char *kopt=res->ZeroSuppOption();
274     if (strstr(fParam.Data(),"file") ) {
275         fD.Set(fNofMaps);
276         fT1.Set(fNofMaps);
277         if (strstr(kopt,"2D")) {
278             fT2.Set(fNofMaps);
279             fTol.Set(0);
280             Init2D();       // desactivate if param change module by module
281         } else if(strstr(kopt,"1D"))  {
282             fT2.Set(2);
283             fTol.Set(2);
284             Init1D();      // desactivate if param change module by module
285         } // end if strstr
286     } else {
287         fD.Set(2);
288         fTol.Set(2);
289         fT1.Set(2);
290         fT2.Set(2);
291         SetCompressParam();
292     } // end if else strstr
293
294     Bool_t write = res->OutputOption();
295     if(write && strstr(kopt,"2D")) MakeTreeB();
296
297     // call here if baseline does not change by module
298     // ReadBaseline();
299
300     fITS       = (AliITS*)gAlice->GetModule("ITS");
301     Int_t size = fNofMaps*fMaxNofSamples;
302     fStream    = new AliITSInStream(size);
303
304     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
305     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
306     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
307     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];  
308 }
309 //______________________________________________________________________
310 AliITSsimulationSDD::~AliITSsimulationSDD() { 
311     // destructor
312
313     //    delete fpList;
314     delete fHitSigMap2;
315     delete fHitNoiMap2;
316     delete fStream;
317     delete fElectronics;
318
319     fITS = 0;
320
321     if (fHis) {
322         fHis->Delete(); 
323         delete fHis;     
324     } // end if fHis
325     if(fTreeB) delete fTreeB;           
326     if(fInZR)  delete [] fInZR;
327     if(fInZI)  delete [] fInZI;        
328     if(fOutZR) delete [] fOutZR;
329     if(fOutZI) delete [] fOutZI;
330     if(fAnodeFire) delete [] fAnodeFire;
331 }
332 //______________________________________________________________________
333 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
334     // create maps to build the lists of tracks for each summable digit
335     fModule = module;
336     fEvent  = event;
337     ClearMaps();
338     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);    
339 }
340 //______________________________________________________________________
341 void AliITSsimulationSDD::ClearMaps() {
342     // clear maps
343     fpList->ClearMap();
344     fHitSigMap2->ClearMap();
345     fHitNoiMap2->ClearMap();
346 }
347 //______________________________________________________________________
348 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
349     // digitize module using the "slow" detector simulator creating
350     // summable digits.
351
352     TObjArray *fHits = mod->GetHits();
353     Int_t nhits      = fHits->GetEntriesFast();
354     if( !nhits ) return;
355
356     InitSimulationModule( md, ev );
357     HitsToAnalogDigits( mod );
358     ChargeToSignal( fModule,kFALSE ); // - Process signal without add noise
359     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
360     ChargeToSignal( fModule,kTRUE );  // - Process only noise
361     fHitMap2 = fHitSigMap2;   // - Return to signal map
362     WriteSDigits();
363     ClearMaps();
364 }
365 //______________________________________________________________________
366 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
367                                                Int_t mask ) {
368     // Add Summable digits to module maps.
369   AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
370     Int_t    nItems = pItemArray->GetEntries();
371     Double_t maxadc = res->MaxAdc();
372     Bool_t sig = kFALSE;
373     
374     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
375     for( Int_t i=0; i<nItems; i++ ) {
376         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
377         if( pItem->GetModule() != fModule ) {
378             Error( "AliITSsimulationSDD","Error reading, SDigits module "
379                    "%d != current module %d: exit",
380                    pItem->GetModule(), fModule );
381             return sig;
382         } // end if
383
384         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
385         
386         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
387         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
388         Double_t sigAE = pItem2->GetSignalAfterElect();
389         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
390         Int_t ia;
391         Int_t it;
392         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
393         fHitMap2->SetHit( ia, it, sigAE );
394         fAnodeFire[ia] = kTRUE;
395     }
396     return sig;
397 }
398 //______________________________________________________________________
399 void AliITSsimulationSDD::FinishSDigitiseModule() {
400     // digitize module using the "slow" detector simulator from
401     // the sum of summable digits.
402     FinishDigits() ;
403     ClearMaps();
404 }
405 //______________________________________________________________________
406 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
407     // create maps to build the lists of tracks for each digit
408
409     TObjArray *fHits = mod->GetHits();
410     Int_t nhits      = fHits->GetEntriesFast();
411
412     InitSimulationModule( md, ev );
413
414     if( !nhits && fCheckNoise ) {
415         ChargeToSignal( fModule,kTRUE );  // process noise
416         GetNoise();
417         ClearMaps();
418         return;
419     } else 
420         if( !nhits ) return;
421         
422     HitsToAnalogDigits( mod );
423     ChargeToSignal( fModule,kTRUE );  // process signal + noise
424
425     for( Int_t i=0; i<fNofMaps; i++ ) {
426         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
427             Int_t jdx = j*fScaleSize;
428             Int_t index = fpList->GetHitIndex( i, j );
429             AliITSpListItem pItemTmp2( fModule, index, 0. );
430             // put the fScaleSize analog digits in only one
431             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
432                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
433                 if( pItemTmp == 0 ) continue;
434                 pItemTmp2.Add( pItemTmp );
435             }
436             fpList->DeleteHit( i, j );
437             fpList->AddItemTo( 0, &pItemTmp2 );
438         }
439     }
440     FinishDigits();
441     ClearMaps();
442 }
443 //______________________________________________________________________
444 void AliITSsimulationSDD::FinishDigits() {
445     // introduce the electronics effects and do zero-suppression if required
446
447     ApplyDeadChannels(fModule);
448     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
449
450     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
451     const char *kopt = res->ZeroSuppOption();
452     ZeroSuppression( kopt );
453 }
454 //______________________________________________________________________
455 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
456     // create maps to build the lists of tracks for each digit
457
458   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
459   AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
460   TObjArray *hits     = mod->GetHits();
461     Int_t      nhits    = hits->GetEntriesFast();
462
463     //    Int_t      arg[6]   = {0,0,0,0,0,0};
464     Int_t    dummy      = 0;
465     Int_t    nofAnodes  = fNofMaps/2;
466     Float_t  sddLength  = seg->Dx();
467     Float_t  sddWidth   = seg->Dz();
468     Float_t  anodePitch = seg->Dpz(dummy);
469     Float_t  timeStep   = seg->Dpx(dummy);
470     Float_t  driftSpeed = res->DriftSpeed();
471     Float_t  maxadc     = res->MaxAdc();    
472     Float_t  topValue   = res->DynamicRange();
473     Float_t  cHloss     = res->ChargeLoss();
474     Float_t  norm       = maxadc/topValue;
475     Double_t dfCoeff, s1; res->DiffCoeff(dfCoeff,s1); // Signal 2d Shape
476     Double_t eVpairs    = res->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
477     Float_t  nsigma     = res->NSigmaIntegration(); //
478     Int_t    nlookups   = res->GausNLookUp();       //
479     Float_t  jitter     = res->JitterError(); // 
480
481     // Piergiorgio's part (apart for few variables which I made float
482     // when i thought that can be done
483     // Fill detector maps with GEANT hits
484     // loop over hits in the module
485
486     const Float_t kconv = 1.0e+6;  // GeV->KeV
487     Int_t    itrack      = 0;
488     Int_t    hitDetector; // detector number (lay,lad,hitDetector)
489     Int_t    iWing;       // which detector wing/side.
490     Int_t    detector;    // 2*(detector-1)+iWing
491     Int_t    ii,kk,ka,kt; // loop indexs
492     Int_t    ia,it,index; // sub-pixel integration indexies
493     Int_t    iAnode;      // anode number.
494     Int_t    timeSample;  // time buckett.
495     Int_t    anodeWindow; // anode direction charge integration width
496     Int_t    timeWindow;  // time direction charge integration width
497     Int_t    jamin,jamax; // anode charge integration window
498     Int_t    jtmin,jtmax; // time charge integration window
499     Int_t    ndiv;        // Anode window division factor.
500     Int_t    nsplit;      // the number of splits in anode and time windows==1.
501     Int_t    nOfSplits;   // number of times track length is split into
502     Float_t  nOfSplitsF;  // Floating point version of nOfSplits.
503     Float_t  kkF;         // Floating point version of loop index kk.
504     Float_t  pathInSDD; // Track length in SDD.
505     Float_t  drPath; // average position of track in detector. in microns
506     Float_t  drTime; // Drift time
507     Float_t  nmul;   // drift time window multiplication factor.
508     Float_t  avDrft;  // x position of path length segment in cm.
509     Float_t  avAnode; // Anode for path length segment in Anode number (float)
510     Float_t  xAnode;  // Floating point anode number.
511     Float_t  driftPath; // avDrft in microns.
512     Float_t  width;     // width of signal at anodes.
513     Double_t  depEnergy; // Energy deposited in this GEANT step.
514     Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
515     Double_t sigA; // sigma of signal at anode.
516     Double_t sigT; // sigma in time/drift direction for track segment
517     Double_t aStep,aConst; // sub-pixel size and offset anode
518     Double_t tStep,tConst; // sub-pixel size and offset time
519     Double_t amplitude; // signal amplitude for track segment in nanoAmpere
520     Double_t chargeloss; // charge loss for track segment.
521     Double_t anodeAmplitude; // signal amplitude in anode direction
522     Double_t aExpo;          // exponent of Gaussian anode direction
523     Double_t timeAmplitude;  // signal amplitude in time direction
524     Double_t tExpo;          // exponent of Gaussian time direction
525     //  Double_t tof;            // Time of flight in ns of this step.    
526
527     for(ii=0; ii<nhits; ii++) {
528         if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
529                               depEnergy,itrack)) continue;
530         xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
531         depEnergy  *= kconv;
532         hitDetector = mod->GetDet();
533         //tof         = 1.E+09*(mod->GetHit(ii)->GetTOF()); // tof in ns.
534         //if(tof>sddLength/driftSpeed) continue; // hit happed too late.
535         
536         // scale path to simulate a perpendicular track
537         // continue if the particle did not lose energy
538         // passing through detector
539         if (!depEnergy) {
540             AliDebug(1,
541              Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
542                   itrack,ii,mod->GetIndex()));
543             continue;
544         } // end if !depEnergy
545
546         pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
547
548         if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
549         drPath = 10000.*(dxL[0]+2.*xL[0])*0.5;
550         if(drPath < 0) drPath = -drPath;
551         drPath = sddLength-drPath;
552         if(drPath < 0) {
553             AliDebug(1, // this should be fixed at geometry level
554               Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
555                    drPath,sddLength,dxL[0],xL[0]));
556             continue;
557         } // end if drPath < 0
558
559         // Compute number of segments to brake step path into
560         drTime = drPath/driftSpeed;  //   Drift Time
561         sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
562         // calcuate the number of time the path length should be split into.
563         nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
564         if(fFlag) nOfSplits = 1;
565
566         // loop over path segments, init. some variables.
567         depEnergy /= nOfSplits;
568         nOfSplitsF = (Float_t) nOfSplits;
569         for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
570             kkF       = (Float_t) kk + 0.5;
571             avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
572             avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
573             driftPath = 10000.*avDrft;
574
575             iWing = 2;  // Assume wing is 2
576             if(driftPath < 0) { // if wing is not 2 it is 1.
577                 iWing     = 1;
578                 driftPath = -driftPath;
579             } // end if driftPath < 0
580             driftPath = sddLength-driftPath;
581             detector  = 2*(hitDetector-1) + iWing;
582             if(driftPath < 0) {
583                 AliDebug(1, // this should be fixed at geometry level
584                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
585                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
586                 continue;
587             } // end if driftPath < 0
588
589             //   Drift Time
590             drTime     = driftPath/driftSpeed; // drift time for segment.
591             timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1);
592             // compute time Sample including tof information. The tof only 
593             // effects the time of the signal is recoreded and not the
594             // the defusion.
595             // timeSample = (Int_t) (fScaleSize*(drTime+tof)/timeStep + 1);
596             if(timeSample > fScaleSize*fMaxNofSamples) {
597                 Warning("HitsToAnalogDigits","Wrong Time Sample: %e",
598                         timeSample);
599                 continue;
600             } // end if timeSample > fScaleSize*fMaxNoofSamples
601
602             //   Anode
603             xAnode = 10000.*(avAnode)/anodePitch + nofAnodes/2;  // +1?
604             if(xAnode*anodePitch > sddWidth || xAnode*anodePitch < 0.) 
605                 Warning("HitsToAnalogDigits",
606                         "Exceedubg sddWidth=%e Z = %e",
607                         sddWidth,xAnode*anodePitch);
608             iAnode = (Int_t) (1.+xAnode); // xAnode?
609             if(iAnode < 1 || iAnode > nofAnodes) {
610                 Warning("HitToAnalogDigits","Wrong iAnode: 1<%d>%d",
611                         iAnode,nofAnodes);
612                 continue;
613             } // end if iAnode < 1 || iAnode > nofAnodes
614
615             // store straight away the particle position in the array
616             // of particles and take idhit=ii only when part is entering (this
617             // requires FillModules() in the macro for analysis) :
618     
619             // Sigma along the anodes for track segment.
620             sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
621             sigT       = sigA/driftSpeed;
622             // Peak amplitude in nanoAmpere
623             amplitude  = fScaleSize*160.*depEnergy/
624                 (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA);
625             amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to 
626             // account for clock variations 
627             // (reference value: 40 MHz)
628             chargeloss = 1.-cHloss*driftPath/1000;
629             amplitude *= chargeloss;
630             width  = 2.*nsigma/(nlookups-1);
631             // Spread the charge 
632             // Pixel index
633             ndiv = 2;
634             nmul = 3.; 
635             if(drTime > 1200.) { 
636                 ndiv = 4;
637                 nmul = 1.5;
638             } // end if drTime > 1200.
639             // Sub-pixel index
640             nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2;
641             // Sub-pixel size see computation of aExpo and tExpo.
642             aStep  = anodePitch/(nsplit*fScaleSize*sigA);
643             aConst = xAnode*anodePitch/sigA;
644             tStep  = timeStep/(nsplit*fScaleSize*sigT);
645             tConst = drTime/sigT;
646             // Define SDD window corresponding to the hit
647             anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1);
648             timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
649             jamin = (iAnode - anodeWindow/ndiv - 1)*fScaleSize*nsplit +1;
650             jamax = (iAnode + anodeWindow/ndiv)*fScaleSize*nsplit;
651             if(jamin <= 0) jamin = 1;
652             if(jamax > fScaleSize*nofAnodes*nsplit) 
653                 jamax = fScaleSize*nofAnodes*nsplit;
654             // jtmin and jtmax are Hard-wired
655             jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1;
656             jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit;
657             if(jtmin <= 0) jtmin = 1;
658             if(jtmax > fScaleSize*fMaxNofSamples*nsplit) 
659                 jtmax = fScaleSize*fMaxNofSamples*nsplit;
660             // Spread the charge in the anode-time window
661             for(ka=jamin; ka <=jamax; ka++) {
662                 ia = (ka-1)/(fScaleSize*nsplit) + 1;
663                 if(ia <= 0) {
664                     Warning("HitsToAnalogDigits","ia < 1: ");
665                     continue;
666                 } // end if
667                 if(ia > nofAnodes) ia = nofAnodes;
668                 aExpo     = (aStep*(ka-0.5)-aConst);
669                 if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
670                 else {
671                     dummy          = (Int_t) ((aExpo+nsigma)/width);
672                     anodeAmplitude = amplitude*res->GausLookUp(dummy);
673                 } // end if TMath::Abs(aEspo) > nsigma
674                 // index starts from 0
675                 index = ((detector+1)%2)*nofAnodes+ia-1;
676                 if(anodeAmplitude) for(kt=jtmin; kt<=jtmax; kt++) {
677                     it = (kt-1)/nsplit+1;  // it starts from 1
678                     if(it<=0){
679                         Warning("HitsToAnalogDigits","it < 1:");
680                         continue;
681                     } // end if 
682                     if(it>fScaleSize*fMaxNofSamples)
683                         it = fScaleSize*fMaxNofSamples;
684                     tExpo    = (tStep*(kt-0.5)-tConst);
685                     if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
686                     else {
687                         dummy         = (Int_t) ((tExpo+nsigma)/width);
688                         timeAmplitude = anodeAmplitude*
689                             res->GausLookUp(dummy);
690                     } // end if TMath::Abs(tExpo) > nsigma
691                     // build the list of Sdigits for this module        
692                     //                    arg[0]     = index;
693                     //                    arg[1]     = it;
694                     //                    arg[2]     = itrack; // track number
695                     //                    arg[3]     = ii-1; // hit number.
696                     timeAmplitude *= norm;
697                     timeAmplitude *= 10;
698                     //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
699                     Double_t charge = timeAmplitude;
700                     charge += fHitMap2->GetSignal(index,it-1);
701                     fHitMap2->SetHit(index, it-1, charge);
702                     fpList->AddSignal(index,it-1,itrack,ii-1,
703                                       mod->GetIndex(),timeAmplitude);
704                     fAnodeFire[index] = kTRUE;                 
705                 } // end if anodeAmplitude and loop over time in window
706             } // loop over anodes in window
707         } // end loop over "sub-hits"
708     } // end loop over hits
709 }
710 /*
711 //______________________________________________________________________
712 void AliITSsimulationSDD::ListOfFiredCells(Int_t *arg,Double_t timeAmplitude,
713                                            TObjArray *alist,TClonesArray *padr){
714     // Returns the list of "fired" cells.
715
716     Int_t index     = arg[0];
717     Int_t ik        = arg[1];
718     Int_t idtrack   = arg[2];
719     Int_t idhit     = arg[3];
720     Int_t counter   = arg[4];
721     Int_t countadr  = arg[5];
722     Double_t charge = timeAmplitude;
723     charge += fHitMap2->GetSignal(index,ik-1);
724     fHitMap2->SetHit(index, ik-1, charge);
725
726     Int_t digits[3];
727     Int_t it = (Int_t)((ik-1)/fScaleSize);
728     digits[0] = index;
729     digits[1] = it;
730     digits[2] = (Int_t)timeAmplitude;
731     Float_t phys;
732     if (idtrack >= 0) phys = (Float_t)timeAmplitude;
733     else phys = 0;
734
735     Double_t cellcharge = 0.;
736     AliITSTransientDigit* pdigit;
737     // build the list of fired cells and update the info
738     if (!fHitMap1->TestHit(index, it)) {
739         new((*padr)[countadr++]) TVector(3);
740         TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
741         trinfo(0) = (Float_t)idtrack;
742         trinfo(1) = (Float_t)idhit;
743         trinfo(2) = (Float_t)timeAmplitude;
744       
745         alist->AddAtAndExpand(new AliITSTransientDigit(phys,digits),counter);
746         fHitMap1->SetHit(index, it, counter);
747         counter++;
748         pdigit=(AliITSTransientDigit*)alist->At(alist->GetLast());
749         // list of tracks
750         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
751         trlist->Add(&trinfo);
752     } else {
753         pdigit = (AliITSTransientDigit*) fHitMap1->GetHit(index, it);
754         for(Int_t kk=0;kk<fScaleSize;kk++) {
755             cellcharge += fHitMap2->GetSignal(index,fScaleSize*it+kk);
756         }  // end for kk
757         // update charge
758         (*pdigit).fSignal = (Int_t)cellcharge;
759         (*pdigit).fPhysics += phys;                        
760         // update list of tracks
761         TObjArray* trlist = (TObjArray*)pdigit->TrackList();
762         Int_t lastentry = trlist->GetLast();
763         TVector *ptrkp = (TVector*)trlist->At(lastentry);
764         TVector &trinfo = *ptrkp;
765         Int_t lasttrack = Int_t(trinfo(0));
766         Float_t lastcharge=(trinfo(2));
767         if (lasttrack==idtrack ) {
768             lastcharge += (Float_t)timeAmplitude;
769             trlist->RemoveAt(lastentry);
770             trinfo(0) = lasttrack;
771             trinfo(1) = idhit;
772             trinfo(2) = lastcharge;
773             trlist->AddAt(&trinfo,lastentry);
774         } else {                  
775             new((*padr)[countadr++]) TVector(3);
776             TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
777             trinfo(0) = (Float_t)idtrack;
778             trinfo(1) = (Float_t)idhit;
779             trinfo(2) = (Float_t)timeAmplitude;
780             trlist->Add(&trinfo);
781         } // end if lasttrack==idtrack
782
783         if(AliDebugLevel()){
784             // check the track list - debugging
785             Int_t trk[20], htrk[20];
786             Float_t chtrk[20];  
787             Int_t nptracks = trlist->GetEntriesFast();
788             if (nptracks > 2) {
789                 Int_t tr;
790                 for (tr=0;tr<nptracks;tr++) {
791                     TVector *pptrkp = (TVector*)trlist->At(tr);
792                     TVector &pptrk  = *pptrkp;
793                     trk[tr]   = Int_t(pptrk(0));
794                     htrk[tr]  = Int_t(pptrk(1));
795                     chtrk[tr] = (pptrk(2));
796                     cout << "nptracks "<<nptracks << endl;
797                 } // end for tr
798             } // end if nptracks
799         } // end if AliDebugLevel()
800     } //  end if pdigit
801
802     // update counter and countadr for next call.
803     arg[4] = counter;
804     arg[5] = countadr;
805 }
806 */
807 //____________________________________________
808 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signal ) {
809     // Adds a Digit.
810     Int_t size = AliITSdigitSPD::GetNTracks();
811     Int_t digits[3];
812     Int_t * tracks = new Int_t[size];
813     Int_t * hits = new Int_t[size];
814     Float_t phys;
815     Float_t * charges = new Float_t[size];
816
817     digits[0] = i;
818     digits[1] = j;
819     digits[2] = signal;
820
821     AliITSpListItem *pItem = fpList->GetpListItem( i, j );
822     if( pItem == 0 ) {
823         phys = 0.0;
824         for( Int_t l=0; l<size; l++ ) {
825             tracks[l]  = 0;
826             hits[l]    = 0;
827             charges[l] = 0.0;
828         }
829     } else {
830         Int_t idtrack =  pItem->GetTrack( 0 );
831         if( idtrack >= 0 ) phys = pItem->GetSignal();  
832         else phys = 0.0;
833
834         for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
835             tracks[l]  = pItem->GetTrack( l );
836             hits[l]    = pItem->GetHit( l );
837             charges[l] = pItem->GetSignal( l );
838         }else{
839             tracks[l]  = -3;
840             hits[l]    = -1;
841             charges[l] = 0.0;
842         }// end for if
843     }
844
845     fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges ); 
846     delete [] tracks;
847     delete [] hits;
848     delete [] charges;
849 }
850 //______________________________________________________________________
851 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise) {
852     // add baseline, noise, electronics and ADC saturation effects
853
854     char opt1[20], opt2[20];
855     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(mod);
856
857     res->ParamOptions(opt1,opt2);
858     char *read = strstr(opt1,"file");
859     Double_t baseline, noise; 
860
861     if (read) {
862         static Bool_t readfile=kTRUE;
863         //read baseline and noise from file
864         if (readfile) ReadBaseline();
865         readfile=kFALSE;
866     } else res->GetNoiseParam(noise,baseline);
867
868     Float_t contrib=0;
869     Int_t i,k,kk;
870     Float_t maxadc = res->MaxAdc();    
871     if(!fDoFFT) {
872         for (i=0;i<fNofMaps;i++) {
873             if( !fAnodeFire[i] ) continue;
874             if (read && i<fNofMaps) GetAnodeBaseline(i,baseline,noise);
875             for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
876                 fInZR[k]  = fHitMap2->GetSignal(i,k);
877                 if( bAddNoise ) {
878                     contrib   = (baseline + noise*gRandom->Gaus());
879                     fInZR[k] += contrib;
880                 }
881             } // end for k
882             for(k=0; k<fMaxNofSamples; k++) {
883                 Double_t newcont = 0.;
884                 Double_t maxcont = 0.;
885                 for(kk=0;kk<fScaleSize;kk++) {
886                     newcont = fInZR[fScaleSize*k+kk];
887                     if(newcont > maxcont) maxcont = newcont;
888                 } // end for kk
889                 newcont = maxcont;
890                 if (newcont >= maxadc) newcont = maxadc -1;
891                 if(newcont >= baseline){
892                     Warning("","newcont=%d>=baseline=%d",newcont,baseline);
893                 } // end if
894                 // back to analog: ?
895                 fHitMap2->SetHit(i,k,newcont);
896             }  // end for k
897         } // end for i loop over anodes
898         return;
899     } // end if DoFFT
900
901     for (i=0;i<fNofMaps;i++) {
902         if( !fAnodeFire[i] ) continue;
903         if  (read && i<fNofMaps) GetAnodeBaseline(i,baseline,noise);
904         for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
905             fInZR[k]  = fHitMap2->GetSignal(i,k);
906             if( bAddNoise ) {
907                 contrib   = (baseline + noise*gRandom->Gaus());
908                 fInZR[k] += contrib;
909             }
910             fInZI[k]  = 0.;
911         } // end for k
912         FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
913         for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
914             Double_t rw = fElectronics->GetTraFunReal(k);
915             Double_t iw = fElectronics->GetTraFunImag(k);
916             fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
917             fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
918         } // end for k
919         FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
920         for(k=0; k<fMaxNofSamples; k++) {
921             Double_t newcont1 = 0.;
922             Double_t maxcont1 = 0.;
923             for(kk=0;kk<fScaleSize;kk++) {
924                 newcont1 = fOutZR[fScaleSize*k+kk];
925                 if(newcont1 > maxcont1) maxcont1 = newcont1;
926             } // end for kk
927             newcont1 = maxcont1;
928             if (newcont1 >= maxadc) newcont1 = maxadc -1;
929             fHitMap2->SetHit(i,k,newcont1);
930         } // end for k
931     } // end for i loop over anodes
932     return;
933 }
934 //____________________________________________________________________
935 void AliITSsimulationSDD::ApplyDeadChannels(Int_t mod) {    
936     // Set dead channel signal to zero
937     AliITSresponseSDD * response = (AliITSresponseSDD *)GetResponseModel(mod);
938     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
939     // nothing to do
940     if( response->IsDead() ||   
941         ( response->GetDeadChips() == 0 &&
942           response->GetDeadChannels() == 0 ) ) return;  
943     
944     // static AliITS *iTS = (AliITS*)gAlice->GetModule( "ITS" );
945
946     Int_t fMaxNofSamples = seg->Npx();    
947     // AliITSgeom *geom = iTS->GetITSgeom();
948     // Int_t firstSDDMod = geom->GetStartDet( 1 );
949     // loop over wings
950     for( Int_t j=0; j<2; j++ ) {
951       // Int_t mod = (fModule-firstSDDMod)*2 + j;
952       for( Int_t u=0; u<response->Chips(); u++ )
953         for( Int_t v=0; v<response->Channels(); v++ ) {
954           Float_t gain = response->Gain(j, u, v );
955           for( Int_t k=0; k<fMaxNofSamples; k++ ) {
956             Int_t i = j*response->Chips()*response->Channels() +
957               u*response->Channels() + 
958               v;
959             Double_t signal =  gain * fHitMap2->GetSignal( i, k );
960             fHitMap2->SetHit( i, k, signal );  ///
961           }
962         }
963     } 
964 }
965 //______________________________________________________________________
966 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
967     // function add the crosstalk effect to signal
968     // temporal function, should be checked...!!!
969     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
970   
971     Int_t fNofMaps = seg->Npz();
972     Int_t fMaxNofSamples = seg->Npx();
973
974     // create and inizialice crosstalk map
975     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
976     if( ctk == NULL ) {
977         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
978         return;
979     }
980     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
981     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(mod);
982
983     Double_t noise, baseline;
984     res->GetNoiseParam( noise, baseline );
985     for( Int_t z=0; z<fNofMaps; z++ ) {
986         Bool_t on = kFALSE;
987         Int_t tstart = 0;
988         Int_t tstop = 0;
989         Int_t nTsteps = 0;
990         
991         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
992             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
993             if( fadc > baseline ) {
994                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
995                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
996                     if( fadc1 < fadc ) continue;
997                     on = kTRUE;
998                     nTsteps = 0;
999                     tstart = l;
1000                 }
1001                 nTsteps++;
1002             }
1003             else { // end fadc > baseline
1004                 if( on == kTRUE ) {        
1005                     if( nTsteps > 2 ) {
1006                         tstop = l;
1007                         // make smooth derivative
1008                         Float_t* dev = new Float_t[fMaxNofSamples+1];
1009                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
1010                         if( ctk == NULL ) {
1011                             Error( "ApplyCrosstalk", 
1012                                    "no memory for temporal array: exit \n" );
1013                             return;
1014                         }
1015                         for( Int_t i=tstart; i<tstop; i++ ) {   
1016                             if( i > 2 && i < fMaxNofSamples-2 )
1017                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 ) 
1018                                     -0.1*fHitMap2->GetSignal( z,i-1 ) 
1019                                     +0.1*fHitMap2->GetSignal( z,i+1 ) 
1020                                     +0.2*fHitMap2->GetSignal( z,i+2 );
1021                         }
1022                         
1023                         // add crosstalk contribution to neibourg anodes  
1024                         for( Int_t i=tstart; i<tstop; i++ ) {
1025                             Int_t anode = z - 1;
1026                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); // 
1027                             Float_t ctktmp =  -dev[i1] * 0.25;
1028                             if( anode > 0 ) {
1029                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
1030                             }
1031                             anode = z + 1;
1032                             if( anode < fNofMaps ) {
1033                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
1034                             }
1035                         }
1036                         delete [] dev;
1037                         
1038                     } // if( nTsteps > 2 )
1039                     on = kFALSE;
1040                 }  // if( on == kTRUE )
1041             }  // else
1042         }
1043     }
1044     
1045     for( Int_t a=0; a<fNofMaps; a++ )
1046         for( Int_t t=0; t<fMaxNofSamples; t++ ) {     
1047             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
1048             fHitMap2->SetHit( a, t, signal );
1049         }
1050
1051     delete [] ctk;
1052 }
1053 //______________________________________________________________________
1054 void AliITSsimulationSDD::GetAnodeBaseline(Int_t i,Double_t &baseline,
1055                                            Double_t &noise) const{
1056     // Returns the Baseline for a particular anode.
1057     baseline = fBaseline[i];
1058     noise    = fNoise[i];
1059 }
1060 //______________________________________________________________________
1061 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
1062                                            Int_t &th){
1063     // Returns the compression alogirthm parameters
1064     Int_t size = fD.GetSize();
1065     if (size > 2 ) {
1066         db=fD[i]; tl=fT1[i]; th=fT2[i];
1067     } else {
1068         if (size <= 2 && i>=fNofMaps/2) {
1069             db=fD[1]; tl=fT1[1]; th=fT2[1];
1070         } else {
1071             db=fD[0]; tl=fT1[0]; th=fT2[0];
1072         } // end if size <=2 && i>=fNofMaps/2
1073     } // end if size >2
1074 }
1075 //______________________________________________________________________
1076 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl){
1077     // returns the compression alogirthm parameters
1078     Int_t size = fD.GetSize();
1079
1080     if (size > 2 ) {
1081         db=fD[i]; tl=fT1[i];
1082     } else {
1083         if (size <= 2 && i>=fNofMaps/2) {
1084             db=fD[1]; tl=fT1[1]; 
1085         } else {
1086             db=fD[0]; tl=fT1[0]; 
1087         } // end if size <=2 && i>=fNofMaps/2
1088         //      Warning("CompressionParam",
1089         //             "Size= %d . Values i=%d ; db= %d ; tl= %d",
1090         //              size,i,db,tl);
1091     } // end if size > 2
1092 }
1093 //______________________________________________________________________
1094 void AliITSsimulationSDD::SetCompressParam(){
1095     // Sets the compression alogirthm parameters  
1096     Int_t cp[8],i;
1097     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1098
1099     res->GiveCompressParam(cp);
1100     for (i=0; i<2; i++) {
1101         fD[i]   = cp[i];
1102         fT1[i]  = cp[i+2];
1103         fT2[i]  = cp[i+4];
1104         fTol[i] = cp[i+6];
1105     } // end for i
1106 }
1107 //______________________________________________________________________
1108 void AliITSsimulationSDD::ReadBaseline(){
1109     // read baseline and noise from file - either a .root file and in this
1110     // case data should be organised in a tree with one entry for each
1111     // module => reading should be done accordingly
1112     // or a classic file and do smth. like this:
1113     // Read baselines and noise for SDD
1114
1115     Int_t na,pos;
1116     Float_t bl,n;
1117     char input[100], base[100], param[100];
1118     char *filtmp;
1119     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1120     res->Filenames(input,base,param);
1121     fFileName=base;
1122     //
1123     filtmp = gSystem->ExpandPathName(fFileName.Data());
1124     FILE *bline = fopen(filtmp,"r");
1125     na = 0;
1126
1127     if(bline) {
1128         while(fscanf(bline,"%d %f %f",&pos, &bl, &n) != EOF) {
1129             if (pos != na+1) {
1130                 Error("ReadBaseline","Anode number not in increasing order!",
1131                       filtmp);
1132                 exit(1);
1133             } // end if pos != na+1
1134             fBaseline[na]=bl;
1135             fNoise[na]=n;
1136             na++;
1137         } // end while
1138     } else {
1139         Error("ReadBaseline"," THE BASELINE FILE %s DOES NOT EXIST !",filtmp);
1140         exit(1);
1141     } // end if(bline)
1142
1143     fclose(bline);
1144     delete [] filtmp;
1145 }
1146 //______________________________________________________________________
1147 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
1148     // To the 10 to 8 bit lossive compression.
1149     // code from Davide C. and Albert W.
1150
1151     if (signal < 128)  return signal;
1152     if (signal < 256)  return (128+((signal-128)>>1));
1153     if (signal < 512)  return (192+((signal-256)>>3));
1154     if (signal < 1024) return (224+((signal-512)>>4));
1155     return 0;
1156 }
1157 /*
1158 //______________________________________________________________________
1159 AliITSMap*   AliITSsimulationSDD::HitMap(Int_t i){
1160     //Return the correct map.
1161
1162     return ((i==0)? fHitMap1 : fHitMap2);
1163 }
1164 */
1165 //______________________________________________________________________
1166 void AliITSsimulationSDD::ZeroSuppression(const char *option) {
1167     // perform the zero suppresion
1168
1169     if (strstr(option,"2D")) {
1170         //Init2D();              // activate if param change module by module
1171         Compress2D();
1172     } else if (strstr(option,"1D")) {
1173         //Init1D();              // activate if param change module by module
1174         Compress1D();  
1175     } else StoreAllDigits();
1176 }
1177 //______________________________________________________________________
1178 void AliITSsimulationSDD::Init2D(){
1179     // read in and prepare arrays: fD, fT1, fT2
1180     //                         savemu[nanodes], savesigma[nanodes] 
1181     // read baseline and noise from file - either a .root file and in this
1182     // case data should be organised in a tree with one entry for each
1183     // module => reading should be done accordingly
1184     // or a classic file and do smth. like this ( code from Davide C. and
1185     // Albert W.) :
1186     // Read 2D zero-suppression parameters for SDD
1187
1188     if (!strstr(fParam.Data(),"file")) return;
1189
1190     Int_t na,pos,tempTh;
1191     Float_t mu,sigma;
1192     Float_t *savemu    = new Float_t [fNofMaps];
1193     Float_t *savesigma = new Float_t [fNofMaps];
1194     char input[100],basel[100],par[100];
1195     char *filtmp;
1196     Double_t tmp1,tmp2;
1197     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1198
1199     res->Thresholds(tmp1,tmp2);
1200     Int_t minval = static_cast<Int_t>(tmp1);
1201
1202     res->Filenames(input,basel,par);
1203     fFileName = par;
1204     //
1205     filtmp = gSystem->ExpandPathName(fFileName.Data());
1206     FILE *param = fopen(filtmp,"r");
1207     na = 0;
1208
1209     if(param) {
1210         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1211             if (pos != na+1) {
1212                 Error("Init2D","Anode number not in increasing order!",filtmp);
1213                 exit(1);
1214             } // end if pos != na+1
1215             savemu[na] = mu;
1216             savesigma[na] = sigma;
1217             if ((2.*sigma) < mu) {
1218                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1219                 mu = 2.0 * sigma;
1220             } else fD[na] = 0;
1221             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1222             if (tempTh < 0) tempTh=0;
1223             fT1[na] = tempTh;
1224             tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
1225             if (tempTh < 0) tempTh=0;
1226             fT2[na] = tempTh;
1227             na++;
1228         } // end while
1229     } else {
1230         Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1231         exit(1);
1232     } // end if(param)
1233
1234     fclose(param);
1235     delete [] filtmp;
1236     delete [] savemu;
1237     delete [] savesigma;
1238 }
1239 //______________________________________________________________________
1240 void AliITSsimulationSDD::Compress2D(){
1241     // simple ITS cluster finder -- online zero-suppression conditions
1242
1243     Int_t db,tl,th; 
1244     Double_t tmp1,tmp2;
1245     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1246
1247     res->Thresholds(tmp1,tmp2); 
1248     Int_t minval   = static_cast<Int_t>(tmp1);
1249     Bool_t write   = res->OutputOption();   
1250     Bool_t do10to8 = res->Do10to8();
1251     Int_t nz, nl, nh, low, i, j; 
1252
1253     for (i=0; i<fNofMaps; i++) {
1254         CompressionParam(i,db,tl,th);
1255         nz  = 0; 
1256         nl  = 0;
1257         nh  = 0;
1258         low = 0;
1259         for (j=0; j<fMaxNofSamples; j++) {
1260             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1261             signal -= db; // if baseline eq. is done here
1262             if (signal <= 0) {nz++; continue;}
1263             if ((signal - tl) < minval) low++;
1264             if ((signal - th) >= minval) {
1265                 nh++;
1266                 Bool_t cond=kTRUE;
1267                 FindCluster(i,j,signal,minval,cond);
1268                 if(cond && j &&
1269                    ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1270                     if(do10to8) signal = Convert10to8(signal);
1271                     AddDigit(i,j,signal);
1272                 } // end if cond&&j&&()
1273             } else if ((signal - tl) >= minval) nl++;
1274         } // end for j loop time samples
1275         if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1276     } //end for i loop anodes
1277
1278     char hname[30];
1279     if (write) {
1280         sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1281         TreeB()->Write(hname);
1282         // reset tree
1283         TreeB()->Reset();
1284     } // end if write
1285 }
1286 //______________________________________________________________________
1287 void  AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
1288                                        Int_t minval,Bool_t &cond){
1289     // Find clusters according to the online 2D zero-suppression algorithm
1290     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1291     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1292   
1293     Bool_t do10to8 = res->Do10to8();
1294     Bool_t high    = kFALSE;
1295
1296     fHitMap2->FlagHit(i,j);
1297     //
1298     //  check the online zero-suppression conditions
1299     //  
1300     const Int_t kMaxNeighbours = 4;
1301     Int_t nn;
1302     Int_t dbx,tlx,thx;  
1303     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1304     seg->Neighbours(i,j,&nn,xList,yList);
1305     Int_t in,ix,iy,qns;
1306     for (in=0; in<nn; in++) {
1307         ix=xList[in];
1308         iy=yList[in];
1309         if (fHitMap2->TestHit(ix,iy)==kUnused) {
1310             CompressionParam(ix,dbx,tlx,thx);
1311             Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1312             qn -= dbx; // if baseline eq. is done here
1313             if ((qn-tlx) < minval) {
1314                 fHitMap2->FlagHit(ix,iy);
1315                 continue;
1316             } else {
1317                 if ((qn - thx) >= minval) high=kTRUE;
1318                 if (cond) {
1319                     if(do10to8) signal = Convert10to8(signal);
1320                     AddDigit(i,j,signal);
1321                 } // end if cond
1322                 if(do10to8) qns = Convert10to8(qn);
1323                 else qns=qn;
1324                 if (!high) AddDigit(ix,iy,qns);
1325                 cond=kFALSE;
1326                 if(!high) fHitMap2->FlagHit(ix,iy);
1327             } // end if qn-tlx < minval
1328         } // end if  TestHit
1329     } // end for in loop over neighbours
1330 }
1331 //______________________________________________________________________
1332 void AliITSsimulationSDD::Init1D(){
1333     // this is just a copy-paste of input taken from 2D algo
1334     // Torino people should give input
1335     // Read 1D zero-suppression parameters for SDD
1336     
1337     if (!strstr(fParam.Data(),"file")) return;
1338
1339     Int_t na,pos,tempTh;
1340     Float_t mu,sigma;
1341     Float_t *savemu    = new Float_t [fNofMaps];
1342     Float_t *savesigma = new Float_t [fNofMaps];
1343     char input[100],basel[100],par[100];
1344     char *filtmp;
1345     Double_t tmp1,tmp2;
1346     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1347
1348     res->Thresholds(tmp1,tmp2);
1349     Int_t minval = static_cast<Int_t>(tmp1);
1350
1351     res->Filenames(input,basel,par);
1352     fFileName=par;
1353
1354     //  set first the disable and tol param
1355     SetCompressParam();
1356     //
1357     filtmp = gSystem->ExpandPathName(fFileName.Data());
1358     FILE *param = fopen(filtmp,"r");
1359     na = 0;
1360  
1361     if (param) {
1362         fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1363         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1364             if (pos != na+1) {
1365                 Error("Init1D","Anode number not in increasing order!",filtmp);
1366                 exit(1);
1367             } // end if pos != na+1
1368             savemu[na]=mu;
1369             savesigma[na]=sigma;
1370             if ((2.*sigma) < mu) {
1371                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1372                 mu = 2.0 * sigma;
1373             } else fD[na] = 0;
1374             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1375             if (tempTh < 0) tempTh=0;
1376             fT1[na] = tempTh;
1377             na++;
1378         } // end while
1379     } else {
1380         Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1381         exit(1);
1382     } // end if(param)
1383
1384     fclose(param);
1385     delete [] filtmp;
1386     delete [] savemu;
1387     delete [] savesigma;
1388
1389 //______________________________________________________________________
1390 void AliITSsimulationSDD::Compress1D(){
1391     // 1D zero-suppression algorithm (from Gianluca A.)
1392     Int_t    dis,tol,thres,decr,diff;
1393     UChar_t *str=fStream->Stream();
1394     Int_t    counter=0;
1395     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1396
1397     Bool_t   do10to8=res->Do10to8();
1398     Int_t    last=0;
1399     Int_t    k,i,j;
1400
1401     for (k=0; k<2; k++) {
1402         tol = Tolerance(k);
1403         dis = Disable(k);  
1404         for (i=0; i<fNofMaps/2; i++) {
1405             Bool_t firstSignal=kTRUE;
1406             Int_t idx=i+k*fNofMaps/2;
1407             if( !fAnodeFire[idx] ) continue;
1408             CompressionParam(idx,decr,thres); 
1409             
1410             for (j=0; j<fMaxNofSamples; j++) {
1411                 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1412                 signal -= decr;  // if baseline eq.
1413                 if(do10to8) signal = Convert10to8(signal);
1414                 if (signal <= thres) {
1415                     signal=0;
1416                     diff=128; 
1417                     last=0; 
1418                     // write diff in the buffer for HuffT
1419                     str[counter]=(UChar_t)diff;
1420                     counter++;
1421                     continue;
1422                 } // end if signal <= thres
1423                 diff=signal-last;
1424                 if (diff > 127) diff=127;
1425                 if (diff < -128) diff=-128;
1426                 if (signal < dis) {
1427                     // tol has changed to 8 possible cases ? - one can write
1428                     // this if(TMath::Abs(diff)<tol) ... else ...
1429                     if(TMath::Abs(diff)<tol) diff=0;
1430                     // or keep it as it was before
1431                     AddDigit(idx,j,last+diff);
1432                 } else {
1433                     AddDigit(idx,j,signal);
1434                 } // end if singal < dis
1435                 diff += 128;
1436                 // write diff in the buffer used to compute Huffman tables
1437                 if (firstSignal) str[counter]=(UChar_t)signal;
1438                 else str[counter]=(UChar_t)diff;
1439                 counter++;
1440                 last=signal;
1441                 firstSignal=kFALSE;
1442             } // end for j loop time samples
1443         } // end for i loop anodes  one half of detector 
1444     } //  end for k
1445
1446     // check
1447     fStream->CheckCount(counter);
1448
1449     // open file and write out the stream of diff's
1450     static Bool_t open=kTRUE;
1451     static TFile *outFile;
1452     Bool_t write = res->OutputOption();
1453     TDirectory *savedir = gDirectory;
1454  
1455     if (write ) {
1456         if(open) {
1457             SetFileName("stream.root");
1458             cout<<"filename "<<fFileName<<endl;
1459             outFile=new TFile(fFileName,"recreate");
1460             cout<<"I have opened "<<fFileName<<" file "<<endl;
1461         } // end if open
1462         open = kFALSE;
1463         outFile->cd();
1464         fStream->Write();
1465     }  // endif write
1466
1467     fStream->ClearStream();
1468
1469     // back to galice.root file
1470     if(savedir) savedir->cd();
1471 }
1472 //______________________________________________________________________
1473 void AliITSsimulationSDD::StoreAllDigits(){
1474     // if non-zero-suppressed data
1475     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1476
1477     Bool_t do10to8 = res->Do10to8();
1478     Int_t i, j, digits[3];
1479
1480     for (i=0; i<fNofMaps; i++) {
1481         for (j=0; j<fMaxNofSamples; j++) {
1482             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1483             if(do10to8) signal = Convert10to8(signal);
1484             digits[0] = i;
1485             digits[1] = j;
1486             digits[2] = signal;
1487             fITS->AddRealDigit(1,digits);
1488         } // end for j
1489     } // end for i
1490
1491 //______________________________________________________________________
1492 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
1493     // Creates histograms of maps for debugging
1494     Int_t i;
1495
1496     fHis=new TObjArray(fNofMaps);
1497     for (i=0;i<fNofMaps;i++) {
1498         TString sddName("sdd_");
1499         Char_t candNum[4];
1500         sprintf(candNum,"%d",i+1);
1501         sddName.Append(candNum);
1502         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1503                              0.,(Float_t) scale*fMaxNofSamples), i);
1504     } // end for i
1505 }
1506 //______________________________________________________________________
1507 void AliITSsimulationSDD::FillHistograms(){
1508     // fill 1D histograms from map
1509
1510     if (!fHis) return;
1511
1512     for( Int_t i=0; i<fNofMaps; i++) {
1513         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1514         Int_t nsamples = hist->GetNbinsX();
1515         for( Int_t j=0; j<nsamples; j++) {
1516             Double_t signal=fHitMap2->GetSignal(i,j);
1517             hist->Fill((Float_t)j,signal);
1518         } // end for j
1519     } // end for i
1520 }
1521 //______________________________________________________________________
1522 void AliITSsimulationSDD::ResetHistograms(){
1523     // Reset histograms for this detector
1524     Int_t i;
1525
1526     for (i=0;i<fNofMaps;i++ ) {
1527         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
1528     } // end for i
1529 }
1530 //______________________________________________________________________
1531 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) { 
1532     // Fills a histogram from a give anode.  
1533
1534     if (!fHis) return 0;
1535
1536     if(wing <=0 || wing > 2) {
1537         Warning("GetAnode","Wrong wing number: %d",wing);
1538         return NULL;
1539     } // end if wing <=0 || wing >2
1540     if(anode <=0 || anode > fNofMaps/2) {
1541         Warning("GetAnode","Wrong anode number: %d",anode);
1542         return NULL;
1543     } // end if ampde <=0 || andoe > fNofMaps/2
1544
1545     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1546     return (TH1F*)(fHis->At(index));
1547 }
1548 //______________________________________________________________________
1549 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
1550     // Writes the histograms to a file
1551
1552     if (!fHis) return;
1553
1554     hfile->cd();
1555     Int_t i;
1556     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
1557     return;
1558 }
1559 //______________________________________________________________________
1560 Float_t AliITSsimulationSDD::GetNoise() {  
1561     // Returns the noise value
1562     //Bool_t do10to8=GetResp()->Do10to8();
1563     //noise will always be in the liniar part of the signal
1564     Int_t decr;
1565     Int_t threshold = fT1[0];
1566     char opt1[20], opt2[20];
1567     AliITSresponseSDD* res = (AliITSresponseSDD*)GetResponseModel(fDetType->GetITSgeom()->GetStartSDD());
1568
1569     res->ParamOptions(opt1,opt2);
1570     fParam=opt2;
1571     char *same = strstr(opt1,"same");
1572     Double_t noise,baseline;
1573     if (same) {
1574         res->GetNoiseParam(noise,baseline);
1575     } else {
1576         static Bool_t readfile=kTRUE;
1577         //read baseline and noise from file
1578         if (readfile) ReadBaseline();
1579         readfile=kFALSE;
1580     } // end if same
1581
1582     TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1583     if(c2) delete c2->GetPrimitive("noisehist");
1584     if(c2) delete c2->GetPrimitive("anode");
1585     else     c2=new TCanvas("c2");
1586     c2->cd();
1587     c2->SetFillColor(0);
1588
1589     TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1590     TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1591                            (float)fMaxNofSamples);
1592     Int_t i,k;
1593     for (i=0;i<fNofMaps;i++) {
1594         CompressionParam(i,decr,threshold); 
1595         if  (!same) GetAnodeBaseline(i,baseline,noise);
1596         anode->Reset();
1597         for (k=0;k<fMaxNofSamples;k++) {
1598             Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1599             //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1600             if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1601             anode->Fill((float)k,signal);
1602         } // end for k
1603         anode->Draw();
1604         c2->Update();
1605     } // end for i
1606     TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1607     noisehist->Fit("gnoise","RQ");
1608     noisehist->Draw();
1609     c2->Update();
1610     Float_t mnoise = gnoise->GetParameter(1);
1611     cout << "mnoise : " << mnoise << endl;
1612     Float_t rnoise = gnoise->GetParameter(2);
1613     cout << "rnoise : " << rnoise << endl;
1614     delete noisehist;
1615     return rnoise;
1616 }
1617 //______________________________________________________________________
1618 void AliITSsimulationSDD::WriteSDigits(){
1619     // Fills the Summable digits Tree
1620     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1621
1622     for( Int_t i=0; i<fNofMaps; i++ ) {
1623         if( !fAnodeFire[i] ) continue;
1624         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1625             Double_t sig = fHitMap2->GetSignal( i, j );
1626             if( sig > 0.2 ) {
1627                 Int_t jdx = j*fScaleSize;
1628                 Int_t index = fpList->GetHitIndex( i, j );
1629                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1630                 // put the fScaleSize analog digits in only one
1631                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1632                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1633                     if( pItemTmp == 0 ) continue;
1634                     pItemTmp2.Add( pItemTmp );
1635                 }
1636                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1637                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1638                 aliITS->AddSumDigit( pItemTmp2 );
1639             } // end if (sig > 0.2)
1640         }
1641     }
1642     return;
1643 }
1644 //______________________________________________________________________
1645 void AliITSsimulationSDD::PrintStatus() const {
1646     // Print SDD simulation Parameters
1647
1648     cout << "**************************************************" << endl;
1649     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1650     cout << "**************************************************" << endl;
1651     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1652     cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1653     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1654     cout << "Number pf Anodes used: " << fNofMaps << endl;
1655     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1656     cout << "Scale size factor: " << fScaleSize << endl;
1657     cout << "**************************************************" << endl;
1658 }