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