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