Cosmetic changes.
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSDD.cxx
CommitLineData
b0f5e3fc 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
b0f5e3fc 16#include <iostream.h>
17#include <stdlib.h>
18#include <stdio.h>
1ca7869b 19#include <string.h>
20
94de3818 21#include <TSystem.h>
22#include <TROOT.h>
608f25d8 23#include <TStopwatch.h>
ece86d9a 24#include <TCanvas.h>
25#include <TF1.h>
26#include <TRandom.h>
1ca7869b 27#include <TH1.h>
28#include <TFile.h>
29#include <TVector.h>
30#include <TArrayI.h>
31#include <TArrayF.h>
ece86d9a 32
b0f5e3fc 33#include "AliRun.h"
e8189707 34#include "AliITS.h"
ece86d9a 35#include "AliITShit.h"
36#include "AliITSdigit.h"
37#include "AliITSmodule.h"
e8189707 38#include "AliITSMapA1.h"
39#include "AliITSMapA2.h"
e8189707 40#include "AliITSetfSDD.h"
41#include "AliITSRawData.h"
b0f5e3fc 42#include "AliITSHuffman.h"
ece86d9a 43#include "AliITSsegmentation.h"
44#include "AliITSresponse.h"
1ca7869b 45#include "AliITSsimulationSDD.h"
b0f5e3fc 46
b0f5e3fc 47ClassImp(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
8a33ae9e 66//______________________________________________________________________
b0f5e3fc 67Int_t power(Int_t b, Int_t e) {
8a33ae9e 68 // compute b to the e power, where both b and e are Int_ts.
69 Int_t power = 1,i;
b0f5e3fc 70
8a33ae9e 71 for(i=0; i<e; i++) power *= b;
72 return power;
73}
74//______________________________________________________________________
b0f5e3fc 75void FastFourierTransform(AliITSetfSDD *alisddetf,Double_t *real,
76 Double_t *imag,Int_t direction) {
8a33ae9e 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
b0f5e3fc 106
8a33ae9e 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
b0f5e3fc 131 return;
132}
8a33ae9e 133//______________________________________________________________________
b0f5e3fc 134AliITSsimulationSDD::AliITSsimulationSDD(){
8a33ae9e 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();
b0f5e3fc 162}
8a33ae9e 163//______________________________________________________________________
164AliITSsimulationSDD::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();
b0f5e3fc 171}
8a33ae9e 172//______________________________________________________________________
173AliITSsimulationSDD& 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 ;
b0f5e3fc 180}
8a33ae9e 181//______________________________________________________________________
182AliITSsimulationSDD::AliITSsimulationSDD(AliITSsegmentation *seg,AliITSresponse *resp){
38867c90 183 // Standard Constructor
ece86d9a 184
38867c90 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();
b0f5e3fc 197
38867c90 198 fHitMap2 = new AliITSMapA2(fSegmentation,fScaleSize,1);
199 fHitMap1 = new AliITSMapA1(fSegmentation);
b0f5e3fc 200
8a33ae9e 201 fNofMaps = fSegmentation->Npz();
202 fMaxNofSamples = fSegmentation->Npx();
203
38867c90 204 Float_t sddLength = fSegmentation->Dx();
8a33ae9e 205 Float_t sddWidth = fSegmentation->Dz();
b0f5e3fc 206
8a33ae9e 207 Int_t dummy = 0;
38867c90 208 Float_t anodePitch = fSegmentation->Dpz(dummy);
8a33ae9e 209 Double_t timeStep = (Double_t)fSegmentation->Dpx(dummy);
210 Float_t driftSpeed = fResponse->DriftSpeed();
b0f5e3fc 211
38867c90 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
b0f5e3fc 217
38867c90 218 if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
219 Error("AliITSsimulationSDD",
220 "Time Interval > Allowed Time Interval: exit\n");
221 return;
222 } // end if
b0f5e3fc 223
38867c90 224 fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
225 fResponse->Electronics());
b0f5e3fc 226
38867c90 227 char opt1[20], opt2[20];
228 fResponse->ParamOptions(opt1,opt2);
8a33ae9e 229 fParam = opt2;
38867c90 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
b0f5e3fc 238
38867c90 239 const char *kopt=fResponse->ZeroSuppOption();
240 if (strstr(fParam,"file") ) {
241 fD.Set(fNofMaps);
242 fT1.Set(fNofMaps);
243 if (strstr(kopt,"2D")) {
b0f5e3fc 244 fT2.Set(fNofMaps);
245 fTol.Set(0);
246 Init2D(); // desactivate if param change module by module
38867c90 247 } else if(strstr(kopt,"1D")) {
b0f5e3fc 248 fT2.Set(2);
249 fTol.Set(2);
250 Init1D(); // desactivate if param change module by module
38867c90 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
b0f5e3fc 259
8a33ae9e 260 Bool_t write = fResponse->OutputOption();
38867c90 261 if(write && strstr(kopt,"2D")) MakeTreeB();
b0f5e3fc 262
38867c90 263 // call here if baseline does not change by module
264 // ReadBaseline();
b0f5e3fc 265
8a33ae9e 266 fITS = (AliITS*)gAlice->GetModule("ITS");
267 Int_t size = fNofMaps*fMaxNofSamples;
268 fStream = new AliITSInStream(size);
b0f5e3fc 269
38867c90 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];
5d18fa90 274
b0f5e3fc 275}
8a33ae9e 276//______________________________________________________________________
b0f5e3fc 277AliITSsimulationSDD::~AliITSsimulationSDD() {
8a33ae9e 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;
b0f5e3fc 302}
8a33ae9e 303//______________________________________________________________________
b0f5e3fc 304void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
8a33ae9e 305 // create maps to build the lists of tracks for each digit
b0f5e3fc 306
307 TObjArray *fHits = mod->GetHits();
8a33ae9e 308 Int_t nhits = fHits->GetEntriesFast();
309 fModule = md;
310 fEvent = ev;
311
ece86d9a 312 if (!nhits && fCheckNoise) {
313 ChargeToSignal();
314 GetNoise();
315 fHitMap2->ClearMap();
316 return;
317 } else if (!nhits) return;
318
8a33ae9e 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(); //
44a312c3 341
b0f5e3fc 342 // Piergiorgio's part (apart for few variables which I made float
343 // when i thought that can be done
b0f5e3fc 344 // Fill detector maps with GEANT hits
345 // loop over hits in the module
346
8a33ae9e 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.
44a312c3 387
b0f5e3fc 388 for(ii=0; ii<nhits; ii++) {
8a33ae9e 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) :
b0f5e3fc 468
8a33ae9e 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"
44a312c3 546 } // end loop over hits
b0f5e3fc 547
8a33ae9e 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();
b0f5e3fc 562}
8a33ae9e 563//______________________________________________________________________
b0f5e3fc 564void AliITSsimulationSDD::ListOfFiredCells(Int_t *arg,Double_t timeAmplitude,
8a33ae9e 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
b0f5e3fc 634
635#ifdef print
8a33ae9e 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);
b0f5e3fc 649 // set printings
8a33ae9e 650 } // end for tr
651 } // end if nptracks
b0f5e3fc 652#endif
8a33ae9e 653 } // end if pdigit
b0f5e3fc 654
8a33ae9e 655 // update counter and countadr for next call.
656 arg[4] = counter;
657 arg[5] = countadr;
b0f5e3fc 658}
b0f5e3fc 659//____________________________________________
660
661void AliITSsimulationSDD::AddDigit(Int_t i, Int_t j, Int_t signal){
8a33ae9e 662 // Adds a Digit.
b0f5e3fc 663 // tag with -1 signals coming from background tracks
664 // tag with -2 signals coming from pure electronic noise
665
e8189707 666 Int_t digits[3], tracks[3], hits[3];
b0f5e3fc 667 Float_t phys, charges[3];
668
e8189707 669 Int_t trk[20], htrk[20];
b0f5e3fc 670 Float_t chtrk[20];
671
ece86d9a 672 Bool_t do10to8=fResponse->Do10to8();
673
674 if(do10to8) signal=Convert8to10(signal);
b0f5e3fc 675 AliITSTransientDigit *obj = (AliITSTransientDigit*)fHitMap1->GetHit(i,j);
8a33ae9e 676 digits[0] = i;
677 digits[1] = j;
678 digits[2] = signal;
b0f5e3fc 679 if (!obj) {
680 phys=0;
681 Int_t k;
e8189707 682 for (k=0;k<3;k++) {
8a33ae9e 683 tracks[k]=-2;
684 charges[k]=0;
685 hits[k]=-1;
686 } // end for k
e8189707 687 fITS->AddSimDigit(1,phys,digits,tracks,hits,charges);
b0f5e3fc 688 } else {
8a33ae9e 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);
b0f5e3fc 727
8a33ae9e 728 } // end if/else !obj
b0f5e3fc 729}
8a33ae9e 730//______________________________________________________________________
731void 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
b0f5e3fc 777}
8a33ae9e 778//______________________________________________________________________
b0f5e3fc 779void AliITSsimulationSDD::ChargeToSignal() {
8a33ae9e 780 // add baseline, noise, electronics and ADC saturation effects
b0f5e3fc 781
8a33ae9e 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
ece86d9a 821
ece86d9a 822 for (i=0;i<fNofMaps;i++) {
8a33ae9e 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);
ece86d9a 831 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
8a33ae9e 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);
ece86d9a 838 for(k=0; k<fMaxNofSamples; k++) {
8a33ae9e 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
ece86d9a 850 return;
b0f5e3fc 851}
8a33ae9e 852//______________________________________________________________________
b0f5e3fc 853void AliITSsimulationSDD::GetAnodeBaseline(Int_t i,Float_t &baseline,
854 Float_t &noise){
8a33ae9e 855 // Returns the Baseline for a particular anode.
856 baseline = fBaseline[i];
857 noise = fNoise[i];
b0f5e3fc 858}
8a33ae9e 859//______________________________________________________________________
b0f5e3fc 860void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
861 Int_t &th){
8a33ae9e 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
b0f5e3fc 873}
8a33ae9e 874//______________________________________________________________________
b0f5e3fc 875void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl){
8a33ae9e 876 // returns the compression alogirthm parameters
877 Int_t size = fD.GetSize();
b0f5e3fc 878
8a33ae9e 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
b0f5e3fc 888}
8a33ae9e 889//______________________________________________________________________
b0f5e3fc 890void AliITSsimulationSDD::SetCompressParam(){
8a33ae9e 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
b0f5e3fc 903}
8a33ae9e 904//______________________________________________________________________
b0f5e3fc 905void AliITSsimulationSDD::ReadBaseline(){
8a33ae9e 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
b0f5e3fc 911
912 Int_t na,pos;
913 Float_t bl,n;
e8189707 914 char input[100], base[100], param[100];
b0f5e3fc 915 char *filtmp;
916
e8189707 917 fResponse->Filenames(input,base,param);
b0f5e3fc 918 fFileName=base;
919//
920 filtmp = gSystem->ExpandPathName(fFileName.Data());
921 FILE *bline = fopen(filtmp,"r");
8a33ae9e 922// printf("filtmp %s\n",filtmp);
b0f5e3fc 923 na = 0;
924
925 if(bline) {
8a33ae9e 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
b0f5e3fc 936 } else {
8a33ae9e 937 Error("ReadBaseline"," THE BASELINE FILE %s DOES NOT EXIST !",filtmp);
938 exit(1);
b0f5e3fc 939 } // end if(bline)
8a33ae9e 940
b0f5e3fc 941 fclose(bline);
942 delete [] filtmp;
b0f5e3fc 943}
8a33ae9e 944//______________________________________________________________________
945Int_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;
b0f5e3fc 954}
8a33ae9e 955//______________________________________________________________________
956Int_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//______________________________________________________________________
b0f5e3fc 978AliITSMap* AliITSsimulationSDD::HitMap(Int_t i){
8a33ae9e 979 //Return the correct map.
980
b0f5e3fc 981 return ((i==0)? fHitMap1 : fHitMap2);
982}
8a33ae9e 983//______________________________________________________________________
e8189707 984void AliITSsimulationSDD::ZeroSuppression(const char *option) {
8a33ae9e 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();
b0f5e3fc 994}
8a33ae9e 995//______________________________________________________________________
b0f5e3fc 996void AliITSsimulationSDD::Init2D(){
8a33ae9e 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
b0f5e3fc 1005
1006 if (!strstr(fParam,"file")) return;
1007
1008 Int_t na,pos,tempTh;
1009 Float_t mu,sigma;
8a33ae9e 1010 Float_t *savemu = new Float_t [fNofMaps];
e8189707 1011 Float_t *savesigma = new Float_t [fNofMaps];
1012 char input[100],basel[100],par[100];
b0f5e3fc 1013 char *filtmp;
b0f5e3fc 1014 Int_t minval = fResponse->MinVal();
1015
e8189707 1016 fResponse->Filenames(input,basel,par);
8a33ae9e 1017 fFileName = par;
b0f5e3fc 1018//
1019 filtmp = gSystem->ExpandPathName(fFileName.Data());
1020 FILE *param = fopen(filtmp,"r");
1021 na = 0;
1022
1023 if(param) {
8a33ae9e 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;
b0f5e3fc 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++;
8a33ae9e 1042 } // end while
b0f5e3fc 1043 } else {
8a33ae9e 1044 Error("Init2D "," THE FILE %s DOES NOT EXIST !",filtmp);
1045 exit(1);
b0f5e3fc 1046 } // end if(param)
8a33ae9e 1047
b0f5e3fc 1048 fclose(param);
1049 delete [] filtmp;
5d18fa90 1050 delete [] savemu;
e8189707 1051 delete [] savesigma;
8a33ae9e 1052}
1053//______________________________________________________________________
b0f5e3fc 1054void AliITSsimulationSDD::Compress2D(){
8a33ae9e 1055 // simple ITS cluster finder -- online zero-suppression conditions
b0f5e3fc 1056
b0f5e3fc 1057 Int_t db,tl,th;
8a33ae9e 1058 Int_t minval = fResponse->MinVal();
1059 Bool_t write = fResponse->OutputOption();
1060 Bool_t do10to8 = fResponse->Do10to8();
b0f5e3fc 1061 Int_t nz, nl, nh, low, i, j;
1062
e8189707 1063 for (i=0; i<fNofMaps; i++) {
b0f5e3fc 1064 CompressionParam(i,db,tl,th);
8a33ae9e 1065 nz = 0;
1066 nl = 0;
1067 nh = 0;
1068 low = 0;
e8189707 1069 for (j=0; j<fMaxNofSamples; j++) {
b0f5e3fc 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;
b0f5e3fc 1077 FindCluster(i,j,signal,minval,cond);
8a33ae9e 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&&()
b0f5e3fc 1082 } else if ((signal - tl) >= minval) nl++;
8a33ae9e 1083 } // end for j loop time samples
1084 if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1085 } //end for i loop anodes
b0f5e3fc 1086
8a33ae9e 1087 char hname[30];
1088 if (write) {
b0f5e3fc 1089 sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1090 TreeB()->Write(hname);
1091 // reset tree
1092 TreeB()->Reset();
8a33ae9e 1093 } // end if write
1094}
1095//______________________________________________________________________
b0f5e3fc 1096void AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
ece86d9a 1097 Int_t minval,Bool_t &cond){
8a33ae9e 1098 // Find clusters according to the online 2D zero-suppression algorithm
1099 Bool_t do10to8 = fResponse->Do10to8();
1100 Bool_t high = kFALSE;
b0f5e3fc 1101
1102 fHitMap2->FlagHit(i,j);
1103//
1104// check the online zero-suppression conditions
1105//
8a33ae9e 1106 const Int_t kMaxNeighbours = 4;
b0f5e3fc 1107 Int_t nn;
1108 Int_t dbx,tlx,thx;
8a33ae9e 1109 Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
e8189707 1110 fSegmentation->Neighbours(i,j,&nn,xList,yList);
ece86d9a 1111 Int_t in,ix,iy,qns;
e8189707 1112 for (in=0; in<nn; in++) {
1113 ix=xList[in];
1114 iy=yList[in];
b0f5e3fc 1115 if (fHitMap2->TestHit(ix,iy)==kUnused) {
8a33ae9e 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
b0f5e3fc 1136}
8a33ae9e 1137//______________________________________________________________________
b0f5e3fc 1138void AliITSsimulationSDD::Init1D(){
8a33ae9e 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
b0f5e3fc 1142
1143 if (!strstr(fParam,"file")) return;
1144
1145 Int_t na,pos,tempTh;
1146 Float_t mu,sigma;
8a33ae9e 1147 Float_t *savemu = new Float_t [fNofMaps];
e8189707 1148 Float_t *savesigma = new Float_t [fNofMaps];
1149 char input[100],basel[100],par[100];
b0f5e3fc 1150 char *filtmp;
b0f5e3fc 1151 Int_t minval = fResponse->MinVal();
8a33ae9e 1152
e8189707 1153 fResponse->Filenames(input,basel,par);
1154 fFileName=par;
b0f5e3fc 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) {
8a33ae9e 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
b0f5e3fc 1181 } else {
8a33ae9e 1182 Error("Init1D "," THE FILE %s DOES NOT EXIST !",filtmp);
1183 exit(1);
b0f5e3fc 1184 } // end if(param)
8a33ae9e 1185
b0f5e3fc 1186 fclose(param);
1187 delete [] filtmp;
749bd21a 1188 delete [] savemu;
e8189707 1189 delete [] savesigma;
8a33ae9e 1190}
1191//______________________________________________________________________
b0f5e3fc 1192void AliITSsimulationSDD::Compress1D(){
1193 // 1D zero-suppression algorithm (from Gianluca A.)
8a33ae9e 1194 Int_t dis,tol,thres,decr,diff;
b0f5e3fc 1195 UChar_t *str=fStream->Stream();
8a33ae9e 1196 Int_t counter=0;
1197 Bool_t do10to8=fResponse->Do10to8();
1198 Int_t last=0;
1199 Int_t k,i,j;
ece86d9a 1200
ece86d9a 1201 for (k=0; k<2; k++) {
8a33ae9e 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 ...
ece86d9a 1227 if(TMath::Abs(diff)<tol) diff=0;
1228 // or keep it as it was before
1229 /*
b0f5e3fc 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;
ece86d9a 1233 */
1234 AddDigit(idx,j,last+diff);
8a33ae9e 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
b0f5e3fc 1248
1249 // check
1250 fStream->CheckCount(counter);
1251
1252 // open file and write out the stream of diff's
b0f5e3fc 1253 static Bool_t open=kTRUE;
e8189707 1254 static TFile *outFile;
b0f5e3fc 1255 Bool_t write = fResponse->OutputOption();
1256
1257 if (write ) {
1258 if(open) {
1259 SetFileName("stream.root");
1260 cout<<"filename "<<fFileName<<endl;
e8189707 1261 outFile=new TFile(fFileName,"recreate");
b0f5e3fc 1262 cout<<"I have opened "<<fFileName<<" file "<<endl;
8a33ae9e 1263 } // end if open
1264 open = kFALSE;
e8189707 1265 outFile->cd();
b0f5e3fc 1266 fStream->Write();
1267 } // endif write
1268
8a33ae9e 1269 fStream->ClearStream();
b0f5e3fc 1270
8a33ae9e 1271 // back to galice.root file
b0f5e3fc 1272
8a33ae9e 1273 TTree *fAli=gAlice->TreeK();
1274 TFile *file = 0;
b0f5e3fc 1275
8a33ae9e 1276 if (fAli) file =fAli->GetCurrentFile();
1277 file->cd();
1278}
1279//______________________________________________________________________
b0f5e3fc 1280void AliITSsimulationSDD::StoreAllDigits(){
8a33ae9e 1281 // if non-zero-suppressed data
1282 Bool_t do10to8 = fResponse->Do10to8();
ece86d9a 1283 Int_t i, j, digits[3];
8a33ae9e 1284
e8189707 1285 for (i=0; i<fNofMaps; i++) {
1286 for (j=0; j<fMaxNofSamples; j++) {
8a33ae9e 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
b0f5e3fc 1296}
8a33ae9e 1297//______________________________________________________________________
ece86d9a 1298void AliITSsimulationSDD::CreateHistograms(Int_t scale){
8a33ae9e 1299 // Creates histograms of maps for debugging
1300 Int_t i;
ece86d9a 1301
1302 fHis=new TObjArray(fNofMaps);
e8189707 1303 for (i=0;i<fNofMaps;i++) {
a14dba92 1304 TString sddName("sdd_");
b0f5e3fc 1305 Char_t candNum[4];
1306 sprintf(candNum,"%d",i+1);
ece86d9a 1307 sddName.Append(candNum);
2682e810 1308 //PH (*fHis)[i] = new TH1F(sddName.Data(),"SDD maps",
8a33ae9e 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
b0f5e3fc 1313}
8a33ae9e 1314//______________________________________________________________________
ece86d9a 1315void AliITSsimulationSDD::FillHistograms(){
8a33ae9e 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
ece86d9a 1328}
8a33ae9e 1329//______________________________________________________________________
b0f5e3fc 1330void AliITSsimulationSDD::ResetHistograms(){
b0f5e3fc 1331 // Reset histograms for this detector
b0f5e3fc 1332 Int_t i;
8a33ae9e 1333
e8189707 1334 for (i=0;i<fNofMaps;i++ ) {
8a33ae9e 1335 //PH if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
2682e810 1336 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();
8a33ae9e 1337 } // end for i
b0f5e3fc 1338}
8a33ae9e 1339//______________________________________________________________________
b0f5e3fc 1340TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
8a33ae9e 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));
b0f5e3fc 1357}
8a33ae9e 1358//______________________________________________________________________
b0f5e3fc 1359void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
8a33ae9e 1360 // Writes the histograms to a file
b0f5e3fc 1361
8a33ae9e 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;
b0f5e3fc 1369}
8a33ae9e 1370//______________________________________________________________________
ece86d9a 1371Float_t AliITSsimulationSDD::GetNoise() {
8a33ae9e 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();
ece86d9a 1419 c2->Update();
8a33ae9e 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;
b0f5e3fc 1426}
8a33ae9e 1427//______________________________________________________________________
44a312c3 1428void AliITSsimulationSDD::Print() {
8a33ae9e 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;
44a312c3 1441}