ITS/AliITSsimulationSSD.cxx

   1
   2 #include <stdio.h>
   3 #include <TObjArray.h>
   4
   5 #include "AliITSsegmentationSSD.h"
   6 #include "AliITSresponseSSD.h"
   7 #include "AliITSsimulationSSD.h"
   8 #include "AliITSdictSSD.h"
   9 #include "AliITSdcsSSD.h"
  10 #include "AliITS.h"
  11 #include "AliITShit.h"
  12 #include "AliITSdigit.h"
  13 #include "AliITSmodule.h"
  14 #include "AliRun.h"
  15
  16
  17 ClassImp(AliITSsimulationSSD);
  18 //------------------------------------------------------------
  19 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
  20                                          AliITSresponse *resp){
  21   // Constructor
  22
  23
  24     fSegmentation = seg;
  25     fResponse = resp;
  26     fDCS = new AliITSdcsSSD(seg,resp);
  27
  28     fNstrips = fSegmentation->Npx();
  29     fPitch = fSegmentation->Dpx(0);
  30
  31     fP = new TArrayF(fNstrips+1);
  32     fN = new TArrayF(fNstrips+1);
  33
  34     fTracksP = new AliITSdictSSD[fNstrips+1];
  35     fTracksN = new AliITSdictSSD[fNstrips+1];
  36
  37
  38     fSteps  = 10;   // still hard-wired - set in SetDetParam and get it via
  39                      // fDCS together with the others eventually
  40
  41 }
  42 //___________________________________________________________________________
  43 AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD
  44                                                                       &source){
  45 // Operator =
  46     if(this==&source) return *this;
  47
  48     this->fDCS = new AliITSdcsSSD(*(source.fDCS));
  49     this->fN   = new TArrayF(*(source.fN));
  50     this->fP   = new TArrayF(*(source.fP));
  51     this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
  52     this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
  53     this->fNstrips = source.fNstrips;
  54     this->fPitch   = source.fPitch;
  55     this->fSteps   = source.fSteps;
  56     return *this;
  57 }
  58 //_____________________________________________________________
  59 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
  60   // copy constructor
  61    *this = source;
  62 }
  63 //____________________________________________________________________________
  64 AliITSsimulationSSD::~AliITSsimulationSSD() {
  65   // anihilator
  66
  67
  68     if(fP) delete fP;
  69     if(fN) delete fN;
  70
  71     if(fTracksP) delete [] fTracksP;
  72     if(fTracksN) delete [] fTracksN;
  73
  74     delete fDCS;
  75
  76 }
  77 //_______________________________________________________________
  78 //
  79 // Hit to digit
  80 //_______________________________________________________________
  81 //
  82 void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
  83                                          Int_t dummy) {
  84   // Digitizes one SSD module of hits.
  85
  86     TObjArray *hits = mod->GetHits();
  87     Int_t nhits = hits->GetEntriesFast();
  88     if (!nhits) return;
  89
  90     //printf("simSSD: module nhits %d %d\n",module,nhits);
  91
  92     Int_t i;
  93     for(i=0; i<fNstrips; i++) {
  94        (*fP)[i] = 0;
  95        (*fN)[i] = 0;
  96        fTracksP[i].ZeroTracks();
  97        fTracksN[i].ZeroTracks();
  98     }
  99
 100     for(i=0; i<nhits; i++) {
 101        Int_t idtrack=mod->GetHitTrackIndex(i);
 102        HitToDigit(i,idtrack,nhits,hits);
 103     }
 104
 105
 106
 107     ApplyNoise();
 108     ApplyCoupling();
 109     ApplyThreshold();
 110     ApplyDAQ();
 111
 112
 113 }
 114
 115 //---------------------------------------------------------------
 116
 117 void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
 118                                      Int_t nhits,TObjArray *hits) {
 119   // Turns one or more hits in an SSD module into one or more digits.
 120
 121     Int_t      stripP, stripN, i;
 122     Float_t    dsP, dsN;
 123     Float_t    sP, sN;
 124     Float_t    eP, eN;
 125     Float_t    arrayEP[10];         // hard-wired number of steps
 126     Float_t    arrayEN[10];
 127     Int_t      track = -3;
 128
 129     Float_t    ionization = 0;
 130     Float_t    signal;
 131
 132     AliITSdictSSD *dict;
 133
 134
 135     // check if this is the right order !!!!!
 136
 137     AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
 138     AliITShit *hitE = (AliITShit*)hits->At(hitNo);
 139
 140
 141     while (!((hitE->StatusExiting()) ||
 142              (hitE->StatusDisappeared()) ||
 143              (hitE->StatusStop()))) {
 144
 145         if (++hitNo<nhits) {
 146            ionization = hitE->GetIonization();
 147            hitE = (AliITShit*)hits->At(hitNo);
 148         }
 149     }
 150
 151
 152     if (hitI->GetTrack() == hitE->GetTrack())
 153        //track = idtrack;
 154        track = hitI->GetTrack();
 155     else
 156        printf("!!! Emergency !!!\n");
 157
 158
 159     ionization += hitE->GetIonization();
 160
 161     const Float_t kconvm=10000.;  // cm -> microns
 162
 163     Float_t xI, yI, zI;
 164     hitI->GetPositionL(xI, yI, zI);
 165
 166     xI *= kconvm;
 167     yI *= kconvm;
 168     zI *= kconvm;
 169
 170     Float_t xE, yE, zE;
 171     hitE->GetPositionL(xE, yE, zE);
 172
 173     xE *= kconvm;
 174     yE *= kconvm;
 175     zE *= kconvm;
 176
 177     Float_t dx = (xE - xI);
 178     Float_t dz = (zE - zI);
 179
 180
 181     // Debuging
 182     /*
 183     fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
 184
 185        printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d  %d\n",
 186              hitNo, xI, zI, dx, dz,
 187              stripP, stripN, track);
 188      printf("%10.5f %10d \n", ionization, hitI->fTrack);
 189     */
 190
 191     // end of debuging
 192
 193
 194     eP=0;
 195     eN=0;
 196
 197     for (i=0; i<fSteps; i++) {
 198
 199       //        arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
 200       //        arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
 201         arrayEP[i] = ionization/fSteps;
 202         arrayEN[i] = ionization/fSteps;
 203
 204         eP += arrayEP[i];
 205         eN += arrayEN[i];
 206     }
 207
 208     const Float_t kconv = 1.0e9 / 3.6;  // GeV -> e-hole pairs
 209
 210     for(i=0; i<fSteps; i++) {
 211
 212         arrayEP[i] = kconv * arrayEP[i] * (ionization / eP);
 213         arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);
 214     }
 215
 216     dx /= fSteps;
 217     dz /= fSteps;
 218
 219     Float_t sigmaP, sigmaN;
 220     fResponse->SigmaSpread(sigmaP,sigmaN);
 221
 222     //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);
 223
 224     Float_t noiseP, noiseN;
 225     fResponse->GetNoiseParam(noiseP,noiseN);
 226
 227     //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);
 228
 229     for(i=0; i<fSteps; i++) {
 230
 231         Int_t j;
 232
 233         fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
 234         //printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
 235         dsP    = Get2Strip(1,stripP,xI, zI); // Between 0-1
 236         dsN    = Get2Strip(0,stripN,xI, zI); // Between 0-1
 237
 238         //sP = sigmaP * sqrt(300. * i / (fSteps));
 239         //sN = sigmaN * sqrt(300. * i /(fSteps-i));
 240
 241         sP = sigmaP * sqrt(300. * (i+1) / (fSteps));
 242         sN = sigmaN * sqrt(300. * (i+1) /(fSteps-i));
 243
 244
 245         sP = (i<2        && dsP>0.3 && dsP<0.7)? 20. : sP;  // square of (microns)
 246         sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN;  // square of (microns)
 247
 248         sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP;  // square of (microns)
 249         sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN;  // square of (microns)
 250
 251
 252         //printf("i=%d SigmaP SigmaN sP sN %f %f %e %e\n",i,sigmaP, sigmaN,sP,sN);
 253
 254         for (j=-1; j<2; j++) {
 255             if (stripP+j<0 || stripP+j>fNstrips) continue;
 256             signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
 257             //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
 258             if (signal > noiseP/fSteps) {
 259                (*fP)[stripP+j] += signal;
 260                dict = (fTracksP+stripP+j);
 261                (*dict).AddTrack(track);
 262             }
 263         }  // end for j loop over neighboring strips
 264
 265         for (j=-1; j<2; j++) {
 266             if (stripN+j<0 || stripN+j>fNstrips) continue;
 267             signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
 268             //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
 269             if (signal > noiseN/fSteps) {
 270                (*fN)[stripN+j] += signal;
 271                dict = (fTracksN+stripN+j);    //co to jest
 272                (*dict).AddTrack(track);
 273             }
 274         }  // end for j loop over neighboring strips
 275
 276         xI += dx;
 277         zI += dz;
 278     }
 279
 280
 281 }
 282
 283
 284 //____________________________________________________________________
 285 //
 286 //  Private Methods for Simulation
 287 //______________________________________________________________________
 288 //
 289
 290 void AliITSsimulationSSD::ApplyNoise() {
 291   // Apply Noise.
 292    Float_t noiseP, noiseN;
 293    fResponse->GetNoiseParam(noiseP,noiseN);
 294
 295     Int_t i;
 296     for(i = 0; i<fNstrips; i++) {
 297        (*fP)[i] += gRandom->Gaus(0,noiseP);
 298        (*fN)[i] += gRandom->Gaus(0,noiseN);
 299     }
 300 }
 301
 302 //_________________________________________________________________________
 303
 304 void AliITSsimulationSSD::ApplyCoupling() {
 305   // Apply the effecto of electronic coupling between channels
 306     Int_t i;
 307     for(i = 1; i<fNstrips-1; i++) {
 308       (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
 309       (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
 310     }
 311 }
 312
 313 //__________________________________________________________________________
 314
 315 void AliITSsimulationSSD::ApplyThreshold() {
 316   // Applies the effect of a threshold on the signals for digitization.
 317    Float_t noiseP, noiseN;
 318    fResponse->GetNoiseParam(noiseP,noiseN);
 319
 320    // or introduce the SetThresholds in fResponse
 321
 322     Int_t i;
 323     for(i=0; i<fNstrips; i++) {
 324        (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
 325        (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0;
 326     }
 327
 328 }
 329
 330 //__________________________________________________________________________
 331
 332 void AliITSsimulationSSD::ApplyDAQ() {
 333   // Converts simulated signals to simulated ADC counts
 334     AliITS *its=(AliITS*)gAlice->GetModule("ITS");
 335
 336     Float_t noiseP, noiseN;
 337     fResponse->GetNoiseParam(noiseP,noiseN);
 338
 339     char opt[30],dummy[20];
 340     fResponse->ParamOptions(opt,dummy);
 341
 342     Int_t i,j;
 343     if (strstr(opt,"SetInvalid")) {
 344       printf("invalid option %s\n",opt);
 345       // Set signal = 0 if invalid strip
 346       for(i=0; i<fNstrips; i++) {
 347          if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
 348          if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
 349       }
 350     }
 351
 352     Int_t digits[3], tracks[3], hits[3];
 353     Float_t charges[3];
 354     Float_t phys=0;
 355     for(i=0;i<3;i++) tracks[i]=-3;
 356     for(i=0; i<fNstrips; i++) {
 357        if( (strstr(opt,"SetInvalid") && (*fP)[i] < noiseP*4) || !(*fP)[i]) continue;
 358           digits[0]=1;
 359           digits[1]=i;
 360           digits[2]=(int)(*fP)[i];
 361           for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
 362             if(j>2) continue;
 363             if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
 364             else tracks[j]=-2;
 365             //printf("P side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
 366             charges[j] = 0;
 367             hits[j] = -1;
 368           }
 369           its->AddSimDigit(2,phys,digits,tracks,hits,charges);
 370
 371           //cout << (fTracksP+i)->GetNTracks();
 372           //
 373           //if ((fTracksP+i)->GetNTracks() == 0) {
 374           //   cout << d.fCoord2 << " " << d.fSignal << "\n";
 375           //}
 376     }
 377
 378
 379     for(i=0; i<fNstrips; i++) {
 380        if( (strstr(opt,"SetInvalid") && (*fN)[i] < noiseN*4)|| !(*fN)[i]) continue;
 381           digits[0]=0;
 382           digits[1]=i;
 383           digits[2]=(int)(*fN)[i];
 384           for( j=0; j<(fTracksN+i)->GetNTracks(); j++) {
 385             if(j>2) continue;
 386             if((fTracksN+i)->GetNTracks()) tracks[j]=(fTracksN+i)->GetTrack(j);
 387             else tracks[j]=-2;
 388             //printf("N side: i,j,tracks[j] %d %d %d\n",i,j,tracks[j]);
 389             charges[j] = 0;
 390             hits[j] = -1;
 391           }
 392           its->AddSimDigit(2,phys,digits,tracks,hits,charges);
 393
 394           //cout << (fTracksN+i)->GetNTracks();
 395           //if ((fTracksN+i)->GetNTracks() == 0) {
 396           //   cout << d.fCoord2 << " " << d.fSignal << "\n";
 397           //}
 398     }
 399
 400 }
 401
 402
 403 //____________________________________________________________________________
 404
 405 Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
 406   // Computes the integral of a gaussian at the mean valuse x with sigma s.
 407     //printf("SDD:F(%e,%e)\n",x,s);
 408
 409     Float_t fval=0;
 410     if(s) fval=0.5*TMath::Erf(x * fPitch / s) ;
 411     else {
 412       Error("SSD simulation: F","sigma is zero!!!",s);
 413     }
 414     return fval;
 415 }
 416
 417 //______________________________________________________________________
 418
 419 Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
 420   // Returns the relative space between two strips.
 421
 422     // flag==1 for Pside, 0 for Nside
 423
 424     Float_t stereoP, stereoN;
 425     fSegmentation->Angles(stereoP,stereoN);
 426
 427     Float_t tanP=TMath::Tan(stereoP);
 428     Float_t tanN=TMath::Tan(stereoN);
 429
 430     Float_t dx = fSegmentation->Dx();
 431     Float_t dz = fSegmentation->Dz();
 432
 433
 434      x += dx/2;
 435      z += dz/2;
 436
 437      if (flag) return (x - z*tanP) / fPitch - iStrip;       // from 0 to 1
 438      else  return (x - tanN*(dz - z)) / fPitch - iStrip;
 439 }
 440 //____________________________________________________________________________