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