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