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