ITS/AliITSsimulationSDD.cxx

   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  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 #include <Riostream.h>
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <string.h>
  22
  23 #include <TCanvas.h>
  24 #include <TF1.h>
  25 #include <TH1.h>
  26 #include <TFile.h>
  27 #include <TRandom.h>
  28 #include <TROOT.h>
  29 #include "AliITS.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"
  40 #include "AliLog.h"
  41 #include "AliRun.h"
  42
  43 ClassImp(AliITSsimulationSDD)
  44 ////////////////////////////////////////////////////////////////////////
  45 // Version: 0                                                         //
  46 // Written by Piergiorgio Cerello                                     //
  47 // November 23 1999                                                   //
  48 //                                                                    //
  49 // AliITSsimulationSDD is the simulation of SDDs.                     //
  50 ////////////////////////////////////////////////////////////////////////
  51
  52 //______________________________________________________________________
  53 AliITSsimulationSDD::AliITSsimulationSDD():
  54 AliITSsimulation(),
  55 fITS(0),
  56 fHitMap2(0),
  57 fHitSigMap2(0),
  58 fHitNoiMap2(0),
  59 fElectronics(0),
  60 fInZR(0),
  61 fInZI(0),
  62 fOutZR(0),
  63 fOutZI(0),
  64 fAnodeFire(0),
  65 fHis(0),
  66 fFlag(kFALSE),
  67 fCrosstalkFlag(kFALSE),
  68 fDoFFT(1),
  69 fNofMaps(0),
  70 fMaxNofSamples(0),
  71 fScaleSize(0){
  72     // Default constructor
  73     SetScaleFourier();
  74     SetPerpendTracksFlag();
  75     SetCrosstalkFlag();
  76     SetDoFFT();
  77 }
  78 //______________________________________________________________________
  79 AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
  80     AliITSsimulation(source),
  81 fITS(source.fITS),
  82 fHitMap2(source.fHitMap2),
  83 fHitSigMap2(source.fHitSigMap2),
  84 fHitNoiMap2(source.fHitNoiMap2),
  85 fElectronics(source.fElectronics),
  86 fInZR(source.fInZR),
  87 fInZI(source.fInZI),
  88 fOutZR(source.fOutZR),
  89 fOutZI(source.fOutZI),
  90 fAnodeFire(source.fAnodeFire),
  91 fHis(source.fHis),
  92 fFlag(source.fFlag),
  93 fCrosstalkFlag(source.fCrosstalkFlag),
  94 fDoFFT(source.fDoFFT),
  95 fNofMaps(source.fNofMaps),
  96 fMaxNofSamples(source.fMaxNofSamples),
  97 fScaleSize(source.fScaleSize){
  98     // Copy constructor to satify Coding roules only.
  99
 100 }
 101 //______________________________________________________________________
 102 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
 103     // Assignment operator to satify Coding roules only.
 104
 105     if(this==&src) return *this;
 106     Error("AliITSsimulationSDD","Not allowed to make a = with "
 107           "AliITSsimulationSDD Using default creater instead");
 108     return *this ;
 109 }
 110 /*
 111 //______________________________________________________________________
 112 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
 113     // Assignment operator to satify Coding roules only.
 114
 115     if(this==&src) return *this;
 116     Error("AliITSsimulationSSD","Not allowed to make a = with "
 117           "AliITSsimulationSDD Using default creater instead");
 118     return *this ;
 119 }
 120 */
 121 //______________________________________________________________________
 122 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
 123 AliITSsimulation(dettyp),
 124 fITS(0),
 125 fHitMap2(0),
 126 fHitSigMap2(0),
 127 fHitNoiMap2(0),
 128 fElectronics(0),
 129 fInZR(0),
 130 fInZI(0),
 131 fOutZR(0),
 132 fOutZI(0),
 133 fAnodeFire(0),
 134 fHis(0),
 135 fFlag(kFALSE),
 136 fCrosstalkFlag(kFALSE),
 137 fDoFFT(1),
 138 fNofMaps(0),
 139 fMaxNofSamples(0),
 140 fScaleSize(0){
 141     // Default Constructor
 142   Init();
 143 }
 144 //______________________________________________________________________
 145 void AliITSsimulationSDD::Init(){
 146     // Standard Constructor
 147
 148     SetScaleFourier();
 149     SetPerpendTracksFlag();
 150     SetCrosstalkFlag();
 151     SetDoFFT();
 152
 153     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 154
 155     AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
 156     fpList = new AliITSpList( seg->Npz(),
 157                               fScaleSize*seg->Npx() );
 158     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
 159     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
 160     fHitMap2 = fHitSigMap2;
 161
 162     fNofMaps = seg->Npz();
 163     fMaxNofSamples = seg->Npx();
 164     fAnodeFire = new Bool_t [fNofMaps];
 165
 166     Float_t sddWidth  = seg->Dz();
 167     Float_t anodePitch = seg->Dpz(0);
 168     Double_t timeStep  = (Double_t)seg->Dpx(0);
 169
 170     if(anodePitch*(fNofMaps/2) > sddWidth) {
 171         Warning("AliITSsimulationSDD",
 172                 "Too many anodes %d or too big pitch %f \n",
 173                 fNofMaps/2,anodePitch);
 174     } // end if
 175
 176
 177     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
 178                                     res->Electronics());
 179
 180     char opt1[20], opt2[20];
 181     res->ParamOptions(opt1,opt2);
 182
 183     fITS       = (AliITS*)gAlice->GetModule("ITS");
 184
 185     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
 186     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
 187     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
 188     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];
 189 }
 190 //______________________________________________________________________
 191 AliITSsimulationSDD::~AliITSsimulationSDD() {
 192     // destructor
 193
 194     //    delete fpList;
 195     delete fHitSigMap2;
 196     delete fHitNoiMap2;
 197     delete fElectronics;
 198
 199     fITS = 0;
 200
 201     if (fHis) {
 202         fHis->Delete();
 203         delete fHis;
 204     } // end if fHis
 205     if(fInZR)  delete [] fInZR;
 206     if(fInZI)  delete [] fInZI;
 207     if(fOutZR) delete [] fOutZR;
 208     if(fOutZI) delete [] fOutZI;
 209     if(fAnodeFire) delete [] fAnodeFire;
 210 }
 211 //______________________________________________________________________
 212 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
 213     // create maps to build the lists of tracks for each summable digit
 214     fModule = module;
 215     fEvent  = event;
 216     ClearMaps();
 217     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);
 218 }
 219 //______________________________________________________________________
 220 void AliITSsimulationSDD::ClearMaps() {
 221     // clear maps
 222     fpList->ClearMap();
 223     fHitSigMap2->ClearMap();
 224     fHitNoiMap2->ClearMap();
 225 }
 226 //______________________________________________________________________
 227 void AliITSsimulationSDD::FastFourierTransform(Double_t *real,
 228                           Double_t *imag,Int_t direction) {
 229     // Do a Fast Fourier Transform
 230
 231     Int_t samples = fElectronics->GetSamples();
 232     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
 233     Int_t m1 = samples;
 234     Int_t m  = samples/2;
 235     Int_t m2 = samples/m1;
 236     Int_t i,j,k;
 237     for(i=1; i<=l; i++) {
 238         for(j=0; j<samples; j += m1) {
 239             Int_t p = 0;
 240             for(k=j; k<= j+m-1; k++) {
 241                 Double_t wsr = fElectronics->GetWeightReal(p);
 242                 Double_t wsi = fElectronics->GetWeightImag(p);
 243                 if(direction == -1) wsi = -wsi;
 244                 Double_t xr = *(real+k+m);
 245                 Double_t xi = *(imag+k+m);
 246                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
 247                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
 248                 *(real+k) += xr;
 249                 *(imag+k) += xi;
 250                 p += m2;
 251             } // end for k
 252         } // end for j
 253         m1 = m;
 254         m /= 2;
 255         m2 += m2;
 256     } // end for i
 257
 258     for(j=0; j<samples; j++) {
 259         Int_t j1 = j;
 260         Int_t p = 0;
 261         Int_t i1;
 262         for(i1=1; i1<=l; i1++) {
 263             Int_t j2 = j1;
 264             j1 /= 2;
 265             p = p + p + j2 - j1 - j1;
 266         } // end for i1
 267         if(p >= j) {
 268             Double_t xr = *(real+j);
 269             Double_t xi = *(imag+j);
 270             *(real+j) = *(real+p);
 271             *(imag+j) = *(imag+p);
 272             *(real+p) = xr;
 273             *(imag+p) = xi;
 274         } // end if p>=j
 275     } // end for j
 276     if(direction == -1) {
 277         for(i=0; i<samples; i++) {
 278             *(real+i) /= samples;
 279             *(imag+i) /= samples;
 280         } // end for i
 281     } // end if direction == -1
 282     return;
 283 }
 284
 285 //______________________________________________________________________
 286 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 287     // digitize module using the "slow" detector simulator creating
 288     // summable digits.
 289
 290     TObjArray *fHits = mod->GetHits();
 291     Int_t nhits      = fHits->GetEntriesFast();
 292     if( !nhits ) return;
 293
 294     InitSimulationModule( md, ev );
 295     HitsToAnalogDigits( mod );  // fills fHitMap2 which is = fHitSigmap2
 296     ChargeToSignal( fModule,kFALSE,kTRUE ); // - Process signal adding gain without adding noise
 297     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
 298     ChargeToSignal( fModule,kTRUE,kFALSE );  // - Process only noise
 299     fHitMap2 = fHitSigMap2;   // - Return to signal map
 300     WriteSDigits();
 301     ClearMaps();
 302 }
 303 //______________________________________________________________________
 304 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
 305                                                Int_t mask ) {
 306     // Add Summable digits to module maps.
 307    AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
 308     Int_t    nItems = pItemArray->GetEntries();
 309     Double_t maxadc = res->MaxAdc();
 310     Bool_t sig = kFALSE;
 311
 312     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
 313     for( Int_t i=0; i<nItems; i++ ) {
 314         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
 315         if( pItem->GetModule() != fModule ) {
 316             Error( "AliITSsimulationSDD","Error reading, SDigits module "
 317                    "%d != current module %d: exit",
 318                    pItem->GetModule(), fModule );
 319             return sig;
 320         } // end if
 321
 322         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
 323
 324         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
 325         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
 326         Double_t sigAE = pItem2->GetSignalAfterElect();
 327         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
 328         Int_t ia;
 329         Int_t it;
 330         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
 331         fHitMap2->SetHit( ia, it, sigAE );
 332         fAnodeFire[ia] = kTRUE;
 333     }
 334     return sig;
 335 }
 336 //______________________________________________________________________
 337 void AliITSsimulationSDD::FinishSDigitiseModule() {
 338     // digitize module using the "slow" detector simulator from
 339     // the sum of summable digits.
 340     FinishDigits() ;
 341     ClearMaps();
 342 }
 343 //______________________________________________________________________
 344 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 345     // create maps to build the lists of tracks for each digit
 346
 347     TObjArray *fHits = mod->GetHits();
 348     Int_t nhits      = fHits->GetEntriesFast();
 349
 350     InitSimulationModule( md, ev );
 351     if( !nhits ) return;
 352
 353     HitsToAnalogDigits( mod );
 354     ChargeToSignal( fModule,kTRUE,kTRUE );  // process signal + noise
 355
 356     for( Int_t i=0; i<fNofMaps; i++ ) {
 357         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 358             Int_t jdx = j*fScaleSize;
 359             Int_t index = fpList->GetHitIndex( i, j );
 360             AliITSpListItem pItemTmp2( fModule, index, 0. );
 361             // put the fScaleSize analog digits in only one
 362             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 363                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
 364                 if( pItemTmp == 0 ) continue;
 365                 pItemTmp2.Add( pItemTmp );
 366             }
 367             fpList->DeleteHit( i, j );
 368             fpList->AddItemTo( 0, &pItemTmp2 );
 369         }
 370     }
 371     FinishDigits();
 372     ClearMaps();
 373 }
 374 //______________________________________________________________________
 375 void AliITSsimulationSDD::FinishDigits() {
 376     // introduce the electronics effects and do zero-suppression if required
 377
 378     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
 379
 380     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 381     const char *kopt = res->GetZeroSuppOption();
 382     if (strstr(kopt,"ZS")) Compress2D();
 383     else StoreAllDigits();
 384 }
 385 //______________________________________________________________________
 386 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
 387     // create maps to build the lists of tracks for each digit
 388   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 389   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 390   TObjArray *hits     = mod->GetHits();
 391     Int_t      nhits    = hits->GetEntriesFast();
 392
 393     //    Int_t      arg[6]   = {0,0,0,0,0,0};
 394     Int_t     nofAnodes  = fNofMaps/2;
 395     Double_t  sddLength  = seg->Dx();
 396     Double_t  sddWidth   = seg->Dz();
 397     Double_t  anodePitch = seg->Dpz(0);
 398     Double_t  timeStep   = seg->Dpx(0);
 399     Double_t  driftSpeed ;  // drift velocity (anode dependent)
 400     //Float_t   maxadc     = res->GetMaxAdc();
 401     //Float_t   topValue   = res->GetDynamicRange();
 402     Double_t  norm       = res->GetMaxAdc()/res->GetDynamicRange(); //   maxadc/topValue;
 403     Double_t  cHloss     = res->GetChargeLoss();
 404     Float_t   dfCoeff, s1; res->DiffCoeff(dfCoeff,s1); // Signal 2d Shape
 405     Double_t  eVpairs    = res->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
 406     Double_t  nsigma     = res->GetNSigmaIntegration(); //
 407     Int_t     nlookups   = res->GetGausNLookUp();       //
 408     Float_t   jitter     = res->GetJitterError(); //
 409
 410     // Piergiorgio's part (apart for few variables which I made float
 411     // when i thought that can be done
 412     // Fill detector maps with GEANT hits
 413     // loop over hits in the module
 414
 415     const Float_t kconv = 1.0e+6;  // GeV->KeV
 416     Int_t     itrack      = 0;
 417     Int_t     iWing;       // which detector wing/side.
 418     Int_t     ii,kk,ka,kt; // loop indexs
 419     Int_t     ia,it,index; // sub-pixel integration indexies
 420     Int_t     iAnode;      // anode number.
 421     Int_t     timeSample;  // time buckett.
 422     Int_t     anodeWindow; // anode direction charge integration width
 423     Int_t     timeWindow;  // time direction charge integration width
 424     Int_t     jamin,jamax; // anode charge integration window
 425     Int_t     jtmin,jtmax; // time charge integration window
 426     Int_t     ndiv;        // Anode window division factor.
 427     Int_t     nsplit;      // the number of splits in anode and time windows==1.
 428     Int_t     nOfSplits;   // number of times track length is split into
 429     Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
 430     Float_t   kkF;         // Floating point version of loop index kk.
 431     Double_t  pathInSDD; // Track length in SDD.
 432     Double_t  drPath; // average position of track in detector. in microns
 433     Double_t  drTime; // Drift time
 434     Double_t  nmul;   // drift time window multiplication factor.
 435     Double_t  avDrft;  // x position of path length segment in cm.
 436     Double_t  avAnode; // Anode for path length segment in Anode number (float)
 437     Double_t  zAnode;  // Floating point anode number.
 438     Double_t  driftPath; // avDrft in microns.
 439     Double_t  width;     // width of signal at anodes.
 440     Double_t  depEnergy; // Energy deposited in this GEANT step.
 441     Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
 442     Double_t  sigA; // sigma of signal at anode.
 443     Double_t  sigT; // sigma in time/drift direction for track segment
 444     Double_t  aStep,aConst; // sub-pixel size and offset anode
 445     Double_t  tStep,tConst; // sub-pixel size and offset time
 446     Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
 447     Double_t  chargeloss; // charge loss for track segment.
 448     Double_t  anodeAmplitude; // signal amplitude in anode direction
 449     Double_t  aExpo;          // exponent of Gaussian anode direction
 450     Double_t  timeAmplitude;  // signal amplitude in time direction
 451     Double_t  tExpo;          // exponent of Gaussian time direction
 452     //  Double_t tof;            // Time of flight in ns of this step.
 453
 454     for(ii=0; ii<nhits; ii++) {
 455       if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
 456                               depEnergy,itrack)) continue;
 457       Float_t xloc=xL[0];
 458       if(xloc>0) iWing=0; // left side, carlos channel 0
 459       else iWing=1; // right side
 460
 461       Float_t zloc=xL[2]+0.5*dxL[2];
 462       zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
 463       driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 464       if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
 465         AliWarning("Time Interval > Allowed Time Interval\n");
 466       }
 467       depEnergy  *= kconv;
 468
 469       // scale path to simulate a perpendicular track
 470       // continue if the particle did not lose energy
 471       // passing through detector
 472       if (!depEnergy) {
 473         AliDebug(1,
 474                  Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
 475                       itrack,ii,mod->GetIndex()));
 476         continue;
 477       } // end if !depEnergy
 478
 479       xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
 480       pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
 481
 482       if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
 483       drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
 484       drPath = sddLength-drPath;
 485       if(drPath < 0) {
 486         AliDebug(1, // this should be fixed at geometry level
 487                  Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
 488                    drPath,sddLength,dxL[0],xL[0]));
 489         continue;
 490       } // end if drPath < 0
 491
 492         // Compute number of segments to brake step path into
 493       drTime = drPath/driftSpeed;  //   Drift Time
 494       sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
 495       // calcuate the number of time the path length should be split into.
 496       nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
 497       if(fFlag) nOfSplits = 1;
 498
 499       // loop over path segments, init. some variables.
 500       depEnergy /= nOfSplits;
 501       nOfSplitsF = (Float_t) nOfSplits;
 502       Float_t theAverage=0.,theSteps=0.;
 503       for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
 504         kkF       = (Float_t) kk + 0.5;
 505         avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
 506         avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
 507         theSteps+=1.;
 508         theAverage+=avAnode;
 509         zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
 510         driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 511         driftPath = TMath::Abs(10000.*avDrft);
 512         driftPath = sddLength-driftPath;
 513         if(driftPath < 0) {
 514           AliDebug(1, // this should be fixed at geometry level
 515                    Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
 516                         driftPath,sddLength,avDrft,dxL[0],xL[0]));
 517           continue;
 518         } // end if driftPath < 0
 519         drTime     = driftPath/driftSpeed; // drift time for segment.
 520         timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1); // time bin in range 1-256 !!!
 521         if(timeSample > fScaleSize*fMaxNofSamples) {
 522           AliWarning(Form("Wrong Time Sample: %e",timeSample));
 523           continue;
 524         } // end if timeSample > fScaleSize*fMaxNoofSamples
 525
 526         if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
 527         if(zAnode*anodePitch > sddWidth || zAnode*anodePitch < 0.)
 528           AliWarning(Form("Exceeding sddWidth=%e Z = %e",sddWidth,zAnode*anodePitch));
 529         iAnode = (Int_t) (1.+zAnode); // iAnode in range 1-256 !!!!
 530         if(iAnode < 1 || iAnode > nofAnodes) {
 531           AliWarning(Form("Wrong iAnode: 1<%d>%d  (xanode=%e)",iAnode,nofAnodes, zAnode));
 532           continue;
 533         } // end if iAnode < 1 || iAnode > nofAnodes
 534
 535         // store straight away the particle position in the array
 536         // of particles and take idhit=ii only when part is entering (this
 537         // requires FillModules() in the macro for analysis) :
 538
 539         // Sigma along the anodes for track segment.
 540         sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
 541         sigT       = sigA/driftSpeed;
 542         // Peak amplitude in nanoAmpere
 543         amplitude  = fScaleSize*160.*depEnergy/
 544           (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA);
 545         amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to
 546         // account for clock variations
 547         // (reference value: 40 MHz)
 548         chargeloss = 1.-cHloss*driftPath/1000.;
 549         amplitude *= chargeloss;
 550         width  = 2.*nsigma/(nlookups-1);
 551         // Spread the charge
 552         // Pixel index
 553         ndiv = 2;
 554         nmul = 3.;
 555         if(drTime > 1200.) {
 556           ndiv = 4;
 557           nmul = 1.5;
 558         } // end if drTime > 1200.
 559         // Sub-pixel index
 560         nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2;
 561         // Sub-pixel size see computation of aExpo and tExpo.
 562         aStep  = anodePitch/(nsplit*fScaleSize*sigA);
 563         aConst = zAnode*anodePitch/sigA;
 564         tStep  = timeStep/(nsplit*fScaleSize*sigT);
 565         tConst = drTime/sigT;
 566         // Define SDD window corresponding to the hit
 567         anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1);
 568         timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
 569         jamin = (iAnode - anodeWindow/ndiv - 2)*fScaleSize*nsplit +1;
 570         jamax = (iAnode + anodeWindow/ndiv + 1)*fScaleSize*nsplit;
 571         if(jamin <= 0) jamin = 1;
 572         if(jamax > fScaleSize*nofAnodes*nsplit)
 573           jamax = fScaleSize*nofAnodes*nsplit;
 574         // jtmin and jtmax are Hard-wired
 575         jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1;
 576         jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit;
 577         if(jtmin <= 0) jtmin = 1;
 578         if(jtmax > fScaleSize*fMaxNofSamples*nsplit)
 579           jtmax = fScaleSize*fMaxNofSamples*nsplit;
 580         // Spread the charge in the anode-time window
 581         for(ka=jamin; ka <=jamax; ka++) {
 582           ia = (ka-1)/(fScaleSize*nsplit) + 1;
 583           if(ia <= 0) {
 584             Warning("HitsToAnalogDigits","ia < 1: ");
 585             continue;
 586           } // end if
 587           if(ia > nofAnodes) ia = nofAnodes;
 588           aExpo     = (aStep*(ka-0.5)-aConst);
 589           if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
 590           else {
 591             Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
 592             anodeAmplitude = amplitude*res->GetGausLookUp(theBin);
 593           } // end if TMath::Abs(aEspo) > nsigma
 594           // index starts from 0
 595           index = iWing*nofAnodes+ia-1;
 596           if(anodeAmplitude){
 597             for(kt=jtmin; kt<=jtmax; kt++) {
 598               it = (kt-1)/nsplit+1;  // it starts from 1
 599               if(it<=0){
 600                 Warning("HitsToAnalogDigits","it < 1:");
 601                 continue;
 602               } // end if
 603               if(it>fScaleSize*fMaxNofSamples)
 604                 it = fScaleSize*fMaxNofSamples;
 605               tExpo    = (tStep*(kt-0.5)-tConst);
 606               if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
 607               else {
 608                 Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
 609                 timeAmplitude = anodeAmplitude*res->GetGausLookUp(theBin);
 610               } // end if TMath::Abs(tExpo) > nsigma
 611               // build the list of Sdigits for this module
 612               //                    arg[0]     = index;
 613               //                    arg[1]     = it;
 614               //                    arg[2]     = itrack; // track number
 615               //                    arg[3]     = ii-1; // hit number.
 616               timeAmplitude *= norm;
 617               timeAmplitude *= 10;
 618               //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
 619               Double_t charge = timeAmplitude;
 620               charge += fHitMap2->GetSignal(index,it-1);
 621               fHitMap2->SetHit(index, it-1, charge);
 622               fpList->AddSignal(index,it-1,itrack,ii-1,
 623                                 mod->GetIndex(),timeAmplitude);
 624               fAnodeFire[index] = kTRUE;
 625             }  // end loop over time in window
 626           } // end if anodeAmplitude
 627         } // loop over anodes in window
 628       } // end loop over "sub-hits"
 629     } // end loop over hits
 630 }
 631
 632 //____________________________________________
 633 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
 634   // Adds a Digit.
 635   Int_t size = AliITSdigit::GetNTracks();
 636
 637   Int_t digits[3];
 638   Int_t * tracks = new Int_t[size];
 639   Int_t * hits = new Int_t[size];
 640   Float_t phys;
 641   Float_t * charges = new Float_t[size];
 642
 643   digits[0] = i;
 644   digits[1] = j;
 645   digits[2] = signalc;
 646
 647   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
 648   if( pItem == 0 ) {
 649     phys = 0.0;
 650     for( Int_t l=0; l<size; l++ ) {
 651       tracks[l]  = 0;
 652       hits[l]    = 0;
 653       charges[l] = 0.0;
 654     }
 655   } else {
 656     Int_t idtrack =  pItem->GetTrack( 0 );
 657     if( idtrack >= 0 ) phys = pItem->GetSignal();
 658     else phys = 0.0;
 659
 660     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
 661       tracks[l]  = pItem->GetTrack( l );
 662       hits[l]    = pItem->GetHit( l );
 663       charges[l] = pItem->GetSignal( l );
 664     }else{
 665       tracks[l]  = -3;
 666       hits[l]    = -1;
 667       charges[l] = 0.0;
 668     }// end for if
 669   }
 670
 671   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale );
 672   delete [] tracks;
 673   delete [] hits;
 674   delete [] charges;
 675 }
 676 //______________________________________________________________________
 677 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
 678   // add baseline, noise, gain, electronics and ADC saturation effects
 679   // apply dead channels
 680
 681   char opt1[20], opt2[20];
 682   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 683   res->GetParamOptions(opt1,opt2);
 684   Double_t baseline=0;
 685   Double_t noise=0;
 686   Double_t gain=0;
 687   Float_t contrib=0;
 688   Int_t i,k,kk;
 689   Float_t maxadc = res->GetMaxAdc();
 690
 691   for (i=0;i<fNofMaps;i++) {
 692     if( !fAnodeFire[i] ) continue;
 693     baseline = res->GetBaseline(i);
 694     noise = res->GetNoise(i);
 695     gain = res->GetChannelGain(i);
 696     if(res->IsBad()) gain=0.;
 697     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
 698     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 699       fInZR[k]  = fHitMap2->GetSignal(i,k);
 700       if(bAddGain) fInZR[k]*=gain;
 701       if( bAddNoise ) {
 702         contrib   = (baseline + noise*gRandom->Gaus());
 703         fInZR[k] += contrib;
 704       }
 705       fInZI[k]  = 0.;
 706     } // end for k
 707     if(!fDoFFT) {
 708       for(k=0; k<fMaxNofSamples; k++) {
 709         Double_t newcont = 0.;
 710         Double_t maxcont = 0.;
 711         for(kk=0;kk<fScaleSize;kk++) {
 712           newcont = fInZR[fScaleSize*k+kk];
 713           if(newcont > maxcont) maxcont = newcont;
 714         } // end for kk
 715         newcont = maxcont;
 716         if (newcont >= maxadc) newcont = maxadc -1;
 717         if(newcont >= baseline){
 718           Warning("","newcont=%d>=baseline=%d",newcont,baseline);
 719         } // end if
 720           // back to analog: ?
 721         fHitMap2->SetHit(i,k,newcont);
 722       }  // end for k
 723     }else{
 724       FastFourierTransform(&fInZR[0],&fInZI[0],1);
 725       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 726         Double_t rw = fElectronics->GetTraFunReal(k);
 727         Double_t iw = fElectronics->GetTraFunImag(k);
 728         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 729         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 730       } // end for k
 731       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
 732       for(k=0; k<fMaxNofSamples; k++) {
 733         Double_t newcont1 = 0.;
 734         Double_t maxcont1 = 0.;
 735         for(kk=0;kk<fScaleSize;kk++) {
 736           newcont1 = fOutZR[fScaleSize*k+kk];
 737           if(newcont1 > maxcont1) maxcont1 = newcont1;
 738         } // end for kk
 739         newcont1 = maxcont1;
 740         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 741         fHitMap2->SetHit(i,k,newcont1);
 742       } // end for k
 743     }
 744   } // end for i loop over anodes
 745   return;
 746 }
 747
 748 //______________________________________________________________________
 749 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 750     // function add the crosstalk effect to signal
 751     // temporal function, should be checked...!!!
 752     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 753
 754     Int_t fNofMaps = seg->Npz();
 755     Int_t fMaxNofSamples = seg->Npx();
 756
 757     // create and inizialice crosstalk map
 758     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 759     if( ctk == NULL ) {
 760         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 761         return;
 762     }
 763     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 764     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 765     for( Int_t z=0; z<fNofMaps; z++ ) {
 766       Double_t baseline = calibr->GetBaseline(z);
 767         Bool_t on = kFALSE;
 768         Int_t tstart = 0;
 769         Int_t tstop = 0;
 770         Int_t nTsteps = 0;
 771
 772         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 773             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 774             if( fadc > baseline ) {
 775                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 776                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 777                     if( fadc1 < fadc ) continue;
 778                     on = kTRUE;
 779                     nTsteps = 0;
 780                     tstart = l;
 781                 }
 782                 nTsteps++;
 783             }
 784             else { // end fadc > baseline
 785                 if( on == kTRUE ) {
 786                     if( nTsteps > 2 ) {
 787                         tstop = l;
 788                         // make smooth derivative
 789                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 790                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 791                         if( ctk == NULL ) {
 792                             Error( "ApplyCrosstalk",
 793                                    "no memory for temporal array: exit \n" );
 794                             return;
 795                         }
 796                         for( Int_t i=tstart; i<tstop; i++ ) {
 797                             if( i > 2 && i < fMaxNofSamples-2 )
 798                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 799                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 800                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 801                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 802                         }
 803
 804                         // add crosstalk contribution to neibourg anodes
 805                         for( Int_t i=tstart; i<tstop; i++ ) {
 806                             Int_t anode = z - 1;
 807                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 808                             Float_t ctktmp =  -dev[i1] * 0.25;
 809                             if( anode > 0 ) {
 810                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 811                             }
 812                             anode = z + 1;
 813                             if( anode < fNofMaps ) {
 814                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 815                             }
 816                         }
 817                         delete [] dev;
 818
 819                     } // if( nTsteps > 2 )
 820                     on = kFALSE;
 821                 }  // if( on == kTRUE )
 822             }  // else
 823         }
 824     }
 825
 826     for( Int_t a=0; a<fNofMaps; a++ )
 827         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 828             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 829             fHitMap2->SetHit( a, t, signal );
 830         }
 831
 832     delete [] ctk;
 833 }
 834
 835 //______________________________________________________________________
 836 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 837     // To the 10 to 8 bit lossive compression.
 838     // code from Davide C. and Albert W.
 839
 840     if (signal < 128)  return signal;
 841     if (signal < 256)  return (128+((signal-128)>>1));
 842     if (signal < 512)  return (192+((signal-256)>>3));
 843     if (signal < 1024) return (224+((signal-512)>>4));
 844     return 0;
 845 }
 846 //______________________________________________________________________
 847 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
 848   // Decompression from 8 to 10 bit
 849
 850   if (signal < 0 || signal > 255) {
 851     AliWarning(Form("Signal value %d out of range",signal));
 852     return 0;
 853   } // end if signal <0 || signal >255
 854
 855   if (signal < 128) return signal;
 856   if (signal < 192) {
 857     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
 858     else  return (128+((signal-128)<<1)+1);
 859   } // end if signal < 192
 860   if (signal < 224) {
 861     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
 862     else  return (256+((signal-192)<<3)+4);
 863   } // end if signal < 224
 864   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
 865   return (512+((signal-224)<<4)+8);
 866 }
 867 //______________________________________________________________________
 868 void AliITSsimulationSDD::Compress2D(){
 869   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
 870   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 871   for (Int_t iWing=0; iWing<2; iWing++) {
 872     Int_t tL=res->GetZSLowThreshold(iWing);
 873     Int_t tH=res->GetZSHighThreshold(iWing);
 874     for (Int_t i=0; i<fNofMaps/2; i++) {
 875       Int_t ian=i+iWing*fNofMaps/2;
 876       if( !fAnodeFire[ian] ) continue;
 877       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
 878         Int_t nLow=0, nHigh=0;
 879         Float_t cC=fHitMap2->GetSignal(ian,itb);
 880         if(cC<=tL) continue;
 881         nLow++; // cC is greater than tL
 882         if(cC>tH) nHigh++;
 883         //                     N
 884         // Get "quintuple":   WCE
 885         //                     S
 886         Float_t wW=0.;
 887         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
 888         if(wW>tL) nLow++;
 889         if(wW>tH) nHigh++;
 890         Float_t eE=0.;
 891         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
 892         if(eE>tL) nLow++;
 893         if(eE>tH) nHigh++;
 894         Float_t nN=0.;
 895         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
 896         if(nN>tL) nLow++;
 897         if(nN>tH) nHigh++;
 898         Float_t sS=0.;
 899         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
 900         if(sS>tL) nLow++;
 901         if(sS>tH) nHigh++;
 902
 903         if(nLow>=3 && nHigh>=1){
 904           Int_t signal=(Int_t)cC;
 905           Int_t signalc = Convert10to8(signal);
 906           Int_t signale = Convert8to10(signalc);
 907           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
 908           AddDigit(ian,itb,signalc,signale);  // store C
 909         }
 910       }
 911     }
 912   }
 913 }
 914
 915
 916 //______________________________________________________________________
 917 void AliITSsimulationSDD::StoreAllDigits(){
 918   // if non-zero-suppressed data
 919   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 920
 921   Int_t i, j, digits[3];
 922
 923   for (i=0; i<fNofMaps; i++) {
 924     for (j=0; j<fMaxNofSamples; j++) {
 925       Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
 926       signal = Convert10to8(signal);
 927       signal = Convert8to10(signal);
 928       digits[0] = i;
 929       digits[1] = j;
 930       digits[2] = signal;
 931       fITS->AddRealDigit(1,digits);
 932     } // end for j
 933   } // end for i
 934 }
 935 //______________________________________________________________________
 936 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
 937     // Creates histograms of maps for debugging
 938     Int_t i;
 939
 940     fHis=new TObjArray(fNofMaps);
 941     for (i=0;i<fNofMaps;i++) {
 942         TString sddName("sdd_");
 943         Char_t candNum[4];
 944         sprintf(candNum,"%d",i+1);
 945         sddName.Append(candNum);
 946         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
 947                              0.,(Float_t) scale*fMaxNofSamples), i);
 948     } // end for i
 949 }
 950 //______________________________________________________________________
 951 void AliITSsimulationSDD::FillHistograms(){
 952     // fill 1D histograms from map
 953
 954     if (!fHis) return;
 955
 956     for( Int_t i=0; i<fNofMaps; i++) {
 957         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
 958         Int_t nsamples = hist->GetNbinsX();
 959         for( Int_t j=0; j<nsamples; j++) {
 960             Double_t signal=fHitMap2->GetSignal(i,j);
 961             hist->Fill((Float_t)j,signal);
 962         } // end for j
 963     } // end for i
 964 }
 965 //______________________________________________________________________
 966 void AliITSsimulationSDD::ResetHistograms(){
 967     // Reset histograms for this detector
 968     Int_t i;
 969
 970     for (i=0;i<fNofMaps;i++ ) {
 971         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 972     } // end for i
 973 }
 974 //______________________________________________________________________
 975 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
 976     // Fills a histogram from a give anode.
 977
 978     if (!fHis) return 0;
 979
 980     if(wing <=0 || wing > 2) {
 981         Warning("GetAnode","Wrong wing number: %d",wing);
 982         return NULL;
 983     } // end if wing <=0 || wing >2
 984     if(anode <=0 || anode > fNofMaps/2) {
 985         Warning("GetAnode","Wrong anode number: %d",anode);
 986         return NULL;
 987     } // end if ampde <=0 || andoe > fNofMaps/2
 988
 989     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
 990     return (TH1F*)(fHis->At(index));
 991 }
 992 //______________________________________________________________________
 993 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
 994     // Writes the histograms to a file
 995
 996     if (!fHis) return;
 997
 998     hfile->cd();
 999     Int_t i;
1000     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
1001     return;
1002 }
1003 //______________________________________________________________________
1004 void AliITSsimulationSDD::WriteSDigits(){
1005     // Fills the Summable digits Tree
1006     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1007
1008     for( Int_t i=0; i<fNofMaps; i++ ) {
1009         if( !fAnodeFire[i] ) continue;
1010         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1011             Double_t sig = fHitMap2->GetSignal( i, j );
1012             if( sig > 0.2 ) {
1013                 Int_t jdx = j*fScaleSize;
1014                 Int_t index = fpList->GetHitIndex( i, j );
1015                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1016                 // put the fScaleSize analog digits in only one
1017                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1018                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1019                     if( pItemTmp == 0 ) continue;
1020                     pItemTmp2.Add( pItemTmp );
1021                 }
1022                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1023                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1024                 aliITS->AddSumDigit( pItemTmp2 );
1025             } // end if (sig > 0.2)
1026         }
1027     }
1028     return;
1029 }
1030 //______________________________________________________________________
1031 void AliITSsimulationSDD::PrintStatus() const {
1032     // Print SDD simulation Parameters
1033
1034     cout << "**************************************************" << endl;
1035     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1036     cout << "**************************************************" << endl;
1037     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1038     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1039     cout << "Number of Anodes used: " << fNofMaps << endl;
1040     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1041     cout << "Scale size factor: " << fScaleSize << endl;
1042     cout << "**************************************************" << endl;
1043 }