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