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