ITS/AliITSsimulationSDD.cxx

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