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