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