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->IsBad()) gain=0.;
 752     if( res->IsChipBad(res->GetChip(i)) ){
 753       printf("Chip bad mod %d  chip %d anode %d\n",mod,res->GetChip(i),i);
 754       gain=0.;
 755     }
 756     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 757       fInZR[k]  = fHitMap2->GetSignal(i,k);
 758       if(bAddGain) fInZR[k]*=gain;
 759       if( bAddNoise ) {
 760         contrib   = (baseline + noise*gRandom->Gaus());
 761         fInZR[k] += contrib;
 762       }
 763       fInZI[k]  = 0.;
 764     } // end for k
 765     if(!fDoFFT) {
 766       for(k=0; k<fMaxNofSamples; k++) {
 767         Double_t newcont = 0.;
 768         Double_t maxcont = 0.;
 769         for(kk=0;kk<fScaleSize;kk++) {
 770           newcont = fInZR[fScaleSize*k+kk];
 771           if(newcont > maxcont) maxcont = newcont;
 772         } // end for kk
 773         newcont = maxcont;
 774         if (newcont >= maxadc) newcont = maxadc -1;
 775         if(newcont >= baseline){
 776           Warning("","newcont=%d>=baseline=%d",newcont,baseline);
 777         } // end if
 778           // back to analog: ?
 779         fHitMap2->SetHit(i,k,newcont);
 780       }  // end for k
 781     }else{
 782       FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
 783       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 784         Double_t rw = fElectronics->GetTraFunReal(k);
 785         Double_t iw = fElectronics->GetTraFunImag(k);
 786         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 787         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 788       } // end for k
 789       FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
 790       for(k=0; k<fMaxNofSamples; k++) {
 791         Double_t newcont1 = 0.;
 792         Double_t maxcont1 = 0.;
 793         for(kk=0;kk<fScaleSize;kk++) {
 794           newcont1 = fOutZR[fScaleSize*k+kk];
 795           if(newcont1 > maxcont1) maxcont1 = newcont1;
 796         } // end for kk
 797         newcont1 = maxcont1;
 798         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 799         fHitMap2->SetHit(i,k,newcont1);
 800       } // end for k
 801     }
 802   } // end for i loop over anodes
 803   return;
 804 }
 805
 806 //______________________________________________________________________
 807 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 808     // function add the crosstalk effect to signal
 809     // temporal function, should be checked...!!!
 810     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 811
 812     Int_t fNofMaps = seg->Npz();
 813     Int_t fMaxNofSamples = seg->Npx();
 814
 815     // create and inizialice crosstalk map
 816     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 817     if( ctk == NULL ) {
 818         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 819         return;
 820     }
 821     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 822     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 823     for( Int_t z=0; z<fNofMaps; z++ ) {
 824       Double_t baseline = calibr->GetBaseline(z);
 825         Bool_t on = kFALSE;
 826         Int_t tstart = 0;
 827         Int_t tstop = 0;
 828         Int_t nTsteps = 0;
 829
 830         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 831             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 832             if( fadc > baseline ) {
 833                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 834                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 835                     if( fadc1 < fadc ) continue;
 836                     on = kTRUE;
 837                     nTsteps = 0;
 838                     tstart = l;
 839                 }
 840                 nTsteps++;
 841             }
 842             else { // end fadc > baseline
 843                 if( on == kTRUE ) {
 844                     if( nTsteps > 2 ) {
 845                         tstop = l;
 846                         // make smooth derivative
 847                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 848                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 849                         if( ctk == NULL ) {
 850                             Error( "ApplyCrosstalk",
 851                                    "no memory for temporal array: exit \n" );
 852                             return;
 853                         }
 854                         for( Int_t i=tstart; i<tstop; i++ ) {
 855                             if( i > 2 && i < fMaxNofSamples-2 )
 856                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 857                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 858                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 859                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 860                         }
 861
 862                         // add crosstalk contribution to neibourg anodes
 863                         for( Int_t i=tstart; i<tstop; i++ ) {
 864                             Int_t anode = z - 1;
 865                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 866                             Float_t ctktmp =  -dev[i1] * 0.25;
 867                             if( anode > 0 ) {
 868                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 869                             }
 870                             anode = z + 1;
 871                             if( anode < fNofMaps ) {
 872                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 873                             }
 874                         }
 875                         delete [] dev;
 876
 877                     } // if( nTsteps > 2 )
 878                     on = kFALSE;
 879                 }  // if( on == kTRUE )
 880             }  // else
 881         }
 882     }
 883
 884     for( Int_t a=0; a<fNofMaps; a++ )
 885         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 886             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 887             fHitMap2->SetHit( a, t, signal );
 888         }
 889
 890     delete [] ctk;
 891 }
 892
 893 //______________________________________________________________________
 894 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
 895                                            Int_t &th) const{
 896     // Returns the compression alogirthm parameters
 897   db=fD[i];
 898   tl=fT1[i];
 899   th=fT2[i];
 900 }
 901 //______________________________________________________________________
 902 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{
 903     // returns the compression alogirthm parameters
 904
 905     db=fD[i];
 906     tl=fT1[i];
 907
 908 }
 909 //______________________________________________________________________
 910 void AliITSsimulationSDD::SetCompressParam(){
 911     // Sets the compression alogirthm parameters
 912     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 913     for(Int_t ian = 0; ian<fNofMaps;ian++){
 914       fD[ian] = (Int_t)(calibr->GetBaseline(ian));
 915       fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5);
 916       fT2[ian] = 0;   // used by 2D clustering - not defined yet
 917       fTol[ian] = 0; // used by 2D clustering - not defined yet
 918     }
 919 }
 920 //______________________________________________________________________
 921 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 922     // To the 10 to 8 bit lossive compression.
 923     // code from Davide C. and Albert W.
 924
 925     if (signal < 128)  return signal;
 926     if (signal < 256)  return (128+((signal-128)>>1));
 927     if (signal < 512)  return (192+((signal-256)>>3));
 928     if (signal < 1024) return (224+((signal-512)>>4));
 929     return 0;
 930 }
 931 /*
 932 //______________________________________________________________________
 933 AliITSMap*   AliITSsimulationSDD::HitMap(Int_t i){
 934     //Return the correct map.
 935
 936     return ((i==0)? fHitMap1 : fHitMap2);
 937 }
 938 */
 939 //______________________________________________________________________
 940 void AliITSsimulationSDD::ZeroSuppression(const char *option) {
 941     // perform the zero suppresion
 942     if (strstr(option,"2D")) {
 943         //Init2D();              // activate if param change module by module
 944         Compress2D();
 945     } else if (strstr(option,"1D")) {
 946         //Init1D();              // activate if param change module by module
 947         Compress1D();
 948     } else StoreAllDigits();
 949 }
 950 //______________________________________________________________________
 951 void AliITSsimulationSDD::Init2D(){
 952     // read in and prepare arrays: fD, fT1, fT2
 953     //                         savemu[nanodes], savesigma[nanodes]
 954     // read baseline and noise from file - either a .root file and in this
 955     // case data should be organised in a tree with one entry for each
 956     // module => reading should be done accordingly
 957     // or a classic file and do smth. like this ( code from Davide C. and
 958     // Albert W.) :
 959     // Read 2D zero-suppression parameters for SDD
 960
 961     if (!strstr(fParam.Data(),"file")) return;
 962
 963     Int_t na,pos,tempTh;
 964     Float_t mu,sigma;
 965     Float_t *savemu    = new Float_t [fNofMaps];
 966     Float_t *savesigma = new Float_t [fNofMaps];
 967     char input[100],basel[100],par[100];
 968     char *filtmp;
 969     Double_t tmp1,tmp2;
 970     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 971
 972     res->Thresholds(tmp1,tmp2);
 973     Int_t minval = static_cast<Int_t>(tmp1);
 974
 975     res->Filenames(input,basel,par);
 976     fFileName = par;
 977     //
 978     filtmp = gSystem->ExpandPathName(fFileName.Data());
 979     FILE *param = fopen(filtmp,"r");
 980     na = 0;
 981
 982     if(param) {
 983         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
 984             if (pos != na+1) {
 985                 Error("Init2D","Anode number not in increasing order!",filtmp);
 986                 exit(1);
 987             } // end if pos != na+1
 988             savemu[na] = mu;
 989             savesigma[na] = sigma;
 990             if ((2.*sigma) < mu) {
 991                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
 992                 mu = 2.0 * sigma;
 993             } else fD[na] = 0;
 994             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
 995             if (tempTh < 0) tempTh=0;
 996             fT1[na] = tempTh;
 997             tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
 998             if (tempTh < 0) tempTh=0;
 999             fT2[na] = tempTh;
1000             na++;
1001         } // end while
1002     } else {
1003         Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1004         exit(1);
1005     } // end if(param)
1006
1007     fclose(param);
1008     delete [] filtmp;
1009     delete [] savemu;
1010     delete [] savesigma;
1011 }
1012 //______________________________________________________________________
1013 void AliITSsimulationSDD::Compress2D(){
1014     // simple ITS cluster finder -- online zero-suppression conditions
1015
1016     Int_t db,tl,th;
1017     Double_t tmp1,tmp2;
1018     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1019
1020     res->Thresholds(tmp1,tmp2);
1021     Int_t minval   = static_cast<Int_t>(tmp1);
1022     Bool_t write   = res->OutputOption();
1023     Bool_t do10to8 = res->Do10to8();
1024     Int_t nz, nl, nh, low, i, j;
1025     SetCompressParam();
1026     for (i=0; i<fNofMaps; i++) {
1027         CompressionParam(i,db,tl,th);
1028         nz  = 0;
1029         nl  = 0;
1030         nh  = 0;
1031         low = 0;
1032         for (j=0; j<fMaxNofSamples; j++) {
1033             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1034             signal -= db; // if baseline eq. is done here
1035             if (signal <= 0) {nz++; continue;}
1036             if ((signal - tl) < minval) low++;
1037             if ((signal - th) >= minval) {
1038                 nh++;
1039                 Bool_t cond=kTRUE;
1040                 FindCluster(i,j,signal,minval,cond);
1041                 if(cond && j &&
1042                    ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1043                     if(do10to8) signal = Convert10to8(signal);
1044                     AddDigit(i,j,signal);
1045                 } // end if cond&&j&&()
1046             } else if ((signal - tl) >= minval) nl++;
1047         } // end for j loop time samples
1048         if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1049     } //end for i loop anodes
1050
1051     char hname[30];
1052     if (write) {
1053         sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1054         TreeB()->Write(hname);
1055         // reset tree
1056         TreeB()->Reset();
1057     } // end if write
1058 }
1059 //______________________________________________________________________
1060 void  AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
1061                                        Int_t minval,Bool_t &cond){
1062     // Find clusters according to the online 2D zero-suppression algorithm
1063     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1064     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1065
1066     Bool_t do10to8 = res->Do10to8();
1067     Bool_t high    = kFALSE;
1068
1069     fHitMap2->FlagHit(i,j);
1070     //
1071     //  check the online zero-suppression conditions
1072     //
1073     const Int_t kMaxNeighbours = 4;
1074     Int_t nn;
1075     Int_t dbx,tlx,thx;
1076     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1077     seg->Neighbours(i,j,&nn,xList,yList);
1078     Int_t in,ix,iy,qns;
1079     for (in=0; in<nn; in++) {
1080         ix=xList[in];
1081         iy=yList[in];
1082         if (fHitMap2->TestHit(ix,iy)==kUnused) {
1083             CompressionParam(ix,dbx,tlx,thx);
1084             Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1085             qn -= dbx; // if baseline eq. is done here
1086             if ((qn-tlx) < minval) {
1087                 fHitMap2->FlagHit(ix,iy);
1088                 continue;
1089             } else {
1090                 if ((qn - thx) >= minval) high=kTRUE;
1091                 if (cond) {
1092                     if(do10to8) signal = Convert10to8(signal);
1093                     AddDigit(i,j,signal);
1094                 } // end if cond
1095                 if(do10to8) qns = Convert10to8(qn);
1096                 else qns=qn;
1097                 if (!high) AddDigit(ix,iy,qns);
1098                 cond=kFALSE;
1099                 if(!high) fHitMap2->FlagHit(ix,iy);
1100             } // end if qn-tlx < minval
1101         } // end if  TestHit
1102     } // end for in loop over neighbours
1103 }
1104 //______________________________________________________________________
1105 void AliITSsimulationSDD::Init1D(){
1106     // this is just a copy-paste of input taken from 2D algo
1107     // Torino people should give input
1108     // Read 1D zero-suppression parameters for SDD
1109
1110     if (!strstr(fParam.Data(),"file")) return;
1111
1112     Int_t na,pos,tempTh;
1113     Float_t mu,sigma;
1114     Float_t *savemu    = new Float_t [fNofMaps];
1115     Float_t *savesigma = new Float_t [fNofMaps];
1116     char input[100],basel[100],par[100];
1117     char *filtmp;
1118     Double_t tmp1,tmp2;
1119     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1120
1121     res->Thresholds(tmp1,tmp2);
1122     Int_t minval = static_cast<Int_t>(tmp1);
1123
1124     res->Filenames(input,basel,par);
1125     fFileName=par;
1126
1127     //  set first the disable and tol param
1128     SetCompressParam();
1129     //
1130     filtmp = gSystem->ExpandPathName(fFileName.Data());
1131     FILE *param = fopen(filtmp,"r");
1132     na = 0;
1133
1134     if (param) {
1135         fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1136         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1137             if (pos != na+1) {
1138                 Error("Init1D","Anode number not in increasing order!",filtmp);
1139                 exit(1);
1140             } // end if pos != na+1
1141             savemu[na]=mu;
1142             savesigma[na]=sigma;
1143             if ((2.*sigma) < mu) {
1144                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1145                 mu = 2.0 * sigma;
1146             } else fD[na] = 0;
1147             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1148             if (tempTh < 0) tempTh=0;
1149             fT1[na] = tempTh;
1150             na++;
1151         } // end while
1152     } else {
1153         Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1154         exit(1);
1155     } // end if(param)
1156
1157     fclose(param);
1158     delete [] filtmp;
1159     delete [] savemu;
1160     delete [] savesigma;
1161 }
1162 //______________________________________________________________________
1163 void AliITSsimulationSDD::Compress1D(){
1164     // 1D zero-suppression algorithm (from Gianluca A.)
1165     Int_t    dis,tol,thres,decr,diff;
1166     UChar_t *str=fStream->Stream();
1167     Int_t    counter=0;
1168     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1169
1170     Bool_t   do10to8=res->Do10to8();
1171     Int_t    last=0;
1172     Int_t    k,i,j;
1173     SetCompressParam();
1174     for (k=0; k<2; k++) {
1175         tol = Tolerance(k);
1176         dis = Disable(k);
1177         for (i=0; i<fNofMaps/2; i++) {
1178             Bool_t firstSignal=kTRUE;
1179             Int_t idx=i+k*fNofMaps/2;
1180             if( !fAnodeFire[idx] ) continue;
1181             CompressionParam(idx,decr,thres);
1182
1183             for (j=0; j<fMaxNofSamples; j++) {
1184                 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1185                 signal -= decr;  // if baseline eq.
1186                 if(do10to8) signal = Convert10to8(signal);
1187                 if (signal <= thres) {
1188                     signal=0;
1189                     diff=128;
1190                     last=0;
1191                     // write diff in the buffer for HuffT
1192                     str[counter]=(UChar_t)diff;
1193                     counter++;
1194                     continue;
1195                 } // end if signal <= thres
1196                 diff=signal-last;
1197                 if (diff > 127) diff=127;
1198                 if (diff < -128) diff=-128;
1199                 if (signal < dis) {
1200                     // tol has changed to 8 possible cases ? - one can write
1201                     // this if(TMath::Abs(diff)<tol) ... else ...
1202                     if(TMath::Abs(diff)<tol) diff=0;
1203                     // or keep it as it was before
1204                     AddDigit(idx,j,last+diff);
1205                 } else {
1206                     AddDigit(idx,j,signal);
1207                 } // end if singal < dis
1208                 diff += 128;
1209                 // write diff in the buffer used to compute Huffman tables
1210                 if (firstSignal) str[counter]=(UChar_t)signal;
1211                 else str[counter]=(UChar_t)diff;
1212                 counter++;
1213                 last=signal;
1214                 firstSignal=kFALSE;
1215             } // end for j loop time samples
1216         } // end for i loop anodes  one half of detector
1217     } //  end for k
1218
1219     // check
1220     fStream->CheckCount(counter);
1221
1222     // open file and write out the stream of diff's
1223     static Bool_t open=kTRUE;
1224     static TFile *outFile;
1225     Bool_t write = res->OutputOption();
1226     TDirectory *savedir = gDirectory;
1227
1228     if (write ) {
1229         if(open) {
1230             SetFileName("stream.root");
1231             cout<<"filename "<<fFileName<<endl;
1232             outFile=new TFile(fFileName,"recreate");
1233             cout<<"I have opened "<<fFileName<<" file "<<endl;
1234         } // end if open
1235         open = kFALSE;
1236         outFile->cd();
1237         fStream->Write();
1238     }  // endif write
1239
1240     fStream->ClearStream();
1241
1242     // back to galice.root file
1243     if(savedir) savedir->cd();
1244 }
1245 //______________________________________________________________________
1246 void AliITSsimulationSDD::StoreAllDigits(){
1247     // if non-zero-suppressed data
1248     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1249
1250     Bool_t do10to8 = res->Do10to8();
1251     Int_t i, j, digits[3];
1252
1253     for (i=0; i<fNofMaps; i++) {
1254         for (j=0; j<fMaxNofSamples; j++) {
1255             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1256             if(do10to8) signal = Convert10to8(signal);
1257             digits[0] = i;
1258             digits[1] = j;
1259             digits[2] = signal;
1260             fITS->AddRealDigit(1,digits);
1261         } // end for j
1262     } // end for i
1263 }
1264 //______________________________________________________________________
1265 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
1266     // Creates histograms of maps for debugging
1267     Int_t i;
1268
1269     fHis=new TObjArray(fNofMaps);
1270     for (i=0;i<fNofMaps;i++) {
1271         TString sddName("sdd_");
1272         Char_t candNum[4];
1273         sprintf(candNum,"%d",i+1);
1274         sddName.Append(candNum);
1275         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1276                              0.,(Float_t) scale*fMaxNofSamples), i);
1277     } // end for i
1278 }
1279 //______________________________________________________________________
1280 void AliITSsimulationSDD::FillHistograms(){
1281     // fill 1D histograms from map
1282
1283     if (!fHis) return;
1284
1285     for( Int_t i=0; i<fNofMaps; i++) {
1286         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1287         Int_t nsamples = hist->GetNbinsX();
1288         for( Int_t j=0; j<nsamples; j++) {
1289             Double_t signal=fHitMap2->GetSignal(i,j);
1290             hist->Fill((Float_t)j,signal);
1291         } // end for j
1292     } // end for i
1293 }
1294 //______________________________________________________________________
1295 void AliITSsimulationSDD::ResetHistograms(){
1296     // Reset histograms for this detector
1297     Int_t i;
1298
1299     for (i=0;i<fNofMaps;i++ ) {
1300         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
1301     } // end for i
1302 }
1303 //______________________________________________________________________
1304 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
1305     // Fills a histogram from a give anode.
1306
1307     if (!fHis) return 0;
1308
1309     if(wing <=0 || wing > 2) {
1310         Warning("GetAnode","Wrong wing number: %d",wing);
1311         return NULL;
1312     } // end if wing <=0 || wing >2
1313     if(anode <=0 || anode > fNofMaps/2) {
1314         Warning("GetAnode","Wrong anode number: %d",anode);
1315         return NULL;
1316     } // end if ampde <=0 || andoe > fNofMaps/2
1317
1318     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1319     return (TH1F*)(fHis->At(index));
1320 }
1321 //______________________________________________________________________
1322 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
1323     // Writes the histograms to a file
1324
1325     if (!fHis) return;
1326
1327     hfile->cd();
1328     Int_t i;
1329     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
1330     return;
1331 }
1332 //______________________________________________________________________
1333 Float_t AliITSsimulationSDD::GetNoise() {
1334     // Returns the noise value
1335     //Bool_t do10to8=GetResp()->Do10to8();
1336     //noise will always be in the liniar part of the signal
1337     Int_t decr;
1338     Int_t threshold = fT1[0];
1339     char opt1[20], opt2[20];
1340     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1341     SetCompressParam();
1342     res->GetParamOptions(opt1,opt2);
1343     fParam=opt2;
1344     Double_t noise,baseline;
1345     //GetBaseline(fModule);
1346
1347     TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1348     if(c2) delete c2->GetPrimitive("noisehist");
1349     if(c2) delete c2->GetPrimitive("anode");
1350     else     c2=new TCanvas("c2");
1351     c2->cd();
1352     c2->SetFillColor(0);
1353
1354     TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1355     TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1356                            (float)fMaxNofSamples);
1357     Int_t i,k;
1358     for (i=0;i<fNofMaps;i++) {
1359         CompressionParam(i,decr,threshold);
1360         baseline = res->GetBaseline(i);
1361         noise = res->GetNoise(i);
1362         anode->Reset();
1363         for (k=0;k<fMaxNofSamples;k++) {
1364             Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1365             //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1366             if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1367             anode->Fill((float)k,signal);
1368         } // end for k
1369         anode->Draw();
1370         c2->Update();
1371     } // end for i
1372     TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1373     noisehist->Fit("gnoise","RQ");
1374     noisehist->Draw();
1375     c2->Update();
1376     Float_t mnoise = gnoise->GetParameter(1);
1377     cout << "mnoise : " << mnoise << endl;
1378     Float_t rnoise = gnoise->GetParameter(2);
1379     cout << "rnoise : " << rnoise << endl;
1380     delete noisehist;
1381     return rnoise;
1382 }
1383 //______________________________________________________________________
1384 void AliITSsimulationSDD::WriteSDigits(){
1385     // Fills the Summable digits Tree
1386     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1387
1388     for( Int_t i=0; i<fNofMaps; i++ ) {
1389         if( !fAnodeFire[i] ) continue;
1390         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1391             Double_t sig = fHitMap2->GetSignal( i, j );
1392             if( sig > 0.2 ) {
1393                 Int_t jdx = j*fScaleSize;
1394                 Int_t index = fpList->GetHitIndex( i, j );
1395                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1396                 // put the fScaleSize analog digits in only one
1397                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1398                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1399                     if( pItemTmp == 0 ) continue;
1400                     pItemTmp2.Add( pItemTmp );
1401                 }
1402                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1403                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1404                 aliITS->AddSumDigit( pItemTmp2 );
1405             } // end if (sig > 0.2)
1406         }
1407     }
1408     return;
1409 }
1410 //______________________________________________________________________
1411 void AliITSsimulationSDD::PrintStatus() const {
1412     // Print SDD simulation Parameters
1413
1414     cout << "**************************************************" << endl;
1415     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1416     cout << "**************************************************" << endl;
1417     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1418     cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1419     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1420     cout << "Number pf Anodes used: " << fNofMaps << endl;
1421     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1422     cout << "Scale size factor: " << fScaleSize << endl;
1423     cout << "**************************************************" << endl;
1424 }