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