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