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