ITS/AliITSsimulationSDD.cxx

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