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