Update for Ds
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
CommitLineData
57817f7c 1/**************************************************************************
2727462b 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 **************************************************************************/
88cb7938 15
57817f7c 16/* $Id$ */
17
b0f5e3fc 18#include <stdio.h>
8a00af9a 19#include <stdlib.h>
4ae5bbc4 20#include <Riostream.h>
b0f5e3fc 21#include <TObjArray.h>
e939a978 22#include <TRandom.h>
23
66b89079 24#include <TGeoGlobalMagField.h>
fd61217e 25#include "AliITSmodule.h"
57817f7c 26#include "AliITSMapA2.h"
c7a4dac0 27#include "AliITSpList.h"
fcf95fc7 28#include "AliITSCalibrationSSD.h"
c7a4dac0 29#include "AliITSsegmentationSSD.h"
14a74335 30//#include "AliITSdcsSSD.h"
b0f5e3fc 31#include "AliITS.h"
3a97c582 32#include "AliITShit.h"
e869281d 33#include "AliITSdigitSSD.h"
b0f5e3fc 34#include "AliRun.h"
66b89079 35#include "AliMagF.h"
569a17d8 36#include "AliITSgeom.h"
57817f7c 37#include "AliITSsimulationSSD.h"
84541af4 38#include "AliITSTableSSD.h"
7d3b8602 39#include <TF1.h>
bb7e41dd 40#include "AliMathBase.h"
b0f5e3fc 41
925e6570 42ClassImp(AliITSsimulationSSD)
fd61217e 43////////////////////////////////////////////////////////////////////////
8ba39da9 44// //
14a74335 45// Author: Enrico Fragiacomo //
46// enrico.fragiacomo@ts.infn.it //
82adfb7d 47// Last revised: june 2008 //
14a74335 48// //
49// AliITSsimulationSSD is the simulation of SSD. //
8ba39da9 50////////////////////////////////////////////////////////////////////////
fd61217e 51
0315d466 52//----------------------------------------------------------------------
aacedc3e 53AliITSsimulationSSD::AliITSsimulationSSD():AliITSsimulation(),
14a74335 54 //fDCS(0),
aacedc3e 55fMapA2(0),
56fIonE(0.0),
57fDifConst(),
e06798d9 58fDriftVel(),
66b89079 59fTimeResponse(NULL),
60fLorentz(kFALSE),
61fTanLorAngP(0),
62fTanLorAngN(0)
63{
2727462b 64 //default Constructor
65 //Inputs:
66 // none.
67 // Outputs:
68 // none.
69 // Return:
70 // A default construction AliITSsimulationSSD class
57817f7c 71}
72//----------------------------------------------------------------------
8ba39da9 73AliITSsimulationSSD::AliITSsimulationSSD(AliITSDetTypeSim* dettyp):
74AliITSsimulation(dettyp),
14a74335 75//fDCS(0),
aacedc3e 76fMapA2(0),
77fIonE(0.0),
78fDifConst(),
7d3b8602 79fDriftVel(),
66b89079 80fTimeResponse(NULL),
81fLorentz(kFALSE),
82fTanLorAngP(0),
83fTanLorAngN(0)
84{
2727462b 85 // Constructor
86 // Input:
87 // AliITSDetTypeSim Pointer to the SSD dettype to be used
88 // Outputs:
89 // none.
90 // Return
91 // A standard constructed AliITSsimulationSSD class
92
7d3b8602 93 fTimeResponse = new TF1("ftimeresponse",".5*x*exp(1.-.5*x)");
2727462b 94 Init();
c7a4dac0 95}
96//----------------------------------------------------------------------
aacedc3e 97void AliITSsimulationSSD::Init(){
2727462b 98 // Inilizer, Inilizes all of the variable as needed in a standard place.
99 // Input:
100 // AliITSsegmentationSSD *seg Pointer to the SSD segmentation to be used
101 // AliITSCalibrationSSD *resp Pointer to the SSD responce class to be used
102 // Outputs:
103 // none.
104 // Return
105 // none.
106 AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
66b89079 107 AliITSSimuParam* simpar = fDetType->GetSimuParam();
2727462b 108
109 SetDriftVelocity(); // use default values in .h file
110 SetIonizeE(); // use default values in .h file
111 SetDiffConst(); // use default values in .h file
112 fpList = new AliITSpList(2,GetNStrips());
113 fMapA2 = new AliITSMapA2(seg);
66b89079 114 SetLorentzDrift(simpar->GetSSDLorentzDrift());
115 if (fLorentz) SetTanLorAngle();
b0f5e3fc 116}
66b89079 117
118//______________________________________________________________________
119Bool_t AliITSsimulationSSD::SetTanLorAngle() {
120 // This function set the Tangent of the Lorentz angles.
121 // output: Bool_t : kTRUE in case of success
122 //
123
124 if(!fDetType) {
125 AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set ");
126 return kFALSE;}
127
128 AliITSSimuParam* simpar = fDetType->GetSimuParam();
129 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
130 if (!fld) AliFatal("The field is not initialized");
131 Double_t bz = fld->SolenoidField();
132
133 fTanLorAngN = TMath::Tan( simpar->LorentzAngleElectron(bz) );
134 fTanLorAngP = TMath::Tan( simpar->LorentzAngleHole(bz) );
135
136 return kTRUE;
137}
138
0315d466 139//______________________________________________________________________
57817f7c 140AliITSsimulationSSD& AliITSsimulationSSD::operator=(
2727462b 141 const AliITSsimulationSSD &s){
142 // Operator =
143
144 if(this==&s) return *this;
145
146 // this->fDCS = new AliITSdcsSSD(*(s.fDCS));
147 this->fMapA2 = s.fMapA2;
148 this->fIonE = s.fIonE;
149 this->fDifConst[0] = s.fDifConst[0];
150 this->fDifConst[1] = s.fDifConst[1];
151 this->fDriftVel[0] = s.fDriftVel[0];
152 this->fDriftVel[1] = s.fDriftVel[1];
7d3b8602 153 this->fTimeResponse = s.fTimeResponse;
66b89079 154 this->fLorentz = s.fLorentz;
155 this->fTanLorAngP = s.fTanLorAngP;
156 this->fTanLorAngN = s.fTanLorAngN;
2727462b 157 return *this;
b0f5e3fc 158}
85f5e9c2 159/*
57817f7c 160//______________________________________________________________________
5402d9ca 161AliITSsimulation& AliITSsimulationSSD::operator=(
2727462b 162 const AliITSsimulation &s){
163 // Operator =
d2f55a22 164
2727462b 165 if(this==&s) return *this;
166 Error("AliITSsimulationSSD","Not allowed to make a = with "
d2f55a22 167 "AliITSsimulationSSD Using default creater instead");
2727462b 168
169 return *this;
d2f55a22 170}
85f5e9c2 171*/
d2f55a22 172//______________________________________________________________________
ac74f489 173AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source):
2727462b 174 AliITSsimulation(source),
7537d03c 175fMapA2(source.fMapA2),
176fIonE(source.fIonE),
177fDifConst(),
7d3b8602 178fDriftVel(),
66b89079 179fTimeResponse(source.fTimeResponse),
180fLorentz(source.fLorentz),
181fTanLorAngP(source.fTanLorAngP),
182fTanLorAngN(source.fTanLorAngN)
183{
2727462b 184 // copy constructor
185 fDifConst[0] = source.fDifConst[0];
186 fDifConst[1] = source.fDifConst[1];
187 fDriftVel[0] = source.fDriftVel[0];
188 fDriftVel[1] = source.fDriftVel[1];
b0f5e3fc 189}
0315d466 190//______________________________________________________________________
b0f5e3fc 191AliITSsimulationSSD::~AliITSsimulationSSD() {
2727462b 192 // destructor
193 delete fMapA2;
7d3b8602 194 delete fTimeResponse;
2727462b 195 //delete fDCS;
0315d466 196}
57817f7c 197//______________________________________________________________________
3a97c582 198void AliITSsimulationSSD::InitSimulationModule(Int_t module,Int_t event){
2727462b 199 // Creates maps to build the list of tracks for each sumable digit
200 // Inputs:
201 // Int_t module // Module number to be simulated
202 // Int_t event // Event number to be simulated
203 // Outputs:
204 // none.
205 // Return
206 // none.
207
208 SetModuleNumber(module);
209 SetEventNumber(event);
210 fMapA2->ClearMap();
211 fpList->ClearMap();
3a97c582 212}
213//______________________________________________________________________
214void AliITSsimulationSSD::FinishSDigitiseModule(){
2727462b 215 // Does the Sdigits to Digits work
216 // Inputs:
217 // none.
218 // Outputs:
219 // none.
220 // Return:
221 // none.
222
223 FillMapFrompList(fpList); // need to check if needed here or not????
224 SDigitToDigit(fModule,fpList);
225 fpList->ClearMap();
226 fMapA2->ClearMap();
3a97c582 227}
228//______________________________________________________________________
aacedc3e 229void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t,Int_t) {
2727462b 230 // Digitizes hits for one SSD module
231 SetModuleNumber(mod->GetIndex());
232
233 HitsToAnalogDigits(mod,fpList);
234 SDigitToDigit(GetModuleNumber(),fpList);
235
236 fpList->ClearMap();
237 fMapA2->ClearMap();
c7a4dac0 238}
239//______________________________________________________________________
aacedc3e 240void AliITSsimulationSSD::SDigitiseModule(AliITSmodule *mod,Int_t,Int_t) {
2727462b 241 // Produces Summable/Analog digits and writes them to the SDigit tree.
c7a4dac0 242
2727462b 243 HitsToAnalogDigits(mod,fpList);
c7a4dac0 244
2727462b 245 WriteSDigits(fpList);
246
247 fpList->ClearMap();
248 fMapA2->ClearMap();
c7a4dac0 249}
250//______________________________________________________________________
251void AliITSsimulationSSD::SDigitToDigit(Int_t module,AliITSpList *pList){
2727462b 252 // Takes the pList and finishes the digitization.
253
254 ApplyNoise(pList,module);
255 ApplyCoupling(pList,module);
256 ApplyDeadChannels(module);
257
258 ChargeToSignal(module,pList);
c7a4dac0 259}
260//______________________________________________________________________
261void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
2727462b 262 AliITSpList *pList){
263 // Loops over all hits to produce Analog/floating point digits. This
264 // is also the first task in producing standard digits.
265 Int_t lasttrack = -2;
266 Int_t idtrack = -2;
267 Double_t x0=0.0, y0=0.0, z0=0.0;
268 Double_t x1=0.0, y1=0.0, z1=0.0;
269 Double_t de=0.0;
270 Int_t module = mod->GetIndex();
7d3b8602 271 Double_t tof = 0.;
272
2727462b 273
274 AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
275
276 TObjArray *hits = mod->GetHits();
277 Int_t nhits = hits->GetEntriesFast();
278 if (nhits<=0) return;
279 AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips());
280 module = mod->GetIndex();
281 if ( mod->GetLayer() == 6 ) seg->SetLayer(6);
282 if ( mod->GetLayer() == 5 ) seg->SetLayer(5);
66b89079 283
2727462b 284 for(Int_t i=0; i<nhits; i++) {
285 // LineSegmentL returns 0 if the hit is entering
286 // If hits is exiting returns positions of entering and exiting hits
287 // Returns also energy loss
288 if(GetDebug(4)){
289 cout << i << " ";
290 cout << mod->GetHit(i)->GetXL() << " "<<mod->GetHit(i)->GetYL();
291 cout << " " << mod->GetHit(i)->GetZL();
292 cout << endl;
293 } // end if
294 if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
7d3b8602 295
296 // Scale down dE/dx according to the hit's TOF wrt to the trigger
297 // Necessary for pileup simulation
298 // EF - 21/04/09
299 tof = mod->GetHit(i)->GetTOF();
300 tof *= 1.E+6; // convert time in microsecond
c9f26a45 301 if(tof<2.) de = de * fTimeResponse->Eval(-1.*tof+2.);
302 else de = 0.;
7d3b8602 303 //
304
2727462b 305 HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav);
306 if (lasttrack != idtrack || i==(nhits-1)) {
82adfb7d 307 GetList(idtrack,i,module,pList,tav);
2727462b 308 } // end if
309 lasttrack=idtrack;
310 } // end if
311 } // end loop over hits
312 delete tav; tav=0;
313 return;
b0f5e3fc 314}
0315d466 315//----------------------------------------------------------------------
fd61217e 316void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
2727462b 317 Double_t z0, Double_t x1, Double_t y1,
318 Double_t z1, Double_t de,
319 AliITSTableSSD *tav) {
320
321 // hit to digit conversion
322
323 AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
324 // Turns hits in SSD module into one or more digits.
325 //Float_t tang[2] = {0.0,0.0};
326 //seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo
327 Double_t x, y, z;
328 Double_t dex=0.0, dey=0.0, dez=0.0;
329 Double_t pairs; // pair generation energy per step.
330 Double_t sigma[2] = {0.,0.};// standard deviation of the diffusion gaussian
331 Double_t tdrift[2] = {0.,0.}; // time of drift
332 Double_t w;
333 Double_t inf[2], sup[2], par0[2];
66b89079 334
335 // Set up corrections for Lorentz drift (ExB)
336 Double_t TanLorAngP = fTanLorAngP;
337 Double_t TanLorAngN = fTanLorAngN;
338 if(seg->GetLayer()==6) {
339 TanLorAngP = -1.*fTanLorAngP;
340 TanLorAngN = -1.*fTanLorAngN;
341 }
342
2727462b 343 // Steps in the module are determined "manually" (i.e. No Geant)
344 // NumOfSteps divide path between entering and exiting hits in steps
345 Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
346 // Enery loss is equally distributed among steps
347 de = de/numOfSteps;
348 pairs = de/GetIonizeE(); // e-h pairs generated
349
350 //-----------------------------------------------------
351 // stepping
352 //-----------------------------------------------------
353 for(Int_t j=0; j<numOfSteps; j++) { // stepping
354
355 x = x0 + (j+0.5)*dex;
356 y = y0 + (j+0.5)*dey;
357 if ( y > (seg->Dy()/2+10)*1.0E-4 ) {
358 // check if particle is within the detector
359 Warning("HitToDigit",
d35ad08f 360 "hit out of detector y0=%e,y=%e,dey=%e,j =%d module=%d, exceed=%e",
361 y0,y,dey,j,module, y-(seg->Dy()/2+10)*1.0E-4);
2727462b 362 return;
363 } // end if
364 z = z0 + (j+0.5)*dez;
82adfb7d 365
2727462b 366 if(GetDebug(4)) cout <<"HitToDigit "<<x<<" "<<y<<" "<<z<< " "
14a74335 367 <<dex<<" "<<dey<<" "<<dez<<endl;
82adfb7d 368
2727462b 369 if(seg->GetLayer()==6) {
370 y=-y; // Lay6 module has sensor up-side-down!!!
371 }
372
2727462b 373 Int_t k;
374 //---------------------------------------------------------
375 // Pside
376 //------------------------------------------------------------
377 k=0;
66b89079 378
379 // w is the coord. perpendicular to the strips
380 // Float_t xp=x*1.e+4,zp=z*1.e+4; // microns
381 Float_t xp=x,zp=z;
382
383 // correction for the Lorentz's angle
384 if(fLorentz) {
385 Float_t deltaxp = (y+(seg->Dy()*1.0E-4)/2)*TanLorAngP;
386 xp+=deltaxp;
387 }
388
389 seg->GetPadTxz(xp,zp);
2727462b 390
391 // calculate drift time
392 // y is the minimum path
393 tdrift[0] = (y+(seg->Dy()*1.0E-4)/2)/GetDriftVelocity(0);
394
395 w = xp; // P side strip number
396
397 if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
398 // this check rejects hits in regions not covered by strips
399 // 0.5 takes into account boundaries
400 if(GetDebug(4)) cout << "Dead SSD region, x,z="<<x<<","<<z<<endl;
401 return; // There are dead region on the SSD sensitive volume!!!
402 } // end if
66b89079 403 // sigma is the standard deviation of the diffusion gaussian
2727462b 404 if(tdrift[k]<0) return;
405
406 sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
407 sigma[k] /= (GetStripPitch()*1.0E-4); //units of Pitch
408
409 if(sigma[k]==0.0) {
410 Error("HitToDigit"," sigma[%d]=0",k);
411 exit(0);
412 } // end if
413
414 par0[k] = pairs;
415 // we integrate the diffusion gaussian from -3sigma to 3sigma
416 inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average
417 sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
418 // IntegrateGaussian does the actual
419 // integration of diffusion gaussian
420 IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
421
422 //------------------------------------------------------
423 // end Pside
424 //-------------------------------------------------------
425
426 //------------------------------------------------------
427 // Nside
428 //-------------------------------------------------------
429 k=1;
66b89079 430
431 xp=x; zp=z;
432
433 // correction for the Lorentz's angle
434 if(fLorentz) {
435 Float_t deltaxn = ((seg->Dy()*1.0E-4)/2-y)*TanLorAngN;
436 xp+=deltaxn;
437 }
438
439
440 seg->GetPadTxz(xp,zp);
441
2727462b 442 tdrift[1] = ((seg->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1);
443
444 //tang[k]=TMath::Tan(tang[k]);
445
446 w = zp; // N side strip number
447
448 if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
449 // this check rejects hits in regions not covered by strips
450 // 0.5 takes into account boundaries
451 if(GetDebug(4)) cout << "Dead SSD region, x,z="<<x<<","<<z<<endl;
452 return; // There are dead region on the SSD sensitive volume.
453 } // end if
454
455 // sigma is the standard deviation of the diffusion gaussian
456 if(tdrift[k]<0) return;
457
458 sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
459 sigma[k] /= (GetStripPitch()*1.0E-4); //units of Pitch
460
461 if(sigma[k]==0.0) {
462 Error("HitToDigit"," sigma[%d]=0",k);
463 exit(0);
464 } // end if
465
466 par0[k] = pairs;
467 // we integrate the diffusion gaussian from -3sigma to 3sigma
468 inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average
469 sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
470 // IntegrateGaussian does the actual
471 // integration of diffusion gaussian
472 IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
473
474 //-------------------------------------------------
475 // end Nside
476 //-------------------------------------------------
477
478
479 } // end stepping
b0f5e3fc 480}
14a74335 481
57817f7c 482//______________________________________________________________________
c7a4dac0 483void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){
2727462b 484 // Apply Noise.
485 Int_t ix;
486 Double_t signal,noise;
487 AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
488
489 // Pside
490 for(ix=0;ix<GetNStrips();ix++){ // loop over strips
491
492 // noise is gaussian
493 noise = (Double_t) gRandom->Gaus(0,res->GetNoiseP(ix));
494
495 // need to calibrate noise
496 // NOTE. noise from the calibration database comes uncalibrated,
497 // it needs to be calibrated in order to be added
498 // to the signal. It will be decalibrated later on together with the noise
499 noise *= (Double_t) res->GetGainP(ix);
500
501 // noise comes in ADC channels from the calibration database
502 // It needs to be converted back to electronVolts
503 noise /= res->GetSSDDEvToADC(1.);
504
505 // Finally, noise is added to the signal
506 signal = noise + fMapA2->GetSignal(0,ix);//get signal from map
507 fMapA2->SetHit(0,ix,signal); // give back signal to map
508 if(signal>0.0) pList->AddNoise(0,ix,module,noise);
509 } // loop over strip
510
511 // Nside
512 for(ix=0;ix<GetNStrips();ix++){ // loop over strips
513 noise = (Double_t) gRandom->Gaus(0,res->GetNoiseN(ix));// give noise to signal
514 noise *= (Double_t) res->GetGainN(ix);
515 noise /= res->GetSSDDEvToADC(1.);
516 signal = noise + fMapA2->GetSignal(1,ix);//get signal from map
517 fMapA2->SetHit(1,ix,signal); // give back signal to map
518 if(signal>0.0) pList->AddNoise(1,ix,module,noise);
519 } // loop over strip
520
b0f5e3fc 521}
0315d466 522//______________________________________________________________________
c7a4dac0 523void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) {
2727462b 524 // Apply the effect of electronic coupling between channels
525 Int_t ix;
526 Double_t signal=0;
527 //AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
528 AliITSSimuParam* res = fDetType->GetSimuParam();
529
530 Double_t *contrLeft = new Double_t[GetNStrips()];
531 Double_t *contrRight = new Double_t[GetNStrips()];
532
533 // P side coupling
534 for(ix=0;ix<GetNStrips();ix++){
535 if(ix>0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetSSDCouplingPL();
536 else contrLeft[ix] = 0.0;
537 if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetSSDCouplingPR();
538 else contrRight[ix] = 0.0;
539 } // loop over strips
540
541 for(ix=0;ix<GetNStrips();ix++){
542 signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingPL() * fMapA2->GetSignal(0,ix)
543 - res->GetSSDCouplingPR() * fMapA2->GetSignal(0,ix);
544 fMapA2->AddSignal(0,ix,signal);
545 if(signal>0.0) pList->AddNoise(0,ix,module,signal);
546 } // loop over strips
547
548 // N side coupling
549 for(ix=0;ix<GetNStrips();ix++){
550 if(ix>0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetSSDCouplingNL();
551 else contrLeft[ix] = 0.0;
552 if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetSSDCouplingNR();
553 else contrRight[ix] = 0.0;
554 } // loop over strips
555
556 for(ix=0;ix<GetNStrips();ix++){
557 signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingNL() * fMapA2->GetSignal(0,ix)
558 - res->GetSSDCouplingNR() * fMapA2->GetSignal(0,ix);
559 fMapA2->AddSignal(1,ix,signal);
560 if(signal>0.0) pList->AddNoise(1,ix,module,signal);
561 } // loop over strips
562
563
564 delete [] contrLeft;
565 delete [] contrRight;
b0f5e3fc 566}
14a74335 567
568//______________________________________________________________________
569void AliITSsimulationSSD::ApplyDeadChannels(Int_t module) {
2727462b 570 // Kill dead channels setting gain to zero
14a74335 571
2727462b 572 AliITSCalibrationSSD* res = (AliITSCalibrationSSD*)GetCalibrationModel(module);
14a74335 573
2727462b 574 for(Int_t i=0;i<GetNStrips();i++){
14a74335 575
2727462b 576 if(res->IsPChannelBad(i)) res->SetGainP(i,0.0);
577 if(res->IsNChannelBad(i)) res->SetGainN(i,0.0);
a64f9843 578
2727462b 579 } // loop over strips
14a74335 580
581}
582
0315d466 583//______________________________________________________________________
fd61217e 584Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
2727462b 585 // Computes the integral of a gaussian using Error Function
586 Float_t sqrt2 = TMath::Sqrt(2.0);
587 Float_t sigm2 = sqrt2*s;
588 Float_t integral;
fd61217e 589
bb7e41dd 590 integral = 0.5 * AliMathBase::ErfFast( (x - av) / sigm2);
2727462b 591 return integral;
0315d466 592}
593//______________________________________________________________________
fd61217e 594void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w,
2727462b 595 Double_t sigma,
596 Double_t inf, Double_t sup,
597 AliITSTableSSD *tav) {
598 // integrate the diffusion gaussian
599 // remind: inf and sup are w-3sigma and w+3sigma
600 // we could define them here instead of passing them
601 // this way we are free to introduce asimmetry
602
603 Double_t a=0.0, b=0.0;
604 Double_t dXCharge1 = 0.0, dXCharge2 = 0.0;
605 // dXCharge1 and 2 are the charge to two neighbouring strips
606 // Watch that we only involve at least two strips
607 // Numbers greater than 2 of strips in a cluster depend on
608 // geometry of the track and delta rays, not charge diffusion!
609
610 Double_t strip = TMath::Floor(w); // closest strip on the left
611
612 if ( TMath::Abs((strip - w)) < 0.5) {
613 // gaussian mean is closer to strip on the left
614 a = inf; // integration starting point
615 if((strip+0.5)<=sup) {
616 // this means that the tail of the gaussian goes beyond
617 // the middle point between strips ---> part of the signal
618 // is given to the strip on the right
619 b = strip + 0.5; // integration stopping point
620 dXCharge1 = F( w, b, sigma) - F(w, a, sigma);
621 dXCharge2 = F( w, sup, sigma) - F(w ,b, sigma);
622 }else {
623 // this means that all the charge is given to the strip on the left
624 b = sup;
625 dXCharge1 = 0.9973; // gaussian integral at 3 sigmas
626 dXCharge2 = 0.0;
627 } // end if
628 dXCharge1 = par * dXCharge1;// normalize by mean of number of carriers
629 dXCharge2 = par * dXCharge2;
630
631 // for the time being, signal is the charge
632 // in ChargeToSignal signal is converted in ADC channel
633 fMapA2->AddSignal(k,(Int_t)strip,dXCharge1);
634 tav->Add(k,(Int_t)strip);
635 if(((Int_t) strip) < (GetNStrips()-1)) {
636 // strip doesn't have to be the last (remind: last=GetNStrips()-1)
637 // otherwise part of the charge is lost
638 fMapA2->AddSignal(k,((Int_t)strip+1),dXCharge2);
639 tav->Add(k,((Int_t)(strip+1)));
640 } // end if
641 }else{
642 // gaussian mean is closer to strip on the right
643 strip++; // move to strip on the rigth
644 b = sup; // now you know where to stop integrating
645 if((strip-0.5)>=inf) {
646 // tail of diffusion gaussian on the left goes left of
647 // middle point between strips
648 a = strip - 0.5; // integration starting point
649 dXCharge1 = F(w, b, sigma) - F(w, a, sigma);
650 dXCharge2 = F(w, a, sigma) - F(w, inf, sigma);
651 }else {
652 a = inf;
653 dXCharge1 = 0.9973; // gaussian integral at 3 sigmas
654 dXCharge2 = 0.0;
655 } // end if
656 dXCharge1 = par * dXCharge1; // normalize by means of carriers
657 dXCharge2 = par * dXCharge2;
658 // for the time being, signal is the charge
659 // in ChargeToSignal signal is converted in ADC channel
660 fMapA2->AddSignal(k,(Int_t)strip,dXCharge1);
661 tav->Add(k,(Int_t)strip);
662 if(((Int_t) strip) > 0) {
663 // strip doesn't have to be the first
664 // otherwise part of the charge is lost
665 fMapA2->AddSignal(k,((Int_t)strip-1),dXCharge2);
666 tav->Add(k,((Int_t)(strip-1)));
667 } // end if
668 } // end if
b0f5e3fc 669}
0315d466 670//______________________________________________________________________
fd61217e 671Int_t AliITSsimulationSSD::NumOfSteps(Double_t x, Double_t y, Double_t z,
2727462b 672 Double_t &dex,Double_t &dey,
673 Double_t &dez){
674 // number of steps
675 // it also returns steps for each coord
676 //AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
677
678 Double_t step = 25E-4;
679 //step = (Double_t) seg->GetStepSize(); // step size (cm)
680 Int_t numOfSteps = (Int_t) (TMath::Sqrt(x*x+y*y+z*z)/step);
681
682 if (numOfSteps < 1) numOfSteps = 1; // one step, at least
683 //numOfSteps=1;
684
685 // we could condition the stepping depending on the incident angle
686 // of the track
687 dex = x/numOfSteps;
688 dey = y/numOfSteps;
689 dez = z/numOfSteps;
690
691 return numOfSteps;
fd61217e 692}
0315d466 693//----------------------------------------------------------------------
c7a4dac0 694void AliITSsimulationSSD::GetList(Int_t label,Int_t hit,Int_t mod,
2727462b 695 AliITSpList *pList,AliITSTableSSD *tav) {
696 // loop over nonzero digits
697 Int_t ix,i;
698 Double_t signal=0.;
699
700 for(Int_t k=0; k<2; k++) {
701 ix=tav->Use(k);
702 while(ix>-1){
703 signal = fMapA2->GetSignal(k,ix);
704 if(signal==0.0) {
705 ix=tav->Use(k);
706 continue;
707 } // end if signal==0.0
708 // check the signal magnitude
709 for(i=0;i<pList->GetNSignals(k,ix);i++){
710 signal -= pList->GetTSignal(k,ix,i);
711 } // end for i
712 // compare the new signal with already existing list
713 if(signal>0)pList->AddSignal(k,ix,label,hit,mod,signal);
714 ix=tav->Use(k);
715 } // end of loop on strips
716 } // end of loop on P/N side
717 tav->Clear();
fd61217e 718}
0315d466 719//----------------------------------------------------------------------
18f63405 720void AliITSsimulationSSD::ChargeToSignal(Int_t module,const AliITSpList *pList) {
2727462b 721 // charge to signal
722 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
723 Float_t threshold = 0.;
724 Int_t size = AliITSdigitSSD::GetNTracks();
725 Int_t * digits = new Int_t[size];
726 Int_t * tracks = new Int_t[size];
727 Int_t * hits = new Int_t[size];
728 Int_t j1;
729 Float_t charges[3] = {0.0,0.0,0.0};
730 Float_t signal;
731 AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
732 AliITSSimuParam* simpar = fDetType->GetSimuParam();
733
734 for(Int_t k=0;k<2;k++){ // both sides (0=Pside, 1=Nside)
735 for(Int_t ix=0;ix<GetNStrips();ix++){ // loop over strips
14a74335 736
737 // if strip is dead -> gain=0
738 if( ((k==0)&&(res->GetGainP(ix)==0)) || ((k==1)&&(res->GetGainN(ix)==0))) continue;
739
3a4139a2 740 signal = fMapA2->GetSignal(k,ix);
14a74335 741 // signal has to be uncalibrated
742 // In real life, gains are supposed to be calculated from calibration runs,
743 // stored in the calibration DB and used in the reconstruction
744 // (see AliITSClusterFinderSSD.cxx)
745 if(k==0) signal /= res->GetGainP(ix);
746 else signal /= res->GetGainN(ix);
747
748 // signal is converted in unit of ADC
8be4e1b1 749 signal = res->GetSSDDEvToADC(signal);
f957a032 750 if(signal>4095.) signal = 4095.;//if exceeding, accumulate last one
14a74335 751
752 // threshold for zero suppression is set on the basis of the noise
753 // A good value is 3*sigma_noise
ced4d9bc 754 if(k==0) threshold = res->GetNoiseP(ix);
755 else threshold = res->GetNoiseN(ix);
756
8be4e1b1 757 threshold *= simpar->GetSSDZSThreshold(); // threshold at 3 sigma noise
ced4d9bc 758
14a74335 759 if(signal < threshold) continue;
ced4d9bc 760 //cout<<signal<<" "<<threshold<<endl;
14a74335 761
762 digits[0] = k;
763 digits[1] = ix;
764 digits[2] = TMath::Nint(signal);
765 for(j1=0;j1<size;j1++)if(j1<pList->GetNEntries()){
766 // only three in digit.
767 tracks[j1] = pList->GetTrack(k,ix,j1);
768 hits[j1] = pList->GetHit(k,ix,j1);
769 }else{
770 tracks[j1] = -3;
771 hits[j1] = -1;
772 } // end for j1
773 // finally add digit
774 aliITS->AddSimDigit(2,0,digits,tracks,hits,charges);
2727462b 775 } // end for ix
776 } // end for k
777 delete [] digits;
778 delete [] tracks;
779 delete [] hits;
c7a4dac0 780}
781//______________________________________________________________________
782void AliITSsimulationSSD::WriteSDigits(AliITSpList *pList){
2727462b 783 // Fills the Summable digits Tree
784 Int_t i,ni,j,nj;
785 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
786
787 pList->GetMaxMapIndex(ni,nj);
788 for(i=0;i<ni;i++)for(j=0;j<nj;j++){
789 if(pList->GetSignalOnly(i,j)>0.0){
790 aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
791 if(GetDebug(4)) cout << "pListSSD: "<<*(pList->GetpListItem(i,j))
792 << endl;
793 } // end if
794 } // end for i,j
795 return;
c7a4dac0 796}
797//______________________________________________________________________
798void AliITSsimulationSSD::FillMapFrompList(AliITSpList *pList){
2727462b 799 // Fills fMap2A from the pList of Summable digits
800 Int_t k,ix;
c7a4dac0 801
2727462b 802 for(k=0;k<2;k++)for(ix=0;ix<GetNStrips();ix++)
803 fMapA2->AddSignal(k,ix,pList->GetSignal(k,ix));
804 return;
b0f5e3fc 805}
57817f7c 806//______________________________________________________________________
807void AliITSsimulationSSD::Print(ostream *os){
2727462b 808 //Standard output format for this class
809
810 //AliITSsimulation::Print(os);
811 *os << fIonE <<",";
812 *os << fDifConst[0] <<","<< fDifConst[1] <<",";
813 *os << fDriftVel[0] <<","<< fDriftVel[1];
814 //*os <<","; fDCS->Print(os);
815 //*os <<","; fMapA2->Print(os);
57817f7c 816}
817//______________________________________________________________________
818void AliITSsimulationSSD::Read(istream *is){
2727462b 819 // Standard output streaming function.
820
821 //AliITSsimulation::Read(is);
822 *is >> fIonE;
823 *is >> fDifConst[0] >> fDifConst[1];
824 *is >> fDriftVel[0] >> fDriftVel[1];
825 //fDCS->Read(is);
826 //fMapA2->Read(is);
57817f7c 827}
828//______________________________________________________________________
829ostream &operator<<(ostream &os,AliITSsimulationSSD &source){
2727462b 830 // Standard output streaming function.
57817f7c 831
2727462b 832 source.Print(&os);
833 return os;
57817f7c 834}
835//______________________________________________________________________
836istream &operator>>(istream &os,AliITSsimulationSSD &source){
2727462b 837 // Standard output streaming function.
57817f7c 838
2727462b 839 source.Read(&os);
840 return os;
57817f7c 841}
c7a4dac0 842//______________________________________________________________________
2727462b 843
844
845