1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include <Riostream.h>
30 #include "AliITSMapA2.h"
31 #include "AliITSRawData.h"
32 #include "AliITSdigitSPD.h"
33 #include "AliITSetfSDD.h"
34 #include "AliITSmodule.h"
35 #include "AliITSpList.h"
36 #include "AliITSresponseSDD.h"
37 #include "AliITSCalibrationSDD.h"
38 #include "AliITSsegmentationSDD.h"
39 #include "AliITSsimulationSDD.h"
43 ClassImp(AliITSsimulationSDD)
44 ////////////////////////////////////////////////////////////////////////
46 // Written by Piergiorgio Cerello //
47 // November 23 1999 //
49 // AliITSsimulationSDD is the simulation of SDDs. //
50 ////////////////////////////////////////////////////////////////////////
52 //______________________________________________________________________
53 Int_t power(Int_t b, Int_t e) {
54 // compute b to the e power, where both b and e are Int_ts.
57 for(i=0; i<e; i++) power *= b;
60 //______________________________________________________________________
61 void FastFourierTransform(AliITSetfSDD *alisddetf,Double_t *real,
62 Double_t *imag,Int_t direction) {
63 // Do a Fast Fourier Transform
65 Int_t samples = alisddetf->GetSamples();
66 Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
69 Int_t m2 = samples/m1;
72 for(j=0; j<samples; j += m1) {
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);
92 for(j=0; j<samples; j++) {
96 for(i1=1; i1<=l; i1++) {
99 p = p + p + j2 - j1 - j1;
102 Double_t xr = *(real+j);
103 Double_t xi = *(imag+j);
104 *(real+j) = *(real+p);
105 *(imag+j) = *(imag+p);
110 if(direction == -1) {
111 for(i=0; i<samples; i++) {
112 *(real+i) /= samples;
113 *(imag+i) /= samples;
115 } // end if direction == -1
118 //______________________________________________________________________
119 AliITSsimulationSDD::AliITSsimulationSDD():
142 fCrosstalkFlag(kFALSE),
147 // Default constructor
149 SetPerpendTracksFlag();
154 //______________________________________________________________________
155 AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
156 AliITSsimulation(source),
158 fHitMap2(source.fHitMap2),
159 fHitSigMap2(source.fHitSigMap2),
160 fHitNoiMap2(source.fHitNoiMap2),
161 fStream(source.fStream),
162 fElectronics(source.fElectronics),
165 fOutZR(source.fOutZR),
166 fOutZI(source.fOutZI),
167 fAnodeFire(source.fAnodeFire),
173 fTreeB(source.fTreeB),
174 fParam(source.fParam),
175 fFileName(source.fFileName),
177 fCheckNoise(source.fCheckNoise),
178 fCrosstalkFlag(source.fCrosstalkFlag),
179 fDoFFT(source.fDoFFT),
180 fNofMaps(source.fNofMaps),
181 fMaxNofSamples(source.fMaxNofSamples),
182 fScaleSize(source.fScaleSize){
183 // Copy constructor to satify Coding roules only.
186 //______________________________________________________________________
187 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
188 // Assignment operator to satify Coding roules only.
190 if(this==&src) return *this;
191 Error("AliITSsimulationSDD","Not allowed to make a = with "
192 "AliITSsimulationSDD Using default creater instead");
195 //______________________________________________________________________
196 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
197 // Assignment operator to satify Coding roules only.
199 if(this==&src) return *this;
200 Error("AliITSsimulationSSD","Not allowed to make a = with "
201 "AliITSsimulationSDD Using default creater instead");
205 //______________________________________________________________________
206 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
207 AliITSsimulation(dettyp),
229 fCrosstalkFlag(kFALSE),
234 // Default Constructor
237 //______________________________________________________________________
238 void AliITSsimulationSDD::Init(){
239 // Standard Constructor
242 SetPerpendTracksFlag();
247 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
249 AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
250 fpList = new AliITSpList( seg->Npz(),
251 fScaleSize*seg->Npx() );
252 fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
253 fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
254 fHitMap2 = fHitSigMap2;
256 fNofMaps = seg->Npz();
257 fMaxNofSamples = seg->Npx();
258 fAnodeFire = new Bool_t [fNofMaps];
260 Float_t sddLength = seg->Dx();
261 Float_t sddWidth = seg->Dz();
264 Float_t anodePitch = seg->Dpz(dummy);
265 Double_t timeStep = (Double_t)seg->Dpx(dummy);
266 Float_t driftSpeed = res->DriftSpeed();
268 if(anodePitch*(fNofMaps/2) > sddWidth) {
269 Warning("AliITSsimulationSDD",
270 "Too many anodes %d or too big pitch %f \n",
271 fNofMaps/2,anodePitch);
274 if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
275 Error("AliITSsimulationSDD",
276 "Time Interval > Allowed Time Interval: exit\n");
280 fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
283 char opt1[20], opt2[20];
284 res->ParamOptions(opt1,opt2);
287 const char *kopt=res->ZeroSuppOption();
293 Bool_t write = res->OutputOption();
294 if(write && strstr(kopt,"2D")) MakeTreeB();
296 fITS = (AliITS*)gAlice->GetModule("ITS");
297 Int_t size = fNofMaps*fMaxNofSamples;
298 fStream = new AliITSInStream(size);
300 fInZR = new Double_t [fScaleSize*fMaxNofSamples];
301 fInZI = new Double_t [fScaleSize*fMaxNofSamples];
302 fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
303 fOutZI = new Double_t [fScaleSize*fMaxNofSamples];
305 //______________________________________________________________________
306 AliITSsimulationSDD::~AliITSsimulationSDD() {
321 if(fTreeB) delete fTreeB;
322 if(fInZR) delete [] fInZR;
323 if(fInZI) delete [] fInZI;
324 if(fOutZR) delete [] fOutZR;
325 if(fOutZI) delete [] fOutZI;
326 if(fAnodeFire) delete [] fAnodeFire;
328 //______________________________________________________________________
329 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
330 // create maps to build the lists of tracks for each summable digit
334 memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);
336 //______________________________________________________________________
337 void AliITSsimulationSDD::ClearMaps() {
340 fHitSigMap2->ClearMap();
341 fHitNoiMap2->ClearMap();
343 //______________________________________________________________________
344 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
345 // digitize module using the "slow" detector simulator creating
348 TObjArray *fHits = mod->GetHits();
349 Int_t nhits = fHits->GetEntriesFast();
352 InitSimulationModule( md, ev );
353 HitsToAnalogDigits( mod );
354 ChargeToSignal( fModule,kFALSE ); // - Process signal without add noise
355 fHitMap2 = fHitNoiMap2; // - Swap to noise map
356 ChargeToSignal( fModule,kTRUE ); // - Process only noise
357 fHitMap2 = fHitSigMap2; // - Return to signal map
361 //______________________________________________________________________
362 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
364 // Add Summable digits to module maps.
365 AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
366 Int_t nItems = pItemArray->GetEntries();
367 Double_t maxadc = res->MaxAdc();
370 // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
371 for( Int_t i=0; i<nItems; i++ ) {
372 AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
373 if( pItem->GetModule() != fModule ) {
374 Error( "AliITSsimulationSDD","Error reading, SDigits module "
375 "%d != current module %d: exit",
376 pItem->GetModule(), fModule );
380 if(pItem->GetSignal()>0.0 ) sig = kTRUE;
382 fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
383 AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
384 Double_t sigAE = pItem2->GetSignalAfterElect();
385 if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
388 fpList->GetMapIndex( pItem->GetIndex(), ia, it );
389 fHitMap2->SetHit( ia, it, sigAE );
390 fAnodeFire[ia] = kTRUE;
394 //______________________________________________________________________
395 void AliITSsimulationSDD::FinishSDigitiseModule() {
396 // digitize module using the "slow" detector simulator from
397 // the sum of summable digits.
401 //______________________________________________________________________
402 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
403 // create maps to build the lists of tracks for each digit
405 TObjArray *fHits = mod->GetHits();
406 Int_t nhits = fHits->GetEntriesFast();
408 InitSimulationModule( md, ev );
410 if( !nhits && fCheckNoise ) {
411 ChargeToSignal( fModule,kTRUE ); // process noise
418 HitsToAnalogDigits( mod );
419 ChargeToSignal( fModule,kTRUE ); // process signal + noise
421 for( Int_t i=0; i<fNofMaps; i++ ) {
422 for( Int_t j=0; j<fMaxNofSamples; j++ ) {
423 Int_t jdx = j*fScaleSize;
424 Int_t index = fpList->GetHitIndex( i, j );
425 AliITSpListItem pItemTmp2( fModule, index, 0. );
426 // put the fScaleSize analog digits in only one
427 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
428 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
429 if( pItemTmp == 0 ) continue;
430 pItemTmp2.Add( pItemTmp );
432 fpList->DeleteHit( i, j );
433 fpList->AddItemTo( 0, &pItemTmp2 );
439 //______________________________________________________________________
440 void AliITSsimulationSDD::FinishDigits() {
441 // introduce the electronics effects and do zero-suppression if required
443 ApplyDeadChannels(fModule);
444 if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
446 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
447 const char *kopt = res->GetZeroSuppOption();
448 ZeroSuppression( kopt );
450 //______________________________________________________________________
451 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
452 // create maps to build the lists of tracks for each digit
454 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
455 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
456 TObjArray *hits = mod->GetHits();
457 Int_t nhits = hits->GetEntriesFast();
459 // Int_t arg[6] = {0,0,0,0,0,0};
461 Int_t nofAnodes = fNofMaps/2;
462 Double_t sddLength = seg->Dx();
463 Double_t sddWidth = seg->Dz();
464 Double_t anodePitch = seg->Dpz(dummy);
465 Double_t timeStep = seg->Dpx(dummy);
466 Double_t driftSpeed = res->GetDriftSpeed();
467 //Float_t maxadc = res->GetMaxAdc();
468 //Float_t topValue = res->GetDynamicRange();
469 Double_t norm = res->GetMaxAdc()/res->GetDynamicRange(); // maxadc/topValue;
470 Double_t cHloss = res->GetChargeLoss();
471 Float_t dfCoeff, s1; res->DiffCoeff(dfCoeff,s1); // Signal 2d Shape
472 Double_t eVpairs = res->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
473 Double_t nsigma = res->GetNSigmaIntegration(); //
474 Int_t nlookups = res->GetGausNLookUp(); //
475 Float_t jitter = res->GetJitterError(); //
477 // Piergiorgio's part (apart for few variables which I made float
478 // when i thought that can be done
479 // Fill detector maps with GEANT hits
480 // loop over hits in the module
482 const Float_t kconv = 1.0e+6; // GeV->KeV
484 Int_t hitDetector; // detector number (lay,lad,hitDetector)
485 Int_t iWing; // which detector wing/side.
486 Int_t detector; // 2*(detector-1)+iWing
487 Int_t ii,kk,ka,kt; // loop indexs
488 Int_t ia,it,index; // sub-pixel integration indexies
489 Int_t iAnode; // anode number.
490 Int_t timeSample; // time buckett.
491 Int_t anodeWindow; // anode direction charge integration width
492 Int_t timeWindow; // time direction charge integration width
493 Int_t jamin,jamax; // anode charge integration window
494 Int_t jtmin,jtmax; // time charge integration window
495 Int_t ndiv; // Anode window division factor.
496 Int_t nsplit; // the number of splits in anode and time windows==1.
497 Int_t nOfSplits; // number of times track length is split into
498 Float_t nOfSplitsF; // Floating point version of nOfSplits.
499 Float_t kkF; // Floating point version of loop index kk.
500 Double_t pathInSDD; // Track length in SDD.
501 Double_t drPath; // average position of track in detector. in microns
502 Double_t drTime; // Drift time
503 Double_t nmul; // drift time window multiplication factor.
504 Double_t avDrft; // x position of path length segment in cm.
505 Double_t avAnode; // Anode for path length segment in Anode number (float)
506 Double_t xAnode; // Floating point anode number.
507 Double_t driftPath; // avDrft in microns.
508 Double_t width; // width of signal at anodes.
509 Double_t depEnergy; // Energy deposited in this GEANT step.
510 Double_t xL[3],dxL[3]; // local hit coordinates and diff.
511 Double_t sigA; // sigma of signal at anode.
512 Double_t sigT; // sigma in time/drift direction for track segment
513 Double_t aStep,aConst; // sub-pixel size and offset anode
514 Double_t tStep,tConst; // sub-pixel size and offset time
515 Double_t amplitude; // signal amplitude for track segment in nanoAmpere
516 Double_t chargeloss; // charge loss for track segment.
517 Double_t anodeAmplitude; // signal amplitude in anode direction
518 Double_t aExpo; // exponent of Gaussian anode direction
519 Double_t timeAmplitude; // signal amplitude in time direction
520 Double_t tExpo; // exponent of Gaussian time direction
521 // Double_t tof; // Time of flight in ns of this step.
523 for(ii=0; ii<nhits; ii++) {
524 if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
525 depEnergy,itrack)) continue;
526 xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
528 hitDetector = mod->GetDet();
529 //tof = 1.E+09*(mod->GetHit(ii)->GetTOF()); // tof in ns.
530 //if(tof>sddLength/driftSpeed) continue; // hit happed too late.
532 // scale path to simulate a perpendicular track
533 // continue if the particle did not lose energy
534 // passing through detector
537 Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
538 itrack,ii,mod->GetIndex()));
540 } // end if !depEnergy
542 pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
544 if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
545 drPath = 10000.*(dxL[0]+2.*xL[0])*0.5;
546 if(drPath < 0) drPath = -drPath;
547 drPath = sddLength-drPath;
549 AliDebug(1, // this should be fixed at geometry level
550 Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
551 drPath,sddLength,dxL[0],xL[0]));
553 } // end if drPath < 0
555 // Compute number of segments to brake step path into
556 drTime = drPath/driftSpeed; // Drift Time
557 sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
558 // calcuate the number of time the path length should be split into.
559 nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
560 if(fFlag) nOfSplits = 1;
562 // loop over path segments, init. some variables.
563 depEnergy /= nOfSplits;
564 nOfSplitsF = (Float_t) nOfSplits;
565 for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
566 kkF = (Float_t) kk + 0.5;
567 avDrft = xL[0]+dxL[0]*kkF/nOfSplitsF;
568 avAnode = xL[2]+dxL[2]*kkF/nOfSplitsF;
569 driftPath = 10000.*avDrft;
571 iWing = 2; // Assume wing is 2
572 if(driftPath < 0) { // if wing is not 2 it is 1.
574 driftPath = -driftPath;
575 } // end if driftPath < 0
576 driftPath = sddLength-driftPath;
577 detector = 2*(hitDetector-1) + iWing;
579 AliDebug(1, // this should be fixed at geometry level
580 Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
581 driftPath,sddLength,avDrft,dxL[0],xL[0]));
583 } // end if driftPath < 0
586 drTime = driftPath/driftSpeed; // drift time for segment.
587 timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1);
588 // compute time Sample including tof information. The tof only
589 // effects the time of the signal is recoreded and not the
591 // timeSample = (Int_t) (fScaleSize*(drTime+tof)/timeStep + 1);
592 if(timeSample > fScaleSize*fMaxNofSamples) {
593 Warning("HitsToAnalogDigits","Wrong Time Sample: %e",
596 } // end if timeSample > fScaleSize*fMaxNoofSamples
599 xAnode = 10000.*(avAnode)/anodePitch + nofAnodes/2; // +1?
600 if(xAnode*anodePitch > sddWidth || xAnode*anodePitch < 0.)
601 Warning("HitsToAnalogDigits",
602 "Exceedubg sddWidth=%e Z = %e",
603 sddWidth,xAnode*anodePitch);
604 iAnode = (Int_t) (1.+xAnode); // xAnode?
605 if(iAnode < 1 || iAnode > nofAnodes) {
606 Warning("HitToAnalogDigits","Wrong iAnode: 1<%d>%d (xanode=%e)",
607 iAnode,nofAnodes, xAnode);
609 } // end if iAnode < 1 || iAnode > nofAnodes
611 // store straight away the particle position in the array
612 // of particles and take idhit=ii only when part is entering (this
613 // requires FillModules() in the macro for analysis) :
615 // Sigma along the anodes for track segment.
616 sigA = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
617 sigT = sigA/driftSpeed;
618 // Peak amplitude in nanoAmpere
619 amplitude = fScaleSize*160.*depEnergy/
620 (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA);
621 amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to
622 // account for clock variations
623 // (reference value: 40 MHz)
624 chargeloss = 1.-cHloss*driftPath/1000;
625 amplitude *= chargeloss;
626 width = 2.*nsigma/(nlookups-1);
634 } // end if drTime > 1200.
636 nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2;
637 // Sub-pixel size see computation of aExpo and tExpo.
638 aStep = anodePitch/(nsplit*fScaleSize*sigA);
639 aConst = xAnode*anodePitch/sigA;
640 tStep = timeStep/(nsplit*fScaleSize*sigT);
641 tConst = drTime/sigT;
642 // Define SDD window corresponding to the hit
643 anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1);
644 timeWindow = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
645 jamin = (iAnode - anodeWindow/ndiv - 1)*fScaleSize*nsplit +1;
646 jamax = (iAnode + anodeWindow/ndiv)*fScaleSize*nsplit;
647 if(jamin <= 0) jamin = 1;
648 if(jamax > fScaleSize*nofAnodes*nsplit)
649 jamax = fScaleSize*nofAnodes*nsplit;
650 // jtmin and jtmax are Hard-wired
651 jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1;
652 jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit;
653 if(jtmin <= 0) jtmin = 1;
654 if(jtmax > fScaleSize*fMaxNofSamples*nsplit)
655 jtmax = fScaleSize*fMaxNofSamples*nsplit;
656 // Spread the charge in the anode-time window
657 for(ka=jamin; ka <=jamax; ka++) {
658 ia = (ka-1)/(fScaleSize*nsplit) + 1;
660 Warning("HitsToAnalogDigits","ia < 1: ");
663 if(ia > nofAnodes) ia = nofAnodes;
664 aExpo = (aStep*(ka-0.5)-aConst);
665 if(TMath::Abs(aExpo) > nsigma) anodeAmplitude = 0.;
667 dummy = (Int_t) ((aExpo+nsigma)/width);
668 anodeAmplitude = amplitude*res->GetGausLookUp(dummy);
669 } // end if TMath::Abs(aEspo) > nsigma
670 // index starts from 0
671 index = ((detector+1)%2)*nofAnodes+ia-1;
672 if(anodeAmplitude) for(kt=jtmin; kt<=jtmax; kt++) {
673 it = (kt-1)/nsplit+1; // it starts from 1
675 Warning("HitsToAnalogDigits","it < 1:");
678 if(it>fScaleSize*fMaxNofSamples)
679 it = fScaleSize*fMaxNofSamples;
680 tExpo = (tStep*(kt-0.5)-tConst);
681 if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
683 dummy = (Int_t) ((tExpo+nsigma)/width);
684 timeAmplitude = anodeAmplitude*
685 res->GetGausLookUp(dummy);
686 } // end if TMath::Abs(tExpo) > nsigma
687 // build the list of Sdigits for this module
690 // arg[2] = itrack; // track number
691 // arg[3] = ii-1; // hit number.
692 timeAmplitude *= norm;
694 // ListOfFiredCells(arg,timeAmplitude,alst,padr);
695 Double_t charge = timeAmplitude;
696 charge += fHitMap2->GetSignal(index,it-1);
697 fHitMap2->SetHit(index, it-1, charge);
698 fpList->AddSignal(index,it-1,itrack,ii-1,
699 mod->GetIndex(),timeAmplitude);
700 fAnodeFire[index] = kTRUE;
701 } // end if anodeAmplitude and loop over time in window
702 } // loop over anodes in window
703 } // end loop over "sub-hits"
704 } // end loop over hits
707 //______________________________________________________________________
708 void AliITSsimulationSDD::ListOfFiredCells(Int_t *arg,Double_t timeAmplitude,
709 TObjArray *alist,TClonesArray *padr){
710 // Returns the list of "fired" cells.
712 Int_t index = arg[0];
714 Int_t idtrack = arg[2];
715 Int_t idhit = arg[3];
716 Int_t counter = arg[4];
717 Int_t countadr = arg[5];
718 Double_t charge = timeAmplitude;
719 charge += fHitMap2->GetSignal(index,ik-1);
720 fHitMap2->SetHit(index, ik-1, charge);
723 Int_t it = (Int_t)((ik-1)/fScaleSize);
726 digits[2] = (Int_t)timeAmplitude;
728 if (idtrack >= 0) phys = (Float_t)timeAmplitude;
731 Double_t cellcharge = 0.;
732 AliITSTransientDigit* pdigit;
733 // build the list of fired cells and update the info
734 if (!fHitMap1->TestHit(index, it)) {
735 new((*padr)[countadr++]) TVector(3);
736 TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
737 trinfo(0) = (Float_t)idtrack;
738 trinfo(1) = (Float_t)idhit;
739 trinfo(2) = (Float_t)timeAmplitude;
741 alist->AddAtAndExpand(new AliITSTransientDigit(phys,digits),counter);
742 fHitMap1->SetHit(index, it, counter);
744 pdigit=(AliITSTransientDigit*)alist->At(alist->GetLast());
746 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
747 trlist->Add(&trinfo);
749 pdigit = (AliITSTransientDigit*) fHitMap1->GetHit(index, it);
750 for(Int_t kk=0;kk<fScaleSize;kk++) {
751 cellcharge += fHitMap2->GetSignal(index,fScaleSize*it+kk);
754 (*pdigit).fSignal = (Int_t)cellcharge;
755 (*pdigit).fPhysics += phys;
756 // update list of tracks
757 TObjArray* trlist = (TObjArray*)pdigit->TrackList();
758 Int_t lastentry = trlist->GetLast();
759 TVector *ptrkp = (TVector*)trlist->At(lastentry);
760 TVector &trinfo = *ptrkp;
761 Int_t lasttrack = Int_t(trinfo(0));
762 Float_t lastcharge=(trinfo(2));
763 if (lasttrack==idtrack ) {
764 lastcharge += (Float_t)timeAmplitude;
765 trlist->RemoveAt(lastentry);
766 trinfo(0) = lasttrack;
768 trinfo(2) = lastcharge;
769 trlist->AddAt(&trinfo,lastentry);
771 new((*padr)[countadr++]) TVector(3);
772 TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
773 trinfo(0) = (Float_t)idtrack;
774 trinfo(1) = (Float_t)idhit;
775 trinfo(2) = (Float_t)timeAmplitude;
776 trlist->Add(&trinfo);
777 } // end if lasttrack==idtrack
780 // check the track list - debugging
781 Int_t trk[20], htrk[20];
783 Int_t nptracks = trlist->GetEntriesFast();
786 for (tr=0;tr<nptracks;tr++) {
787 TVector *pptrkp = (TVector*)trlist->At(tr);
788 TVector &pptrk = *pptrkp;
789 trk[tr] = Int_t(pptrk(0));
790 htrk[tr] = Int_t(pptrk(1));
791 chtrk[tr] = (pptrk(2));
792 cout << "nptracks "<<nptracks << endl;
795 } // end if AliDebugLevel()
798 // update counter and countadr for next call.
803 //____________________________________________
804 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signal ) {
806 Int_t size = AliITSdigit::GetNTracks();
809 Int_t * tracks = new Int_t[size];
810 Int_t * hits = new Int_t[size];
812 Float_t * charges = new Float_t[size];
818 AliITSpListItem *pItem = fpList->GetpListItem( i, j );
821 for( Int_t l=0; l<size; l++ ) {
827 Int_t idtrack = pItem->GetTrack( 0 );
828 if( idtrack >= 0 ) phys = pItem->GetSignal();
831 for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
832 tracks[l] = pItem->GetTrack( l );
833 hits[l] = pItem->GetHit( l );
834 charges[l] = pItem->GetSignal( l );
842 fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges );
847 //______________________________________________________________________
848 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise) {
849 // add baseline, noise, electronics and ADC saturation effects
851 char opt1[20], opt2[20];
852 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
853 res->GetParamOptions(opt1,opt2);
859 Float_t maxadc = res->GetMaxAdc();
861 for (i=0;i<fNofMaps;i++) {
862 if( !fAnodeFire[i] ) continue;
863 baseline = res->GetBaseline(i);
864 noise = res->GetNoise(i);
866 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
867 fInZR[k] = fHitMap2->GetSignal(i,k);
869 contrib = (baseline + noise*gRandom->Gaus());
873 for(k=0; k<fMaxNofSamples; k++) {
874 Double_t newcont = 0.;
875 Double_t maxcont = 0.;
876 for(kk=0;kk<fScaleSize;kk++) {
877 newcont = fInZR[fScaleSize*k+kk];
878 if(newcont > maxcont) maxcont = newcont;
881 if (newcont >= maxadc) newcont = maxadc -1;
882 if(newcont >= baseline){
883 Warning("","newcont=%d>=baseline=%d",newcont,baseline);
886 fHitMap2->SetHit(i,k,newcont);
888 } // end for i loop over anodes
892 for (i=0;i<fNofMaps;i++) {
893 if( !fAnodeFire[i] ) continue;
894 baseline = res->GetBaseline(i);
895 noise = res->GetNoise(i);
896 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
897 fInZR[k] = fHitMap2->GetSignal(i,k);
899 contrib = (baseline + noise*gRandom->Gaus());
904 FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
905 for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
906 Double_t rw = fElectronics->GetTraFunReal(k);
907 Double_t iw = fElectronics->GetTraFunImag(k);
908 fOutZR[k] = fInZR[k]*rw - fInZI[k]*iw;
909 fOutZI[k] = fInZR[k]*iw + fInZI[k]*rw;
911 FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
912 for(k=0; k<fMaxNofSamples; k++) {
913 Double_t newcont1 = 0.;
914 Double_t maxcont1 = 0.;
915 for(kk=0;kk<fScaleSize;kk++) {
916 newcont1 = fOutZR[fScaleSize*k+kk];
917 if(newcont1 > maxcont1) maxcont1 = newcont1;
920 if (newcont1 >= maxadc) newcont1 = maxadc -1;
921 fHitMap2->SetHit(i,k,newcont1);
923 } // end for i loop over anodes
926 //____________________________________________________________________
927 void AliITSsimulationSDD::ApplyDeadChannels(Int_t mod) {
928 // Set dead channel signal to zero
929 AliITSCalibrationSDD * calibr = (AliITSCalibrationSDD *)GetCalibrationModel(mod);
930 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
932 if( calibr->IsDead() ||
933 ( calibr->GetDeadChips() == 0 &&
934 calibr->GetDeadChannels() == 0 ) ) return;
936 // static AliITS *iTS = (AliITS*)gAlice->GetModule( "ITS" );
938 Int_t fMaxNofSamples = seg->Npx();
939 // AliITSgeom *geom = iTS->GetITSgeom();
940 // Int_t firstSDDMod = geom->GetStartDet( 1 );
942 for( Int_t j=0; j<2; j++ ) {
943 // Int_t mod = (fModule-firstSDDMod)*2 + j;
944 for( Int_t u=0; u<calibr->Chips(); u++ )
945 for( Int_t v=0; v<calibr->Channels(); v++ ) {
946 Float_t gain = calibr->Gain(j, u, v );
947 for( Int_t k=0; k<fMaxNofSamples; k++ ) {
948 Int_t i = j*calibr->Chips()*calibr->Channels() +
949 u*calibr->Channels() +
951 Double_t signal = gain * fHitMap2->GetSignal( i, k );
952 fHitMap2->SetHit( i, k, signal ); ///
957 //______________________________________________________________________
958 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
959 // function add the crosstalk effect to signal
960 // temporal function, should be checked...!!!
961 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
963 Int_t fNofMaps = seg->Npz();
964 Int_t fMaxNofSamples = seg->Npx();
966 // create and inizialice crosstalk map
967 Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
969 Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
972 memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
973 AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
974 for( Int_t z=0; z<fNofMaps; z++ ) {
975 Double_t baseline = calibr->GetBaseline(z);
981 for( Int_t l=0; l<fMaxNofSamples; l++ ) {
982 Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
983 if( fadc > baseline ) {
984 if( on == kFALSE && l<fMaxNofSamples-4 ) {
985 Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
986 if( fadc1 < fadc ) continue;
993 else { // end fadc > baseline
997 // make smooth derivative
998 Float_t* dev = new Float_t[fMaxNofSamples+1];
999 memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
1001 Error( "ApplyCrosstalk",
1002 "no memory for temporal array: exit \n" );
1005 for( Int_t i=tstart; i<tstop; i++ ) {
1006 if( i > 2 && i < fMaxNofSamples-2 )
1007 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
1008 -0.1*fHitMap2->GetSignal( z,i-1 )
1009 +0.1*fHitMap2->GetSignal( z,i+1 )
1010 +0.2*fHitMap2->GetSignal( z,i+2 );
1013 // add crosstalk contribution to neibourg anodes
1014 for( Int_t i=tstart; i<tstop; i++ ) {
1015 Int_t anode = z - 1;
1016 Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
1017 Float_t ctktmp = -dev[i1] * 0.25;
1019 ctk[anode*fMaxNofSamples+i] += ctktmp;
1022 if( anode < fNofMaps ) {
1023 ctk[anode*fMaxNofSamples+i] += ctktmp;
1028 } // if( nTsteps > 2 )
1030 } // if( on == kTRUE )
1035 for( Int_t a=0; a<fNofMaps; a++ )
1036 for( Int_t t=0; t<fMaxNofSamples; t++ ) {
1037 Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
1038 fHitMap2->SetHit( a, t, signal );
1044 //______________________________________________________________________
1045 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
1047 // Returns the compression alogirthm parameters
1052 //______________________________________________________________________
1053 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{
1054 // returns the compression alogirthm parameters
1060 //______________________________________________________________________
1061 void AliITSsimulationSDD::SetCompressParam(){
1062 // Sets the compression alogirthm parameters
1063 AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1064 for(Int_t ian = 0; ian<fNofMaps;ian++){
1065 fD[ian] = (Int_t)(calibr->GetBaseline(ian));
1066 fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5);
1067 fT2[ian] = 0; // used by 2D clustering - not defined yet
1068 fTol[ian] = 0; // used by 2D clustering - not defined yet
1071 //______________________________________________________________________
1072 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
1073 // To the 10 to 8 bit lossive compression.
1074 // code from Davide C. and Albert W.
1076 if (signal < 128) return signal;
1077 if (signal < 256) return (128+((signal-128)>>1));
1078 if (signal < 512) return (192+((signal-256)>>3));
1079 if (signal < 1024) return (224+((signal-512)>>4));
1083 //______________________________________________________________________
1084 AliITSMap* AliITSsimulationSDD::HitMap(Int_t i){
1085 //Return the correct map.
1087 return ((i==0)? fHitMap1 : fHitMap2);
1090 //______________________________________________________________________
1091 void AliITSsimulationSDD::ZeroSuppression(const char *option) {
1092 // perform the zero suppresion
1093 if (strstr(option,"2D")) {
1094 //Init2D(); // activate if param change module by module
1096 } else if (strstr(option,"1D")) {
1097 //Init1D(); // activate if param change module by module
1099 } else StoreAllDigits();
1101 //______________________________________________________________________
1102 void AliITSsimulationSDD::Init2D(){
1103 // read in and prepare arrays: fD, fT1, fT2
1104 // savemu[nanodes], savesigma[nanodes]
1105 // read baseline and noise from file - either a .root file and in this
1106 // case data should be organised in a tree with one entry for each
1107 // module => reading should be done accordingly
1108 // or a classic file and do smth. like this ( code from Davide C. and
1110 // Read 2D zero-suppression parameters for SDD
1112 if (!strstr(fParam.Data(),"file")) return;
1114 Int_t na,pos,tempTh;
1116 Float_t *savemu = new Float_t [fNofMaps];
1117 Float_t *savesigma = new Float_t [fNofMaps];
1118 char input[100],basel[100],par[100];
1121 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1123 res->Thresholds(tmp1,tmp2);
1124 Int_t minval = static_cast<Int_t>(tmp1);
1126 res->Filenames(input,basel,par);
1129 filtmp = gSystem->ExpandPathName(fFileName.Data());
1130 FILE *param = fopen(filtmp,"r");
1134 while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1136 Error("Init2D","Anode number not in increasing order!",filtmp);
1138 } // end if pos != na+1
1140 savesigma[na] = sigma;
1141 if ((2.*sigma) < mu) {
1142 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1145 tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1146 if (tempTh < 0) tempTh=0;
1148 tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
1149 if (tempTh < 0) tempTh=0;
1154 Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1161 delete [] savesigma;
1163 //______________________________________________________________________
1164 void AliITSsimulationSDD::Compress2D(){
1165 // simple ITS cluster finder -- online zero-suppression conditions
1169 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1171 res->Thresholds(tmp1,tmp2);
1172 Int_t minval = static_cast<Int_t>(tmp1);
1173 Bool_t write = res->OutputOption();
1174 Bool_t do10to8 = res->Do10to8();
1175 Int_t nz, nl, nh, low, i, j;
1177 for (i=0; i<fNofMaps; i++) {
1178 CompressionParam(i,db,tl,th);
1183 for (j=0; j<fMaxNofSamples; j++) {
1184 Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1185 signal -= db; // if baseline eq. is done here
1186 if (signal <= 0) {nz++; continue;}
1187 if ((signal - tl) < minval) low++;
1188 if ((signal - th) >= minval) {
1191 FindCluster(i,j,signal,minval,cond);
1193 ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1194 if(do10to8) signal = Convert10to8(signal);
1195 AddDigit(i,j,signal);
1196 } // end if cond&&j&&()
1197 } else if ((signal - tl) >= minval) nl++;
1198 } // end for j loop time samples
1199 if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1200 } //end for i loop anodes
1204 sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1205 TreeB()->Write(hname);
1210 //______________________________________________________________________
1211 void AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
1212 Int_t minval,Bool_t &cond){
1213 // Find clusters according to the online 2D zero-suppression algorithm
1214 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1215 AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1217 Bool_t do10to8 = res->Do10to8();
1218 Bool_t high = kFALSE;
1220 fHitMap2->FlagHit(i,j);
1222 // check the online zero-suppression conditions
1224 const Int_t kMaxNeighbours = 4;
1227 Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1228 seg->Neighbours(i,j,&nn,xList,yList);
1230 for (in=0; in<nn; in++) {
1233 if (fHitMap2->TestHit(ix,iy)==kUnused) {
1234 CompressionParam(ix,dbx,tlx,thx);
1235 Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1236 qn -= dbx; // if baseline eq. is done here
1237 if ((qn-tlx) < minval) {
1238 fHitMap2->FlagHit(ix,iy);
1241 if ((qn - thx) >= minval) high=kTRUE;
1243 if(do10to8) signal = Convert10to8(signal);
1244 AddDigit(i,j,signal);
1246 if(do10to8) qns = Convert10to8(qn);
1248 if (!high) AddDigit(ix,iy,qns);
1250 if(!high) fHitMap2->FlagHit(ix,iy);
1251 } // end if qn-tlx < minval
1253 } // end for in loop over neighbours
1255 //______________________________________________________________________
1256 void AliITSsimulationSDD::Init1D(){
1257 // this is just a copy-paste of input taken from 2D algo
1258 // Torino people should give input
1259 // Read 1D zero-suppression parameters for SDD
1261 if (!strstr(fParam.Data(),"file")) return;
1263 Int_t na,pos,tempTh;
1265 Float_t *savemu = new Float_t [fNofMaps];
1266 Float_t *savesigma = new Float_t [fNofMaps];
1267 char input[100],basel[100],par[100];
1270 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1272 res->Thresholds(tmp1,tmp2);
1273 Int_t minval = static_cast<Int_t>(tmp1);
1275 res->Filenames(input,basel,par);
1278 // set first the disable and tol param
1281 filtmp = gSystem->ExpandPathName(fFileName.Data());
1282 FILE *param = fopen(filtmp,"r");
1286 fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1287 while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1289 Error("Init1D","Anode number not in increasing order!",filtmp);
1291 } // end if pos != na+1
1293 savesigma[na]=sigma;
1294 if ((2.*sigma) < mu) {
1295 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1298 tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1299 if (tempTh < 0) tempTh=0;
1304 Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1311 delete [] savesigma;
1313 //______________________________________________________________________
1314 void AliITSsimulationSDD::Compress1D(){
1315 // 1D zero-suppression algorithm (from Gianluca A.)
1316 Int_t dis,tol,thres,decr,diff;
1317 UChar_t *str=fStream->Stream();
1319 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1321 Bool_t do10to8=res->Do10to8();
1325 for (k=0; k<2; k++) {
1328 for (i=0; i<fNofMaps/2; i++) {
1329 Bool_t firstSignal=kTRUE;
1330 Int_t idx=i+k*fNofMaps/2;
1331 if( !fAnodeFire[idx] ) continue;
1332 CompressionParam(idx,decr,thres);
1334 for (j=0; j<fMaxNofSamples; j++) {
1335 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1336 signal -= decr; // if baseline eq.
1337 if(do10to8) signal = Convert10to8(signal);
1338 if (signal <= thres) {
1342 // write diff in the buffer for HuffT
1343 str[counter]=(UChar_t)diff;
1346 } // end if signal <= thres
1348 if (diff > 127) diff=127;
1349 if (diff < -128) diff=-128;
1351 // tol has changed to 8 possible cases ? - one can write
1352 // this if(TMath::Abs(diff)<tol) ... else ...
1353 if(TMath::Abs(diff)<tol) diff=0;
1354 // or keep it as it was before
1355 AddDigit(idx,j,last+diff);
1357 AddDigit(idx,j,signal);
1358 } // end if singal < dis
1360 // write diff in the buffer used to compute Huffman tables
1361 if (firstSignal) str[counter]=(UChar_t)signal;
1362 else str[counter]=(UChar_t)diff;
1366 } // end for j loop time samples
1367 } // end for i loop anodes one half of detector
1371 fStream->CheckCount(counter);
1373 // open file and write out the stream of diff's
1374 static Bool_t open=kTRUE;
1375 static TFile *outFile;
1376 Bool_t write = res->OutputOption();
1377 TDirectory *savedir = gDirectory;
1381 SetFileName("stream.root");
1382 cout<<"filename "<<fFileName<<endl;
1383 outFile=new TFile(fFileName,"recreate");
1384 cout<<"I have opened "<<fFileName<<" file "<<endl;
1391 fStream->ClearStream();
1393 // back to galice.root file
1394 if(savedir) savedir->cd();
1396 //______________________________________________________________________
1397 void AliITSsimulationSDD::StoreAllDigits(){
1398 // if non-zero-suppressed data
1399 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1401 Bool_t do10to8 = res->Do10to8();
1402 Int_t i, j, digits[3];
1404 for (i=0; i<fNofMaps; i++) {
1405 for (j=0; j<fMaxNofSamples; j++) {
1406 Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1407 if(do10to8) signal = Convert10to8(signal);
1411 fITS->AddRealDigit(1,digits);
1415 //______________________________________________________________________
1416 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
1417 // Creates histograms of maps for debugging
1420 fHis=new TObjArray(fNofMaps);
1421 for (i=0;i<fNofMaps;i++) {
1422 TString sddName("sdd_");
1424 sprintf(candNum,"%d",i+1);
1425 sddName.Append(candNum);
1426 fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1427 0.,(Float_t) scale*fMaxNofSamples), i);
1430 //______________________________________________________________________
1431 void AliITSsimulationSDD::FillHistograms(){
1432 // fill 1D histograms from map
1436 for( Int_t i=0; i<fNofMaps; i++) {
1437 TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1438 Int_t nsamples = hist->GetNbinsX();
1439 for( Int_t j=0; j<nsamples; j++) {
1440 Double_t signal=fHitMap2->GetSignal(i,j);
1441 hist->Fill((Float_t)j,signal);
1445 //______________________________________________________________________
1446 void AliITSsimulationSDD::ResetHistograms(){
1447 // Reset histograms for this detector
1450 for (i=0;i<fNofMaps;i++ ) {
1451 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();
1454 //______________________________________________________________________
1455 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
1456 // Fills a histogram from a give anode.
1458 if (!fHis) return 0;
1460 if(wing <=0 || wing > 2) {
1461 Warning("GetAnode","Wrong wing number: %d",wing);
1463 } // end if wing <=0 || wing >2
1464 if(anode <=0 || anode > fNofMaps/2) {
1465 Warning("GetAnode","Wrong anode number: %d",anode);
1467 } // end if ampde <=0 || andoe > fNofMaps/2
1469 Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1470 return (TH1F*)(fHis->At(index));
1472 //______________________________________________________________________
1473 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
1474 // Writes the histograms to a file
1480 for(i=0; i<fNofMaps; i++) fHis->At(i)->Write(); //fAdcs[i]->Write();
1483 //______________________________________________________________________
1484 Float_t AliITSsimulationSDD::GetNoise() {
1485 // Returns the noise value
1486 //Bool_t do10to8=GetResp()->Do10to8();
1487 //noise will always be in the liniar part of the signal
1489 Int_t threshold = fT1[0];
1490 char opt1[20], opt2[20];
1491 AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1493 res->GetParamOptions(opt1,opt2);
1495 Double_t noise,baseline;
1496 //GetBaseline(fModule);
1498 TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1499 if(c2) delete c2->GetPrimitive("noisehist");
1500 if(c2) delete c2->GetPrimitive("anode");
1501 else c2=new TCanvas("c2");
1503 c2->SetFillColor(0);
1505 TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1506 TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1507 (float)fMaxNofSamples);
1509 for (i=0;i<fNofMaps;i++) {
1510 CompressionParam(i,decr,threshold);
1511 baseline = res->GetBaseline(i);
1512 noise = res->GetNoise(i);
1514 for (k=0;k<fMaxNofSamples;k++) {
1515 Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1516 //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1517 if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1518 anode->Fill((float)k,signal);
1523 TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1524 noisehist->Fit("gnoise","RQ");
1527 Float_t mnoise = gnoise->GetParameter(1);
1528 cout << "mnoise : " << mnoise << endl;
1529 Float_t rnoise = gnoise->GetParameter(2);
1530 cout << "rnoise : " << rnoise << endl;
1534 //______________________________________________________________________
1535 void AliITSsimulationSDD::WriteSDigits(){
1536 // Fills the Summable digits Tree
1537 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1539 for( Int_t i=0; i<fNofMaps; i++ ) {
1540 if( !fAnodeFire[i] ) continue;
1541 for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1542 Double_t sig = fHitMap2->GetSignal( i, j );
1544 Int_t jdx = j*fScaleSize;
1545 Int_t index = fpList->GetHitIndex( i, j );
1546 AliITSpListItem pItemTmp2( fModule, index, 0. );
1547 // put the fScaleSize analog digits in only one
1548 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1549 AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1550 if( pItemTmp == 0 ) continue;
1551 pItemTmp2.Add( pItemTmp );
1553 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1554 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1555 aliITS->AddSumDigit( pItemTmp2 );
1556 } // end if (sig > 0.2)
1561 //______________________________________________________________________
1562 void AliITSsimulationSDD::PrintStatus() const {
1563 // Print SDD simulation Parameters
1565 cout << "**************************************************" << endl;
1566 cout << " Silicon Drift Detector Simulation Parameters " << endl;
1567 cout << "**************************************************" << endl;
1568 cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1569 cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1570 cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1571 cout << "Number pf Anodes used: " << fNofMaps << endl;
1572 cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1573 cout << "Scale size factor: " << fScaleSize << endl;
1574 cout << "**************************************************" << endl;