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