]> 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++ corrections
[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.;
753 if( res->IsChipBad(res->GetChip(i)) ){
754 printf("Chip bad mod %d chip %d anode %d\n",mod,res->GetChip(i),i);
755 gain=0.;
756 }
5683bd96 757 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
758 fInZR[k] = fHitMap2->GetSignal(i,k);
759 if(bAddGain) fInZR[k]*=gain;
760 if( bAddNoise ) {
761 contrib = (baseline + noise*gRandom->Gaus());
762 fInZR[k] += contrib;
763 }
764 fInZI[k] = 0.;
765 } // end for k
aacedc3e 766 if(!fDoFFT) {
5683bd96 767 for(k=0; k<fMaxNofSamples; k++) {
768 Double_t newcont = 0.;
769 Double_t maxcont = 0.;
770 for(kk=0;kk<fScaleSize;kk++) {
771 newcont = fInZR[fScaleSize*k+kk];
772 if(newcont > maxcont) maxcont = newcont;
773 } // end for kk
774 newcont = maxcont;
775 if (newcont >= maxadc) newcont = maxadc -1;
776 if(newcont >= baseline){
777 Warning("","newcont=%d>=baseline=%d",newcont,baseline);
778 } // end if
779 // back to analog: ?
780 fHitMap2->SetHit(i,k,newcont);
781 } // end for k
782 }else{
783 FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
784 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
785 Double_t rw = fElectronics->GetTraFunReal(k);
786 Double_t iw = fElectronics->GetTraFunImag(k);
787 fOutZR[k] = fInZR[k]*rw - fInZI[k]*iw;
788 fOutZI[k] = fInZR[k]*iw + fInZI[k]*rw;
789 } // end for k
790 FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
791 for(k=0; k<fMaxNofSamples; k++) {
792 Double_t newcont1 = 0.;
793 Double_t maxcont1 = 0.;
794 for(kk=0;kk<fScaleSize;kk++) {
795 newcont1 = fOutZR[fScaleSize*k+kk];
796 if(newcont1 > maxcont1) maxcont1 = newcont1;
797 } // end for kk
798 newcont1 = maxcont1;
799 if (newcont1 >= maxadc) newcont1 = maxadc -1;
800 fHitMap2->SetHit(i,k,newcont1);
801 } // end for k
802 }
803 } // end for i loop over anodes
804 return;
50d05d7b 805}
5683bd96 806
50d05d7b 807//______________________________________________________________________
8ba39da9 808void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
aacedc3e 809 // function add the crosstalk effect to signal
810 // temporal function, should be checked...!!!
8ba39da9 811 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
812
813 Int_t fNofMaps = seg->Npz();
814 Int_t fMaxNofSamples = seg->Npx();
aacedc3e 815
816 // create and inizialice crosstalk map
817 Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
818 if( ctk == NULL ) {
819 Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
820 return;
821 }
822 memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
f45f6658 823 AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
aacedc3e 824 for( Int_t z=0; z<fNofMaps; z++ ) {
f45f6658 825 Double_t baseline = calibr->GetBaseline(z);
aacedc3e 826 Bool_t on = kFALSE;
827 Int_t tstart = 0;
828 Int_t tstop = 0;
829 Int_t nTsteps = 0;
50d05d7b 830
aacedc3e 831 for( Int_t l=0; l<fMaxNofSamples; l++ ) {
832 Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
833 if( fadc > baseline ) {
834 if( on == kFALSE && l<fMaxNofSamples-4 ) {
835 Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
836 if( fadc1 < fadc ) continue;
837 on = kTRUE;
838 nTsteps = 0;
839 tstart = l;
840 }
841 nTsteps++;
842 }
843 else { // end fadc > baseline
844 if( on == kTRUE ) {
845 if( nTsteps > 2 ) {
846 tstop = l;
847 // make smooth derivative
848 Float_t* dev = new Float_t[fMaxNofSamples+1];
849 memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
850 if( ctk == NULL ) {
851 Error( "ApplyCrosstalk",
852 "no memory for temporal array: exit \n" );
853 return;
854 }
855 for( Int_t i=tstart; i<tstop; i++ ) {
856 if( i > 2 && i < fMaxNofSamples-2 )
857 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
858 -0.1*fHitMap2->GetSignal( z,i-1 )
859 +0.1*fHitMap2->GetSignal( z,i+1 )
860 +0.2*fHitMap2->GetSignal( z,i+2 );
861 }
50d05d7b 862
aacedc3e 863 // add crosstalk contribution to neibourg anodes
864 for( Int_t i=tstart; i<tstop; i++ ) {
865 Int_t anode = z - 1;
866 Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
867 Float_t ctktmp = -dev[i1] * 0.25;
868 if( anode > 0 ) {
869 ctk[anode*fMaxNofSamples+i] += ctktmp;
870 }
871 anode = z + 1;
872 if( anode < fNofMaps ) {
873 ctk[anode*fMaxNofSamples+i] += ctktmp;
874 }
875 }
876 delete [] dev;
50d05d7b 877
aacedc3e 878 } // if( nTsteps > 2 )
879 on = kFALSE;
880 } // if( on == kTRUE )
881 } // else
882 }
3d2c9d72 883 }
50d05d7b 884
aacedc3e 885 for( Int_t a=0; a<fNofMaps; a++ )
886 for( Int_t t=0; t<fMaxNofSamples; t++ ) {
887 Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
888 fHitMap2->SetHit( a, t, signal );
889 }
890
891 delete [] ctk;
50d05d7b 892}
f45f6658 893
8a33ae9e 894//______________________________________________________________________
b0f5e3fc 895void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
f45f6658 896 Int_t &th) const{
aacedc3e 897 // Returns the compression alogirthm parameters
f45f6658 898 db=fD[i];
899 tl=fT1[i];
900 th=fT2[i];
b0f5e3fc 901}
8a33ae9e 902//______________________________________________________________________
f45f6658 903void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{
aacedc3e 904 // returns the compression alogirthm parameters
aacedc3e 905
f45f6658 906 db=fD[i];
907 tl=fT1[i];
908
b0f5e3fc 909}
8a33ae9e 910//______________________________________________________________________
f45f6658 911void AliITSsimulationSDD::SetCompressParam(){
aacedc3e 912 // Sets the compression alogirthm parameters
f45f6658 913 AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
914 for(Int_t ian = 0; ian<fNofMaps;ian++){
915 fD[ian] = (Int_t)(calibr->GetBaseline(ian));
916 fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5);
917 fT2[ian] = 0; // used by 2D clustering - not defined yet
918 fTol[ian] = 0; // used by 2D clustering - not defined yet
919 }
b0f5e3fc 920}
8a33ae9e 921//______________________________________________________________________
922Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
aacedc3e 923 // To the 10 to 8 bit lossive compression.
924 // code from Davide C. and Albert W.
925
926 if (signal < 128) return signal;
927 if (signal < 256) return (128+((signal-128)>>1));
928 if (signal < 512) return (192+((signal-256)>>3));
929 if (signal < 1024) return (224+((signal-512)>>4));
930 return 0;
b0f5e3fc 931}
50d05d7b 932/*
8a33ae9e 933//______________________________________________________________________
b0f5e3fc 934AliITSMap* AliITSsimulationSDD::HitMap(Int_t i){
aacedc3e 935 //Return the correct map.
8a33ae9e 936
aacedc3e 937 return ((i==0)? fHitMap1 : fHitMap2);
3d2c9d72 938}
939*/
8a33ae9e 940//______________________________________________________________________
e8189707 941void AliITSsimulationSDD::ZeroSuppression(const char *option) {
aacedc3e 942 // perform the zero suppresion
aacedc3e 943 if (strstr(option,"2D")) {
944 //Init2D(); // activate if param change module by module
945 Compress2D();
946 } else if (strstr(option,"1D")) {
947 //Init1D(); // activate if param change module by module
948 Compress1D();
949 } else StoreAllDigits();
b0f5e3fc 950}
8a33ae9e 951//______________________________________________________________________
b0f5e3fc 952void AliITSsimulationSDD::Init2D(){
aacedc3e 953 // read in and prepare arrays: fD, fT1, fT2
954 // savemu[nanodes], savesigma[nanodes]
955 // read baseline and noise from file - either a .root file and in this
956 // case data should be organised in a tree with one entry for each
957 // module => reading should be done accordingly
958 // or a classic file and do smth. like this ( code from Davide C. and
959 // Albert W.) :
960 // Read 2D zero-suppression parameters for SDD
961
962 if (!strstr(fParam.Data(),"file")) return;
963
964 Int_t na,pos,tempTh;
965 Float_t mu,sigma;
966 Float_t *savemu = new Float_t [fNofMaps];
967 Float_t *savesigma = new Float_t [fNofMaps];
968 char input[100],basel[100],par[100];
969 char *filtmp;
970 Double_t tmp1,tmp2;
f45f6658 971 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 972
973 res->Thresholds(tmp1,tmp2);
aacedc3e 974 Int_t minval = static_cast<Int_t>(tmp1);
975
8ba39da9 976 res->Filenames(input,basel,par);
aacedc3e 977 fFileName = par;
978 //
979 filtmp = gSystem->ExpandPathName(fFileName.Data());
980 FILE *param = fopen(filtmp,"r");
981 na = 0;
982
983 if(param) {
984 while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
985 if (pos != na+1) {
986 Error("Init2D","Anode number not in increasing order!",filtmp);
987 exit(1);
988 } // end if pos != na+1
989 savemu[na] = mu;
990 savesigma[na] = sigma;
991 if ((2.*sigma) < mu) {
992 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
993 mu = 2.0 * sigma;
994 } else fD[na] = 0;
995 tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
996 if (tempTh < 0) tempTh=0;
997 fT1[na] = tempTh;
998 tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
999 if (tempTh < 0) tempTh=0;
1000 fT2[na] = tempTh;
1001 na++;
1002 } // end while
1003 } else {
1004 Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1005 exit(1);
1006 } // end if(param)
1007
1008 fclose(param);
1009 delete [] filtmp;
1010 delete [] savemu;
1011 delete [] savesigma;
8a33ae9e 1012}
1013//______________________________________________________________________
b0f5e3fc 1014void AliITSsimulationSDD::Compress2D(){
aacedc3e 1015 // simple ITS cluster finder -- online zero-suppression conditions
1016
1017 Int_t db,tl,th;
1018 Double_t tmp1,tmp2;
f45f6658 1019 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 1020
1021 res->Thresholds(tmp1,tmp2);
aacedc3e 1022 Int_t minval = static_cast<Int_t>(tmp1);
8ba39da9 1023 Bool_t write = res->OutputOption();
1024 Bool_t do10to8 = res->Do10to8();
aacedc3e 1025 Int_t nz, nl, nh, low, i, j;
f45f6658 1026 SetCompressParam();
aacedc3e 1027 for (i=0; i<fNofMaps; i++) {
1028 CompressionParam(i,db,tl,th);
1029 nz = 0;
1030 nl = 0;
1031 nh = 0;
1032 low = 0;
1033 for (j=0; j<fMaxNofSamples; j++) {
1034 Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1035 signal -= db; // if baseline eq. is done here
1036 if (signal <= 0) {nz++; continue;}
1037 if ((signal - tl) < minval) low++;
1038 if ((signal - th) >= minval) {
1039 nh++;
1040 Bool_t cond=kTRUE;
1041 FindCluster(i,j,signal,minval,cond);
1042 if(cond && j &&
1043 ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1044 if(do10to8) signal = Convert10to8(signal);
1045 AddDigit(i,j,signal);
1046 } // end if cond&&j&&()
1047 } else if ((signal - tl) >= minval) nl++;
1048 } // end for j loop time samples
1049 if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1050 } //end for i loop anodes
1051
1052 char hname[30];
1053 if (write) {
1054 sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1055 TreeB()->Write(hname);
1056 // reset tree
1057 TreeB()->Reset();
1058 } // end if write
8a33ae9e 1059}
1060//______________________________________________________________________
b0f5e3fc 1061void AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
ece86d9a 1062 Int_t minval,Bool_t &cond){
aacedc3e 1063 // Find clusters according to the online 2D zero-suppression algorithm
f45f6658 1064 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 1065 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1066
1067 Bool_t do10to8 = res->Do10to8();
aacedc3e 1068 Bool_t high = kFALSE;
1069
1070 fHitMap2->FlagHit(i,j);
1071 //
1072 // check the online zero-suppression conditions
1073 //
1074 const Int_t kMaxNeighbours = 4;
1075 Int_t nn;
1076 Int_t dbx,tlx,thx;
1077 Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
8ba39da9 1078 seg->Neighbours(i,j,&nn,xList,yList);
aacedc3e 1079 Int_t in,ix,iy,qns;
1080 for (in=0; in<nn; in++) {
1081 ix=xList[in];
1082 iy=yList[in];
1083 if (fHitMap2->TestHit(ix,iy)==kUnused) {
1084 CompressionParam(ix,dbx,tlx,thx);
1085 Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1086 qn -= dbx; // if baseline eq. is done here
1087 if ((qn-tlx) < minval) {
1088 fHitMap2->FlagHit(ix,iy);
1089 continue;
1090 } else {
1091 if ((qn - thx) >= minval) high=kTRUE;
1092 if (cond) {
1093 if(do10to8) signal = Convert10to8(signal);
1094 AddDigit(i,j,signal);
1095 } // end if cond
1096 if(do10to8) qns = Convert10to8(qn);
1097 else qns=qn;
1098 if (!high) AddDigit(ix,iy,qns);
1099 cond=kFALSE;
1100 if(!high) fHitMap2->FlagHit(ix,iy);
1101 } // end if qn-tlx < minval
1102 } // end if TestHit
1103 } // end for in loop over neighbours
b0f5e3fc 1104}
8a33ae9e 1105//______________________________________________________________________
b0f5e3fc 1106void AliITSsimulationSDD::Init1D(){
aacedc3e 1107 // this is just a copy-paste of input taken from 2D algo
1108 // Torino people should give input
1109 // Read 1D zero-suppression parameters for SDD
1110
1111 if (!strstr(fParam.Data(),"file")) return;
1112
1113 Int_t na,pos,tempTh;
1114 Float_t mu,sigma;
1115 Float_t *savemu = new Float_t [fNofMaps];
1116 Float_t *savesigma = new Float_t [fNofMaps];
1117 char input[100],basel[100],par[100];
1118 char *filtmp;
1119 Double_t tmp1,tmp2;
f45f6658 1120 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 1121
1122 res->Thresholds(tmp1,tmp2);
aacedc3e 1123 Int_t minval = static_cast<Int_t>(tmp1);
1124
8ba39da9 1125 res->Filenames(input,basel,par);
aacedc3e 1126 fFileName=par;
1127
1128 // set first the disable and tol param
1129 SetCompressParam();
1130 //
1131 filtmp = gSystem->ExpandPathName(fFileName.Data());
1132 FILE *param = fopen(filtmp,"r");
1133 na = 0;
1134
1135 if (param) {
1136 fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1137 while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1138 if (pos != na+1) {
1139 Error("Init1D","Anode number not in increasing order!",filtmp);
1140 exit(1);
1141 } // end if pos != na+1
1142 savemu[na]=mu;
1143 savesigma[na]=sigma;
1144 if ((2.*sigma) < mu) {
1145 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1146 mu = 2.0 * sigma;
1147 } else fD[na] = 0;
1148 tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1149 if (tempTh < 0) tempTh=0;
1150 fT1[na] = tempTh;
1151 na++;
1152 } // end while
1153 } else {
1154 Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1155 exit(1);
1156 } // end if(param)
1157
1158 fclose(param);
1159 delete [] filtmp;
1160 delete [] savemu;
1161 delete [] savesigma;
8a33ae9e 1162}
1163//______________________________________________________________________
b0f5e3fc 1164void AliITSsimulationSDD::Compress1D(){
aacedc3e 1165 // 1D zero-suppression algorithm (from Gianluca A.)
1166 Int_t dis,tol,thres,decr,diff;
1167 UChar_t *str=fStream->Stream();
1168 Int_t counter=0;
f45f6658 1169 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 1170
1171 Bool_t do10to8=res->Do10to8();
aacedc3e 1172 Int_t last=0;
1173 Int_t k,i,j;
f45f6658 1174 SetCompressParam();
aacedc3e 1175 for (k=0; k<2; k++) {
1176 tol = Tolerance(k);
1177 dis = Disable(k);
1178 for (i=0; i<fNofMaps/2; i++) {
1179 Bool_t firstSignal=kTRUE;
1180 Int_t idx=i+k*fNofMaps/2;
1181 if( !fAnodeFire[idx] ) continue;
1182 CompressionParam(idx,decr,thres);
1183
1184 for (j=0; j<fMaxNofSamples; j++) {
1185 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1186 signal -= decr; // if baseline eq.
1187 if(do10to8) signal = Convert10to8(signal);
1188 if (signal <= thres) {
1189 signal=0;
1190 diff=128;
1191 last=0;
1192 // write diff in the buffer for HuffT
1193 str[counter]=(UChar_t)diff;
1194 counter++;
1195 continue;
1196 } // end if signal <= thres
1197 diff=signal-last;
1198 if (diff > 127) diff=127;
1199 if (diff < -128) diff=-128;
1200 if (signal < dis) {
1201 // tol has changed to 8 possible cases ? - one can write
1202 // this if(TMath::Abs(diff)<tol) ... else ...
1203 if(TMath::Abs(diff)<tol) diff=0;
1204 // or keep it as it was before
1205 AddDigit(idx,j,last+diff);
1206 } else {
1207 AddDigit(idx,j,signal);
1208 } // end if singal < dis
1209 diff += 128;
1210 // write diff in the buffer used to compute Huffman tables
1211 if (firstSignal) str[counter]=(UChar_t)signal;
1212 else str[counter]=(UChar_t)diff;
1213 counter++;
1214 last=signal;
1215 firstSignal=kFALSE;
1216 } // end for j loop time samples
1217 } // end for i loop anodes one half of detector
1218 } // end for k
b0f5e3fc 1219
1220 // check
aacedc3e 1221 fStream->CheckCount(counter);
b0f5e3fc 1222
aacedc3e 1223 // open file and write out the stream of diff's
1224 static Bool_t open=kTRUE;
1225 static TFile *outFile;
8ba39da9 1226 Bool_t write = res->OutputOption();
aacedc3e 1227 TDirectory *savedir = gDirectory;
b0f5e3fc 1228
aacedc3e 1229 if (write ) {
1230 if(open) {
1231 SetFileName("stream.root");
1232 cout<<"filename "<<fFileName<<endl;
1233 outFile=new TFile(fFileName,"recreate");
1234 cout<<"I have opened "<<fFileName<<" file "<<endl;
1235 } // end if open
1236 open = kFALSE;
1237 outFile->cd();
1238 fStream->Write();
1239 } // endif write
1240
1241 fStream->ClearStream();
1242
1243 // back to galice.root file
1244 if(savedir) savedir->cd();
8a33ae9e 1245}
1246//______________________________________________________________________
b0f5e3fc 1247void AliITSsimulationSDD::StoreAllDigits(){
aacedc3e 1248 // if non-zero-suppressed data
f45f6658 1249 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
8ba39da9 1250
1251 Bool_t do10to8 = res->Do10to8();
aacedc3e 1252 Int_t i, j, digits[3];
1253
1254 for (i=0; i<fNofMaps; i++) {
1255 for (j=0; j<fMaxNofSamples; j++) {
1256 Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1257 if(do10to8) signal = Convert10to8(signal);
1258 digits[0] = i;
1259 digits[1] = j;
1260 digits[2] = signal;
1261 fITS->AddRealDigit(1,digits);
0599a018 1262 } // end for j
aacedc3e 1263 } // end for i
b0f5e3fc 1264}
8a33ae9e 1265//______________________________________________________________________
ece86d9a 1266void AliITSsimulationSDD::CreateHistograms(Int_t scale){
aacedc3e 1267 // Creates histograms of maps for debugging
1268 Int_t i;
1269
1270 fHis=new TObjArray(fNofMaps);
1271 for (i=0;i<fNofMaps;i++) {
1272 TString sddName("sdd_");
1273 Char_t candNum[4];
1274 sprintf(candNum,"%d",i+1);
1275 sddName.Append(candNum);
1276 fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1277 0.,(Float_t) scale*fMaxNofSamples), i);
1278 } // end for i
b0f5e3fc 1279}
8a33ae9e 1280//______________________________________________________________________
ece86d9a 1281void AliITSsimulationSDD::FillHistograms(){
aacedc3e 1282 // fill 1D histograms from map
8a33ae9e 1283
aacedc3e 1284 if (!fHis) return;
8a33ae9e 1285
aacedc3e 1286 for( Int_t i=0; i<fNofMaps; i++) {
1287 TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1288 Int_t nsamples = hist->GetNbinsX();
1289 for( Int_t j=0; j<nsamples; j++) {
1290 Double_t signal=fHitMap2->GetSignal(i,j);
1291 hist->Fill((Float_t)j,signal);
1292 } // end for j
1293 } // end for i
ece86d9a 1294}
8a33ae9e 1295//______________________________________________________________________
b0f5e3fc 1296void AliITSsimulationSDD::ResetHistograms(){
aacedc3e 1297 // Reset histograms for this detector
1298 Int_t i;
8a33ae9e 1299
aacedc3e 1300 for (i=0;i<fNofMaps;i++ ) {
1301 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();
1302 } // end for i
b0f5e3fc 1303}
8a33ae9e 1304//______________________________________________________________________
b0f5e3fc 1305TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
aacedc3e 1306 // Fills a histogram from a give anode.
8a33ae9e 1307
aacedc3e 1308 if (!fHis) return 0;
8a33ae9e 1309
aacedc3e 1310 if(wing <=0 || wing > 2) {
1311 Warning("GetAnode","Wrong wing number: %d",wing);
1312 return NULL;
1313 } // end if wing <=0 || wing >2
1314 if(anode <=0 || anode > fNofMaps/2) {
1315 Warning("GetAnode","Wrong anode number: %d",anode);
1316 return NULL;
1317 } // end if ampde <=0 || andoe > fNofMaps/2
8a33ae9e 1318
aacedc3e 1319 Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1320 return (TH1F*)(fHis->At(index));
b0f5e3fc 1321}
8a33ae9e 1322//______________________________________________________________________
b0f5e3fc 1323void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
aacedc3e 1324 // Writes the histograms to a file
b0f5e3fc 1325
aacedc3e 1326 if (!fHis) return;
8a33ae9e 1327
aacedc3e 1328 hfile->cd();
1329 Int_t i;
1330 for(i=0; i<fNofMaps; i++) fHis->At(i)->Write(); //fAdcs[i]->Write();
1331 return;
b0f5e3fc 1332}
8a33ae9e 1333//______________________________________________________________________
ece86d9a 1334Float_t AliITSsimulationSDD::GetNoise() {
aacedc3e 1335 // Returns the noise value
1336 //Bool_t do10to8=GetResp()->Do10to8();
1337 //noise will always be in the liniar part of the signal
1338 Int_t decr;
1339 Int_t threshold = fT1[0];
1340 char opt1[20], opt2[20];
f45f6658 1341 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1342 SetCompressParam();
fcf95fc7 1343 res->GetParamOptions(opt1,opt2);
aacedc3e 1344 fParam=opt2;
aacedc3e 1345 Double_t noise,baseline;
f45f6658 1346 //GetBaseline(fModule);
aacedc3e 1347
1348 TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1349 if(c2) delete c2->GetPrimitive("noisehist");
1350 if(c2) delete c2->GetPrimitive("anode");
1351 else c2=new TCanvas("c2");
1352 c2->cd();
1353 c2->SetFillColor(0);
1354
1355 TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1356 TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1357 (float)fMaxNofSamples);
1358 Int_t i,k;
1359 for (i=0;i<fNofMaps;i++) {
1360 CompressionParam(i,decr,threshold);
f45f6658 1361 baseline = res->GetBaseline(i);
1362 noise = res->GetNoise(i);
aacedc3e 1363 anode->Reset();
1364 for (k=0;k<fMaxNofSamples;k++) {
1365 Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1366 //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1367 if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1368 anode->Fill((float)k,signal);
1369 } // end for k
1370 anode->Draw();
1371 c2->Update();
1372 } // end for i
1373 TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1374 noisehist->Fit("gnoise","RQ");
1375 noisehist->Draw();
ece86d9a 1376 c2->Update();
aacedc3e 1377 Float_t mnoise = gnoise->GetParameter(1);
1378 cout << "mnoise : " << mnoise << endl;
1379 Float_t rnoise = gnoise->GetParameter(2);
1380 cout << "rnoise : " << rnoise << endl;
1381 delete noisehist;
1382 return rnoise;
50d05d7b 1383}
1384//______________________________________________________________________
1385void AliITSsimulationSDD::WriteSDigits(){
aacedc3e 1386 // Fills the Summable digits Tree
1387 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1388
1389 for( Int_t i=0; i<fNofMaps; i++ ) {
1390 if( !fAnodeFire[i] ) continue;
f6b6d58e 1391 for( Int_t j=0; j<fMaxNofSamples; j++ ) {
aacedc3e 1392 Double_t sig = fHitMap2->GetSignal( i, j );
1393 if( sig > 0.2 ) {
1394 Int_t jdx = j*fScaleSize;
1395 Int_t index = fpList->GetHitIndex( i, j );
1396 AliITSpListItem pItemTmp2( fModule, index, 0. );
1397 // put the fScaleSize analog digits in only one
1398 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1399 AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1400 if( pItemTmp == 0 ) continue;
1401 pItemTmp2.Add( pItemTmp );
1402 }
1403 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1404 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1405 aliITS->AddSumDigit( pItemTmp2 );
1406 } // end if (sig > 0.2)
1407 }
48058160 1408 }
aacedc3e 1409 return;
b0f5e3fc 1410}
8a33ae9e 1411//______________________________________________________________________
d2f55a22 1412void AliITSsimulationSDD::PrintStatus() const {
aacedc3e 1413 // Print SDD simulation Parameters
1414
1415 cout << "**************************************************" << endl;
1416 cout << " Silicon Drift Detector Simulation Parameters " << endl;
1417 cout << "**************************************************" << endl;
1418 cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1419 cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1420 cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1421 cout << "Number pf Anodes used: " << fNofMaps << endl;
1422 cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1423 cout << "Scale size factor: " << fScaleSize << endl;
1424 cout << "**************************************************" << endl;
44a312c3 1425}