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