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