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