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