ITS/AliITSsimulationSSD.cxx

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