ITS/AliITSsimulationSSD.cxx

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <iostream.h>
   4 #include <TObjArray.h>
   5 #include <TRandom.h>
   6 #include <TMath.h>
   7
   8 #include "AliITSmodule.h"
   9 #include "AliITSMapA2.h"
  10 #include "AliITSsegmentationSSD.h"
  11 #include "AliITSresponseSSD.h"
  12 #include "AliITSsimulationSSD.h"
  13 //#include "AliITSdictSSD.h"
  14 #include "AliITSdcsSSD.h"
  15 #include "AliITS.h"
  16 #include "AliRun.h"
  17
  18
  19 ClassImp(AliITSsimulationSSD);
  20 ////////////////////////////////////////////////////////////////////////
  21 // Version: 0
  22 // Written by Enrico Fragiacomo
  23 // July 2000
  24 //
  25 // AliITSsimulationSSD is the simulation of SSDs.
  26
  27 //------------------------------------------------------------
  28 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
  29                                          AliITSresponse *resp){
  30   // Constructor
  31
  32     fSegmentation = seg;
  33     fResponse = resp;
  34   Float_t noise[2] = {0.,0.};
  35   fResponse->GetNoiseParam(noise[0],noise[1]);   // retrieves noise parameters
  36   cout<<"nois1,2 ="<<noise[0]<<","<<noise[1]<<endl;
  37   fDCS = new AliITSdcsSSD(seg,resp);
  38
  39     fNstrips = fSegmentation->Npx();
  40     fPitch = fSegmentation->Dpx(0);
  41     cout<<" Dx,Dz ="<<fSegmentation->Dx()<<","<<fSegmentation->Dz()<<endl;
  42     cout<<"fNstrips="<<fNstrips<<" fPitch="<<fPitch<<endl;
  43
  44     //fP = new TArrayF(fNstrips+1);
  45     //fN = new TArrayF(fNstrips+1);
  46
  47     fMapA2 = new AliITSMapA2(fSegmentation);
  48
  49     //fTracksP = new AliITSdictSSD[fNstrips+1];
  50     //fTracksN = new AliITSdictSSD[fNstrips+1];
  51
  52     fSteps  = 100;   // still hard-wired - set in SetDetParam and get it via
  53                      // fDCS together with the others eventually
  54 }
  55
  56 //___________________________________________________________________________
  57 AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD
  58                                                                       &source){
  59 // Operator =
  60
  61     if(this==&source) return *this;
  62
  63     this->fDCS = new AliITSdcsSSD(*(source.fDCS));
  64     //this->fN   = new TArrayF(*(source.fN));
  65     //this->fP   = new TArrayF(*(source.fP));
  66     this->fMapA2 = source.fMapA2;
  67     //this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
  68     //this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
  69     this->fNstrips = source.fNstrips;
  70     this->fPitch   = source.fPitch;
  71     this->fSteps   = source.fSteps;
  72     return *this;
  73 }
  74 //_____________________________________________________________
  75 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
  76   // copy constructor
  77
  78    *this = source;
  79 }
  80 //____________________________________________________________________________
  81 AliITSsimulationSSD::~AliITSsimulationSSD() {
  82   // destructor
  83   //if(fP) delete fP;
  84   //if(fN) delete fN;
  85     delete fMapA2;
  86     //if(fTracksP) delete [] fTracksP;
  87     //if(fTracksN) delete [] fTracksN;
  88     delete fDCS;
  89 }
  90
  91 //_______________________________________________________________
  92 void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
  93                                          Int_t dummy) {
  94   // Digitizes hits for one SSD module
  95
  96   TObjArray *hits = mod->GetHits();
  97   Int_t nhits = hits->GetEntriesFast();
  98   if (!nhits) return;
  99   //cout<<"!! module, nhits ="<<module<<","<<nhits<<endl; //b.b.
 100
 101   Double_t x0=0.0, y0=0.0, z0=0.0;
 102   Double_t x1=0.0, y1=0.0, z1=0.0;
 103   Double_t de=0.0;
 104   Int_t maxNdigits = 2*fNstrips;
 105   Float_t  **pList = new Float_t* [maxNdigits];
 106   memset(pList,0,sizeof(Float_t*)*maxNdigits);
 107   Int_t indexRange[4] = {0,0,0,0};
 108   static Bool_t first=kTRUE;
 109   Int_t lasttrack = -2;
 110   Int_t idtrack = -2;
 111
 112   for(Int_t i=0; i<nhits; i++) {
 113     // LineSegmentL returns 0 if the hit is entering
 114     // If hits is exiting returns positions of entering and exiting hits
 115     // Returns also energy loss
 116     if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
 117       //cout<<"!! befor HitToDigit: hit ="<<i<<endl; //b.b.
 118       HitToDigit(module, x0, y0, z0, x1, y1, z1, de, indexRange, first);
 119
 120       if (lasttrack != idtrack || i==(nhits-1)) {
 121         GetList(idtrack,pList,indexRange);
 122         first=kTRUE;
 123       }
 124       lasttrack=idtrack;
 125     }
 126   }  // end loop over hits
 127
 128   ApplyNoise();
 129   ApplyCoupling();
 130
 131   ChargeToSignal(pList);
 132
 133   fMapA2->ClearMap();
 134 }
 135
 136 //---------------------------------------------------------------
 137 void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
 138                                      Double_t z0, Double_t x1, Double_t y1,
 139                                      Double_t z1, Double_t de,
 140                                      Int_t *indexRange, Bool_t first) {
 141   // Turns hits in SSD module into one or more digits.
 142
 143   AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
 144   // AliITSresponseSSD *res = new AliITSresponseSSD();
 145
 146   Float_t tang[2] = {0.0,0.0};
 147   seg->Angles(tang[0], tang[1]); // stereo<<  -> tan(stereo)~=stereo
 148   //fSegmentation->Angles(tang[0], tang[1]); // stereo<<  -> tan(stereo)~=stereo
 149   Double_t x, y, z;
 150   Double_t dex=0.0, dey=0.0, dez=0.0;
 151   Double_t pairs;
 152   Double_t  ionE = 3.62E-9;   // ionization energy of Si (GeV)
 153   //ionE = (Double_t)  res->GetIonE();
 154
 155   Double_t sigma[2] = {0.,0.};  // standard deviation of the diffusion gaussian
 156
 157   Double_t D[2] = {11.,30.};                // diffusion constant {h,e} (cm**2/sec)
 158   //D[0] = (Double_t) res->GetDiffusionConstantP();
 159   //D[1] = (Double_t) res->GetDiffusionConstantN();
 160
 161   Double_t tdrift[2] = {0.,0.}; // time of drift
 162   Double_t vdrift[2] = {0.86E6,2.28E6};          // drift velocity (cm/sec)
 163   //vdrift[0] = (Double_t) res->GetDriftVelocityP();
 164   //vdrift[1] = (Double_t) res->GetDriftVelocityN();
 165
 166   Double_t w;
 167   Double_t inf[2], sup[2], par0[2];
 168   // Steps in the module are determined "manually" (i.e. No Geant)
 169   // NumOfSteps divide path between entering and exiting hits in steps
 170   Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
 171
 172   // Enery loss is equally distributed among steps
 173   de = de/numOfSteps;
 174   pairs = de/ionE;             // e-h pairs generated
 175
 176   //cout<<"Dy ="<<seg->Dy()<<endl;
 177   //cout<<"numOfSteps ="<<numOfSteps<<endl;
 178   //cout<<"dex,dey,dez ="<<dex<<","<<dey<<","<<dez<<endl;
 179   //cout<<"y0,y1 ="<<y0<<","<<y1<<endl;
 180   for(Int_t j=0; j<numOfSteps; j++) {     // stepping
 181     //    cout<<"step number ="<<j<<endl;
 182     x = x0 + (j+0.5)*dex;
 183     y = y0 + (j+0.5)*dey;
 184     if ( y > (seg->Dy()/2 +10)*1.0E-4 ) {
 185       //if ( y > (seg->Dy()*1.0E-4/2) ) {
 186       //if ( y > (fSegmentation->Dy()*1.0E-4/2) ) {
 187       // check if particle is within the detector
 188       cout<<"AliITSsimulationSSD::HitToDigit: Warning: hit out of detector y0,y,dey,j ="<<y0<<","<<y<<","<<dey<<","<<j<<endl;
 189       return;
 190     };
 191     z = z0 + (j+0.5)*dez;
 192
 193     // calculate drift time
 194     tdrift[0] = (y+(seg->Dy()*1.0E-4)/2) / vdrift[0]; // y is the minimum path
 195     tdrift[1] = ((seg->Dy()*1.0E-4)/2-y) / vdrift[1]; // y is the minimum path
 196     //tdrift[0] = (y+(fSegmentation->Dy()*1.0E-4)/2) / vdrift[0]; // y is the minimum path
 197     //tdrift[1] = ((fSegmentation->Dy()*1.0E-4)/2-y) / vdrift[1]; // y is the minimum path
 198
 199     for(Int_t k=0; k<2; k++) {   // both sides    remember: 0=Pside 1=Nside
 200
 201       tang[k]=TMath::Tan(tang[k]);
 202
 203       // w is the coord. perpendicular to the strips
 204       if(k==0) {
 205         w = (x+(seg->Dx()*1.0E-4)/2) - (z+(seg->Dz()*1.0E-4)/2)*tang[k];
 206         //w = (x+(fSegmentation->Dx()*1.0E-4)/2) - (z+(fSegmentation->Dz()*1.0E-4)/2)*tang[k];
 207         //cout<<"k,x,z,w ="<<k<<","<<x<<","<<z<<","<<w<<endl;
 208       }
 209       else {
 210         w = (x+(seg->Dx()*1.0E-4)/2) + (z-(seg->Dz()*1.0E-4)/2)*tang[k];
 211         //w = (x+(fSegmentation->Dx()*1.0E-4)/2) + (z-(fSegmentation->Dz()*1.0E-4)/2)*tang[k];
 212         //cout<<"k,x,z,w ="<<k<<","<<x<<","<<z<<","<<w<<endl;
 213       }
 214       w = w / (fPitch*1.0E-4); // w is converted in units of pitch
 215
 216       if((w<(-0.5)) || (w>(fNstrips-0.5))) {
 217         // this check rejects hits in regions not covered by strips
 218         // 0.5 takes into account boundaries
 219         if(k==0) cout<<"AliITSsimulationSSD::HitToDigit: Warning: no strip in this region of P side"<<endl;
 220         else cout<<"AliITSsimulationSSD::HitToDigit: Warning: no strip in this region of N side"<<endl;
 221         return;
 222       }
 223
 224       // sigma is the standard deviation of the diffusion gaussian
 225
 226       if(tdrift[k]<0) return;
 227
 228       sigma[k] = TMath::Sqrt(2*D[k]*tdrift[k]);
 229       sigma[k] = sigma[k] /(fPitch*1.0E-4);  //units of Pitch
 230       if(sigma[k]==0.0) {
 231         cout<<"AliITSsimulationSSD::DigitiseModule: Error: sigma=0"<<endl;
 232         exit(0);
 233       }
 234
 235       par0[k] = pairs;
 236       // we integrate the diffusion gaussian from -3sigma to 3sigma
 237       inf[k] = w - 3*sigma[k];      // 3 sigma from the gaussian average
 238       sup[k] = w + 3*sigma[k];      // 3 sigma from the gaussian average
 239       // IntegrateGaussian does the actual integration of diffusion gaussian
 240       IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],
 241                         indexRange, first);
 242     }  // end for loop over side (0=Pside, 1=Nside)
 243   } // end stepping
 244   delete seg;
 245 }
 246
 247 //____________________________________________________________________
 248
 249 void AliITSsimulationSSD::ApplyNoise() {
 250   // Apply Noise.
 251
 252   Float_t signal;
 253   Float_t noise[2] = {0.,0.};
 254   fResponse->GetNoiseParam(noise[0],noise[1]);   // retrieves noise parameters
 255
 256   for(Int_t k=0;k<2;k++){                        // both sides (0=Pside, 1=Nside)
 257     for(Int_t ix=0;ix<fNstrips;ix++){            // loop over strips
 258       signal = (Float_t) fMapA2->GetSignal(k,ix);// retrieves signal from map
 259
 260       signal += gRandom->Gaus(0,noise[k]);       // add noise to signal
 261       if(signal<0.) signal=0.0;                  // in case noise is negative...
 262
 263       fMapA2->SetHit(k,ix,(Double_t)signal);     // give back signal to map
 264     } // loop over strip
 265   } // loop over k (P or N side)
 266 }
 267
 268 //_________________________________________________________________________
 269
 270 void AliITSsimulationSSD::ApplyCoupling() {
 271   // Apply the effect of electronic coupling between channels
 272   Float_t signal, signalLeft=0, signalRight=0;
 273
 274   for(Int_t ix=0;ix<fNstrips;ix++){
 275     if(ix>0.) signalLeft  = (Float_t) fMapA2->GetSignal(0,ix-1)*fDCS->GetCouplingPL();
 276     else signalLeft = 0.0;
 277     if(ix<(fNstrips-1)) signalRight = (Float_t) fMapA2->GetSignal(0,ix+1)*fDCS->GetCouplingPR();
 278     else signalRight = 0.0;
 279     signal = (Float_t) fMapA2->GetSignal(0,ix);
 280     signal += signalLeft + signalRight;
 281     fMapA2->SetHit(0,ix,(Double_t)signal);
 282
 283     if(ix>0.) signalLeft  = (Float_t) fMapA2->GetSignal(1,ix-1)*fDCS->GetCouplingNL();
 284     else signalLeft = 0.0;
 285     if(ix<(fNstrips-1)) signalRight = (Float_t) fMapA2->GetSignal(1,ix+1)*fDCS->GetCouplingNR();
 286     else signalRight = 0.0;
 287     signal = (Float_t) fMapA2->GetSignal(1,ix);
 288     signal += signalLeft + signalRight;
 289     fMapA2->SetHit(1,ix,(Double_t)signal);
 290   } // loop over strips
 291 }
 292
 293 //____________________________________________________________________________
 294 Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
 295   // Computes the integral of a gaussian using Error Function
 296   Float_t sqrt2 = TMath::Sqrt(2.0);
 297   Float_t sigm2 = sqrt2*s;
 298   Float_t integral;
 299
 300   integral = 0.5 * TMath::Erf( (x - av) / sigm2);
 301   return integral;
 302 }
 303
 304 //_________________________________________________________________________
 305 void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w,
 306                                             Double_t sigma,
 307                                             Double_t inf, Double_t sup,
 308                                             Int_t *indexRange, Bool_t first) {
 309   // integrate the diffusion gaussian
 310   // remind: inf and sup are w-3sigma and w+3sigma
 311   //         we could define them here instead of passing them
 312   //         this way we are free to introduce asimmetry
 313
 314   Double_t a=0.0, b=0.0;
 315   Double_t signal = 0.0, dXCharge1 = 0.0, dXCharge2 = 0.0;
 316   // dXCharge1 and 2 are the charge to two neighbouring strips
 317   // Watch that we only involve at least two strips
 318   // Numbers greater than 2 of strips in a cluster depend on
 319   //  geometry of the track and delta rays, not charge diffusion!
 320
 321   Double_t strip = TMath::Floor(w);         // clostest strip on the left
 322
 323   if ( TMath::Abs((strip - w)) < 0.5) {
 324     // gaussian mean is closer to strip on the left
 325     a = inf;                         // integration starting point
 326     if((strip+0.5)<=sup) {
 327       // this means that the tail of the gaussian goes beyond
 328       // the middle point between strips ---> part of the signal
 329       // is given to the strip on the right
 330       b = strip + 0.5;               // integration stopping point
 331       dXCharge1 = F( w, b, sigma) - F(w, a, sigma);
 332       dXCharge2 = F( w, sup, sigma) - F(w ,b, sigma);
 333     }
 334     else {
 335       // this means that all the charge is given to the strip on the left
 336       b = sup;
 337       dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
 338       dXCharge2 = 0.0;
 339     }
 340
 341     dXCharge1 = par * dXCharge1;  // normalize by mean of number of carriers
 342     dXCharge2 = par * dXCharge2;
 343
 344     // for the time being, signal is the charge
 345     // in ChargeToSignal signal is converted in ADC channel
 346     signal = fMapA2->GetSignal(k,strip);
 347     signal += dXCharge1;
 348
 349     fMapA2->SetHit(k,strip,(Double_t)signal);
 350     if(((Int_t) strip) < (fNstrips-1)) {
 351       // strip doesn't have to be the last (remind: last=fNstrips-1)
 352       // otherwise part of the charge is lost
 353       signal = fMapA2->GetSignal(k,(strip+1));
 354       signal += dXCharge2;
 355       fMapA2->SetHit(k,(strip+1),(Double_t)signal);
 356     }
 357
 358     if(dXCharge1 > 1.) {
 359       if (first) {
 360         indexRange[k*2+0]=indexRange[k*2+1]=(Int_t) strip;
 361         first=kFALSE;
 362       }
 363
 364       indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
 365       indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
 366     }      // dXCharge > 1 e-
 367
 368   }
 369   else {
 370     // gaussian mean is closer to strip on the right
 371     strip++;     // move to strip on the rigth
 372     b = sup;     // now you know where to stop integrating
 373     if((strip-0.5)>=inf) {
 374       // tail of diffusion gaussian on the left goes left of
 375       // middle point between strips
 376       a = strip - 0.5;        // integration starting point
 377       dXCharge1 = F(w, b, sigma) - F(w, a, sigma);
 378       dXCharge2 = F(w, a, sigma) - F(w, inf, sigma);
 379     }
 380     else {
 381       a = inf;
 382       dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
 383       dXCharge2 = 0.0;
 384     }
 385
 386     dXCharge1 = par * dXCharge1;    // normalize by means of carriers
 387     dXCharge2 = par * dXCharge2;
 388
 389     // for the time being, signal is the charge
 390     // in ChargeToSignal signal is converted in ADC channel
 391     signal = fMapA2->GetSignal(k,strip);
 392     signal += dXCharge1;
 393     fMapA2->SetHit(k,strip,(Double_t)signal);
 394     if(((Int_t) strip) > 0) {
 395       // strip doesn't have to be the first
 396       // otherwise part of the charge is lost
 397       signal = fMapA2->GetSignal(k,(strip-1));
 398       signal += dXCharge2;
 399       fMapA2->SetHit(k,(strip-1),(Double_t)signal);
 400     }
 401
 402     if(dXCharge1 > 1.) {
 403       if (first) {
 404         indexRange[k*2+0]=indexRange[k*2+1]=(Int_t) strip;
 405         first=kFALSE;
 406       }
 407
 408       indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
 409       indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
 410     }      // dXCharge > 1 e-
 411   }
 412 }
 413
 414
 415   //_________________________________________________________________________
 416
 417 Int_t AliITSsimulationSSD::NumOfSteps(Double_t x, Double_t y, Double_t z,
 418                                  Double_t & dex,Double_t & dey,Double_t & dez) {
 419   // number of steps
 420   // it also returns steps for each coord
 421   //AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
 422
 423   Double_t step = 25E-4;
 424   //step = (Double_t) seg->GetStepSize();  // step size (cm)
 425   Int_t numOfSteps = (Int_t) (TMath::Sqrt(x*x+y*y+z*z)/step);
 426
 427   if (numOfSteps < 1) numOfSteps = 1;       // one step, at least
 428
 429   // we could condition the stepping depending on the incident angle
 430   // of the track
 431   dex = x/numOfSteps;
 432   dey = y/numOfSteps;
 433   dez = z/numOfSteps;
 434
 435   return numOfSteps;
 436
 437 }
 438
 439 //---------------------------------------------
 440 void AliITSsimulationSSD::GetList(Int_t label,Float_t **pList,Int_t *indexRange) {
 441   // loop over nonzero digits
 442   Int_t ix,globalIndex;
 443   Float_t signal=0.;
 444   Float_t highest,middle,lowest;
 445   // printf("SPD-GetList: indexRange[0] indexRange[1] indexRange[2] indexRange[3] %d %d %d %d\n",indexRange[0], indexRange[1], indexRange[2], indexRange[3]);
 446
 447   for(Int_t k=0; k<2; k++) {
 448     for(ix=indexRange[k*2+0];ix<indexRange[k*2+1]+1;ix++){
 449       if(indexRange[k*2+0]<indexRange[k*2+1])
 450         signal=fMapA2->GetSignal(k,ix);
 451
 452       globalIndex = k*fNstrips+ix; // globalIndex starts from 0!
 453       if(!pList[globalIndex]){
 454
 455         //
 456         // Create new list (6 elements - 3 signals and 3 tracks + total sig)
 457         //
 458
 459         pList[globalIndex] = new Float_t [6];
 460         // set list to -1
 461
 462         *pList[globalIndex] = -2.;
 463         *(pList[globalIndex]+1) = -2.;
 464         *(pList[globalIndex]+2) = -2.;
 465         *(pList[globalIndex]+3) =  0.;
 466         *(pList[globalIndex]+4) =  0.;
 467         *(pList[globalIndex]+5) =  0.;
 468
 469         *pList[globalIndex] = (float)label;
 470         *(pList[globalIndex]+3) = signal;
 471
 472       }
 473       else{
 474
 475         // check the signal magnitude
 476
 477         highest = *(pList[globalIndex]+3);
 478         middle = *(pList[globalIndex]+4);
 479         lowest = *(pList[globalIndex]+5);
 480
 481         signal -= (highest+middle+lowest);
 482
 483         //
 484         //  compare the new signal with already existing list
 485         //
 486
 487         if(signal<lowest) continue; // neglect this track
 488
 489         if (signal>highest){
 490           *(pList[globalIndex]+5) = middle;
 491           *(pList[globalIndex]+4) = highest;
 492           *(pList[globalIndex]+3) = signal;
 493
 494           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 495           *(pList[globalIndex]+1) = *pList[globalIndex];
 496           *pList[globalIndex] = label;
 497         }
 498         else if (signal>middle){
 499           *(pList[globalIndex]+5) = middle;
 500           *(pList[globalIndex]+4) = signal;
 501
 502           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
 503           *(pList[globalIndex]+1) = label;
 504         }
 505         else{
 506           *(pList[globalIndex]+5) = signal;
 507           *(pList[globalIndex]+2) = label;
 508         }
 509       }
 510     } // end of loop pixels in x
 511   } // end of loop over pixels in z
 512 }
 513
 514 //---------------------------------------------
 515 void AliITSsimulationSSD::ChargeToSignal(Float_t **pList) {
 516   // charge to signal
 517
 518   AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 519
 520   Float_t threshold = 0.;
 521
 522   Int_t digits[3], tracks[3],hits[3],gi,j1;
 523   Float_t charges[3];
 524   Float_t signal,phys;
 525   Float_t noise[2] = {0.,0.};
 526   fResponse->GetNoiseParam(noise[0],noise[1]);
 527
 528   for(Int_t k=0;k<2;k++){         // both sides (0=Pside, 1=Nside)
 529
 530     // Threshold for zero-suppression
 531     // It can be defined in AliITSresponseSSD
 532     //             threshold = (Float_t)fResponse->MinVal(k);
 533     // I prefer to think adjusting the threshold "manually", looking
 534     // at the scope, and considering noise standard deviation
 535     threshold = 4.0*noise[k];      // 4 times noise is a choice
 536     //cout<<"SSD: k,thresh ="<<k<<","<<threshold<<endl;
 537     for(Int_t ix=0;ix<fNstrips;ix++){         // loop over strips
 538
 539       signal = (Float_t) fMapA2->GetSignal(k,ix);
 540
 541       gi =k*fNstrips+ix; // global index
 542       if (signal > threshold) {
 543         digits[0]=k;
 544         digits[1]=ix;
 545
 546         // convert to ADC signal
 547         // conversion factor are rather arbitrary (need tuning)
 548         // minimum ionizing particle --> ~30000 pairs --> ADC channel 50
 549         signal = signal*50.0/30000.0;
 550         //cout<<"SSD: 1 signal ="<<signal<<endl;
 551         if(signal>1000.) signal = 1000.0; // if exceeding, accumulate last one
 552         //cout<<"SSD: 2 signal ="<<signal<<endl;
 553         digits[2]=(Int_t) signal;
 554
 555         //gi =k*fNstrips+ix; // global index
 556         for(j1=0;j1<3;j1++){
 557           if (pList[gi]) {
 558             tracks[j1] = (Int_t)(*(pList[gi]+j1));
 559           }
 560           else {
 561             tracks[j1]=-2; //noise
 562           }
 563           charges[j1] = 0;
 564         }
 565
 566         phys=0;
 567
 568         hits[0]=0;
 569         hits[1]=0;
 570         hits[2]=0;
 571         aliITS->AddSimDigit(2,phys,digits,tracks,hits,charges);  // finally add digit
 572
 573         //if(pList[gi]) delete [] pList[gi];
 574       }
 575       if(pList[gi]) delete [] pList[gi];
 576     }
 577   }
 578   delete [] pList;
 579 }
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593