ITS/AliITSsimulationSDD.cxx

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