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