]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSsimulationSPD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
1) Flag all track with kITSupg bit.
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPD.cxx
CommitLineData
c7a4dac0 1/**************************************************************************
20eef074 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**************************************************************************/
c7a4dac0 15
5bfe44ce 16/*
85f5e9c2 17$Id$
5bfe44ce 18*/
f8d9a5b8 19
5bfe44ce 20#include <Riostream.h>
f7a1cc68 21#include <TGeoGlobalMagField.h>
b0f5e3fc 22#include <TH1.h>
9f033001 23#include <TString.h>
e8189707 24#include "AliITS.h"
e869281d 25#include "AliITSdigitSPD.h"
5bfe44ce 26#include "AliITShit.h"
9f033001 27#include "AliITSmodule.h"
c7a4dac0 28#include "AliITSpList.h"
5bfe44ce 29#include "AliITSCalibrationSPD.h"
30#include "AliITSsegmentationSPD.h"
b0f5e3fc 31#include "AliITSsimulationSPD.h"
f77f13c8 32#include "AliLog.h"
5bfe44ce 33#include "AliRun.h"
a4005be7 34#include "AliMagF.h"
bb7e41dd 35#include "AliMathBase.h"
9f033001 36
bb6b39bf 37//#define DEBUG
e8189707 38
fe7d86eb 39using std::endl;
40using std::cout;
b0f5e3fc 41ClassImp(AliITSsimulationSPD)
42////////////////////////////////////////////////////////////////////////
5bfe44ce 43// Version: 1
590d15ee 44// Modified by D. Elia, G.E. Bruno, H. Tydesjo
45// Fast diffusion code by Bjorn S. Nilsen
5bfe44ce 46// March-April 2006
8e2f60b2 47// October 2007: GetCalibrationObjects() removed
b0f5e3fc 48//
5bfe44ce 49// Version: 0
50// Written by Boris Batyunya
51// December 20 1999
21b825a4 52//
5bfe44ce 53//
54// AliITSsimulationSPD is to do the simulation of SPDs.
55//
56////////////////////////////////////////////////////////////////////////
57
c7a4dac0 58//______________________________________________________________________
5bfe44ce 59AliITSsimulationSPD::AliITSsimulationSPD():
60AliITSsimulation(),
61fHis(0),
62fSPDname(),
a4005be7 63fCoupling(),
64fLorentz(kFALSE),
20eef074 65fTanLorAng(0),
66fStrobe(kTRUE),
67fStrobeLenght(4),
68fStrobePhase(-12.5e-9){
69 // Default constructor.
70 // Inputs:
71 // none.
72 // Outputs:
73 // none.
74 // Return:
75 // A default constructed AliITSsimulationSPD class.
76
77 AliDebug(1,Form("Calling default constructor"));
5bfe44ce 78// Init();
b0f5e3fc 79}
c7a4dac0 80//______________________________________________________________________
8ba39da9 81AliITSsimulationSPD::AliITSsimulationSPD(AliITSDetTypeSim *dettyp):
82AliITSsimulation(dettyp),
5bfe44ce 83fHis(0),
84fSPDname(),
a4005be7 85fCoupling(),
86fLorentz(kFALSE),
20eef074 87fTanLorAng(0),
88fStrobe(kTRUE),
89fStrobeLenght(4),
90fStrobePhase(-12.5e-9){
91 // standard constructor
92 // Inputs:
93 // AliITSsegmentation *seg A pointer to the segmentation class
94 // to be used for this simulation
95 // AliITSCalibration *resp A pointer to the responce class to
96 // be used for this simulation
97 // Outputs:
98 // none.
99 // Return:
100 // A default constructed AliITSsimulationSPD class.
101
102 AliDebug(1,Form("Calling standard constructor "));
103 Init();
aacedc3e 104}
105//______________________________________________________________________
5bfe44ce 106void AliITSsimulationSPD::Init(){
20eef074 107 // Initilization
108 // Inputs:
109 // none.
110 // Outputs:
111 // none.
112 // Return:
113 // none.
114 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
115
116 SetModuleNumber(0);
117 SetEventNumber(0);
118 SetMap(new AliITSpList(GetNPixelsZ(),GetNPixelsX()));
119 AliITSSimuParam* simpar = fDetType->GetSimuParam();
120 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
121 Double_t bias = simpar->GetSPDBiasVoltage();
5bfe44ce 122// cout << "Bias Voltage --> " << bias << endl; // dom
20eef074 123 simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),bias);
5bfe44ce 124// set kind of coupling ("old" or "new")
20eef074 125 char opt[20];
126 simpar->GetSPDCouplingOption(opt);
127 char *old = strstr(opt,"old");
128 if (old) {
129 fCoupling=2;
130 } else {
131 fCoupling=1;
132 } // end if
a0a6914c 133 SetLorentzDrift(simpar->GetSPDLorentzDrift());
134 if (fLorentz) SetTanLorAngle(simpar->GetSPDLorentzHoleWeight());
20eef074 135 //SetStrobeGeneration(kFALSE);
136 if (fStrobe) GenerateStrobePhase();
a4005be7 137}
138//______________________________________________________________________
139Bool_t AliITSsimulationSPD::SetTanLorAngle(Double_t WeightHole) {
20eef074 140 // This function set the Tangent of the Lorentz angle.
141 // A weighted average is used for electrons and holes
142 // Input: Double_t WeightHole: wheight for hole: it should be in the range [0,1]
143 // output: Bool_t : kTRUE in case of success
144 //
145 if(!fDetType) {
146 AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set ");
147 return kFALSE;}
148 if(WeightHole<0) {
149 WeightHole=0.;
150 AliWarning("AliITSsimulationSPD::SetTanLorAngle: You have asked for negative Hole weight");
151 AliWarning("AliITSsimulationSPD::SetTanLorAngle: I'm going to use only electrons");
152 }
153 if(WeightHole>1) {
154 WeightHole=1.;
155 AliWarning("AliITSsimulationSPD::SetTanLorAngle: You have asked for weight > 1");
156 AliWarning("AliITSsimulationSPD::SetTanLorAngle: I'm going to use only holes");
157 }
158 Double_t WeightEle=1.-WeightHole;
159 AliITSSimuParam* simpar = fDetType->GetSimuParam();
a0a6914c 160 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
161 if (!fld) AliFatal("The field is not initialized");
162 Double_t bz = fld->SolenoidField();
163 fTanLorAng = TMath::Tan(WeightHole*simpar->LorentzAngleHole(bz) +
164 WeightEle*simpar->LorentzAngleElectron(bz));
20eef074 165 return kTRUE;
c7a4dac0 166}
167//______________________________________________________________________
5bfe44ce 168AliITSsimulationSPD::~AliITSsimulationSPD(){
20eef074 169 // destructor
170 // Inputs:
171 // none.
172 // Outputs:
173 // none.
174 // Return:
175 // none.
176
177 if (fHis) {
178 fHis->Delete();
179 delete fHis;
180 } // end if fHis
b0f5e3fc 181}
c7a4dac0 182//______________________________________________________________________
5bfe44ce 183AliITSsimulationSPD::AliITSsimulationSPD(const
184 AliITSsimulationSPD
7537d03c 185 &s) : AliITSsimulation(s),
186fHis(s.fHis),
187fSPDname(s.fSPDname),
a4005be7 188fCoupling(s.fCoupling),
189fLorentz(s.fLorentz),
20eef074 190fTanLorAng(s.fTanLorAng),
191fStrobe(s.fStrobe),
192fStrobeLenght(s.fStrobeLenght),
193fStrobePhase(s.fStrobePhase){
194 // Copy Constructor
195 // Inputs:
196 // AliITSsimulationSPD &s The original class for which
197 // this class is a copy of
198 // Outputs:
199 // none.
200 // Return:
c7a4dac0 201
b0f5e3fc 202}
c7a4dac0 203//______________________________________________________________________
5bfe44ce 204AliITSsimulationSPD& AliITSsimulationSPD::operator=(const
20eef074 205 AliITSsimulationSPD &s){
206 // Assignment operator
207 // Inputs:
208 // AliITSsimulationSPD &s The original class for which
209 // this class is a copy of
210 // Outputs:
211 // none.
212 // Return:
213
214 if(&s == this) return *this;
215 this->fHis = s.fHis;
216 fCoupling = s.fCoupling;
217 fSPDname = s.fSPDname;
218 fLorentz = s.fLorentz;
219 fTanLorAng = s.fTanLorAng;
220 fStrobe = s.fStrobe;
221 fStrobeLenght = s.fStrobeLenght;
222 fStrobePhase = s.fStrobePhase;
223 return *this;
5bfe44ce 224}
85f5e9c2 225/*
d2f55a22 226//______________________________________________________________________
5bfe44ce 227AliITSsimulation& AliITSsimulationSPD::operator=(const
20eef074 228 AliITSsimulation &s){
229 // Assignment operator
230 // Inputs:
231 // AliITSsimulationSPD &s The original class for which
232 // this class is a copy of
233 // Outputs:
234 // none.
235 // Return:
236
237 if(&s == this) return *this;
238 Error("AliITSsimulationSPD","Not allowed to make a = with "
239 "AliITSsimulationSPD","Using default creater instead");
240
241 return *this;
5bfe44ce 242}
85f5e9c2 243*/
5bfe44ce 244//______________________________________________________________________
245void AliITSsimulationSPD::InitSimulationModule(Int_t module, Int_t event){
20eef074 246 // This function creates maps to build the list of tracks for each
247 // summable digit. Inputs defined by base class.
248 // Inputs:
249 // Int_t module // Module number to be simulated
250 // Int_t event // Event number to be simulated
251 // Outputs:
252 // none
253 // Returns:
254 // none
255
256 AliDebug(1,Form("(module=%d,event=%d)",module,event));
257 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
258 AliITSSimuParam* simpar = fDetType->GetSimuParam();
259 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
260 SetModuleNumber(module);
261 SetEventNumber(event);
262 simpar->SetDistanceOverVoltage(kmictocm*seg->Dy(),simpar->GetSPDBiasVoltage(module));
263 ClearMap();
5bfe44ce 264}
265//_____________________________________________________________________
266void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod,Int_t,
20eef074 267 Int_t event){
268 // This function begins the work of creating S-Digits. Inputs defined
269 // by base class.
270 // Inputs:
271 // AliITSmodule *mod // module
272 // Int_t // not used
273 // Int_t event // Event number
274 // Outputs:
275 // none
276 // Return:
277 // test // test returns kTRUE if the module contained hits
278 // // test returns kFALSE if it did not contain hits
279
280 AliDebug(1,Form("(mod=%p, ,event=%d)",mod,event));
281 if(!(mod->GetNhits())){
282 AliDebug(1,Form("In event %d module %d there are %d hits returning.",
5bfe44ce 283 event, mod->GetIndex(),mod->GetNhits()));
20eef074 284 return;// if module has no hits don't create Sdigits
285 } // end if
286 SetModuleNumber(mod->GetIndex());
287 if (fStrobe) if(event != GetEventNumber()) GenerateStrobePhase();
288 SetEventNumber(event);
289 InitSimulationModule( GetModuleNumber() , event );
290 // HitToSDigit(mod);
291 HitToSDigitFast(mod);
ad7f2bfa 292 if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag()) AddNoisyPixels();
293 if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels();
294
590d15ee 295// cout << "After Remove in SDigitiseModule !!!!!" << endl; // dom
296// cout << "Module " << mod->GetIndex() << " Event " << event << endl; // dom
20eef074 297 WriteSDigits();
298 ClearMap();
5bfe44ce 299}
300//______________________________________________________________________
301void AliITSsimulationSPD::WriteSDigits(){
20eef074 302 // This function adds each S-Digit to pList
303 // Inputs:
304 // none.
305 // Outputs:
306 // none.
307 // Return:
308 // none
309 Int_t ix, nix, iz, niz;
310 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
311
312 AliDebug(1,Form("Writing SDigits for module %d",GetModuleNumber()));
590d15ee 313// cout << "WriteSDigits for module " << GetModuleNumber() << endl; // dom
20eef074 314 GetMap()->GetMaxMapIndex(niz, nix);
315 for(iz=0; iz<niz; iz++)for(ix=0; ix<nix; ix++){
316 if(GetMap()->GetSignalOnly(iz,ix)>0.0){
590d15ee 317// cout << " Signal gt 0 iz ix " << iz << ix << " Module " << GetModuleNumber() << endl; // dom
20eef074 318 aliITS->AddSumDigit(*(GetMap()->GetpListItem(iz,ix)));
5bfe44ce 319 if(AliDebugLevel()>0) {
320 AliDebug(1,Form("%d, %d",iz,ix));
321 cout << *(GetMap()->GetpListItem(iz,ix)) << endl;
20eef074 322 } // end if GetDebug
323 } // end if GetMap()->GetSignalOnly(iz,ix)>0.0
324 } // end for iz,ix
325 return;
3a97c582 326}
327//______________________________________________________________________
328void AliITSsimulationSPD::FinishSDigitiseModule(){
20eef074 329 // This function calls SDigitsToDigits which creates Digits from SDigits
330 // Inputs:
331 // none
332 // Outputs:
333 // none
334 // Return
335 // none
336
337 AliDebug(1,"()");
590d15ee 338// cout << "FinishSDigitiseModule for module " << GetModuleNumber() << endl; // dom
20eef074 339 FrompListToDigits(); // Charge To Signal both adds noise and
340 ClearMap();
341 return;
5bfe44ce 342}
343//______________________________________________________________________
344void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod,Int_t,
20eef074 345 Int_t event){
346 // This function creates Digits straight from the hits and then adds
347 // electronic noise to the digits before adding them to pList
348 // Each of the input variables is passed along to HitToSDigit
349 // Inputs:
350 // AliITSmodule *mod module
351 // Int_t Dummy.
352 // Int_t Dummy
353 // Outputs:
354 // none.
355 // Return:
356 // none.
357
358 if (fStrobe) if(event != GetEventNumber()) GenerateStrobePhase();
1a429a6b 359 AliDebug(1,Form("(mod=%p,,0)",mod));
20eef074 360 // HitToSDigit(mod);
361 InitSimulationModule( mod->GetIndex(), event );
362 HitToSDigitFast(mod);
ad7f2bfa 363
364 if (fDetType->GetSimuParam()->GetSPDAddNoisyFlag()) AddNoisyPixels();
365 if (fDetType->GetSimuParam()->GetSPDRemoveDeadFlag()) RemoveDeadPixels();
590d15ee 366// cout << "After Remove in DigitiseModule in module " << mod->GetIndex() << endl; // dom
20eef074 367 FrompListToDigits();
368 ClearMap();
c7a4dac0 369}
370//______________________________________________________________________
5bfe44ce 371void AliITSsimulationSPD::HitToSDigit(AliITSmodule *mod){
20eef074 372 // Does the charge distributions using Gaussian diffusion charge charing.
373 // Inputs:
374 // AliITSmodule *mod Pointer to this module
375 // Output:
376 // none.
377 // Return:
378 // none.
379 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
380 const Double_t kBunchLenght = 25e-9; // LHC clock
381 TObjArray *hits = mod->GetHits();
382 Int_t nhits = hits->GetEntriesFast();
383 Int_t h,ix,iz,i;
384 Int_t idtrack;
385 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0;
386 Double_t x,y,z,t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
387 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
388 AliITSSimuParam *simpar = fDetType->GetSimuParam();
389 Double_t thick = 0.5*kmictocm*seg->Dy(); // Half Thickness
390 simpar->GetSPDSigmaDiffusionAsymmetry(fda); //
391
392 AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling));
393 if(nhits<=0) return;
394 for(h=0;h<nhits;h++){
395 if(AliDebugLevel()>0) {
5bfe44ce 396 AliDebug(1,Form("Hits, %d", h));
397 cout << *(mod->GetHit(h)) << endl;
20eef074 398 } // end if GetDebug
399 // Check if the hit is inside readout window
400 if (fStrobe)
401 if ((mod->GetHit(h)->GetTOF() < fStrobePhase) ||
402 (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue;
403 if(!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
404 st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
405 if(st>0.0){
406 st = (Double_t)((Int_t)(st/kmictocm)); // number of microns
407 if(st<=1.0) st = 1.0;
408 dt = 1.0/st;
409 for(t=0.0;t<1.0;t+=dt){ // Integrate over t
410 tp = t+0.5*dt;
411 x = x0+x1*tp;
412 y = y0+y1*tp;
413 z = z0+z1*tp;
414 if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
415 //el = res->GeVToCharge((Double_t)(dt*de));
416 el = dt * de / simpar->GetGeVToCharge();
417 if(GetDebug(1)){
418 if(el<=0.0) cout<<"el="<<el<<" dt="<<dt
419 <<" de="<<de<<endl;
420 } // end if GetDebug
421 sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
422 sigx=sig;
423 sigz=sig*fda;
424 if (fLorentz) ld=(y+thick)*fTanLorAng;
425 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
426 } // end for t
427 } else { // st == 0.0 deposit it at this point
428 x = x0;
429 y = y0;
430 z = z0;
431 if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
432 //el = res->GeVToCharge((Double_t)de);
433 el = de / simpar->GetGeVToCharge();
434 sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
435 sigx=sig;
436 sigz=sig*fda;
437 if (fLorentz) ld=(y+thick)*fTanLorAng;
438 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
439 } // end if st>0.0
440
441 } // Loop over all hits h
442
443 // Coupling
444 switch (fCoupling) {
445 default:
446 break;
447 case 1: //case 3:
448 for(i=0;i<GetMap()->GetEntries();i++)
449 if(GetMap()->GetpListItem(i)==0) continue;
450 else{
451 GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
2cad3dbf 452 SetCoupling(iz,ix);
20eef074 453 } // end for i
454 break;
455 case 2: // case 4:
456 for(i=0;i<GetMap()->GetEntries();i++)
457 if(GetMap()->GetpListItem(i)==0) continue;
458 else{
459 GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
2cad3dbf 460 SetCouplingOld(iz,ix);
20eef074 461 } // end for i
462 break;
463 } // end switch
464 if(GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling);
c7a4dac0 465}
466//______________________________________________________________________
5bfe44ce 467void AliITSsimulationSPD::HitToSDigitFast(AliITSmodule *mod){
20eef074 468 // Does the charge distributions using Gaussian diffusion charge charing. // Inputs:
469 // AliITSmodule *mod Pointer to this module
470 // Output:
471 // none.
472 // Return:
473 // none.
474 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
475 const Int_t kn10=10;
476 const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
477 4.325316833e-1,4.869532643e-1,5.130467358e-1,
478 5.674683167e-1,6.602952159e-1,7.833023029e-1,
479 9.255628306e-1};
480 const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
481 7.472567455e-2,3.333567215e-2,3.333567215e-2,
482 7.472567455e-2,1.095431813e-1,1.346333597e-1,
483 1.477621124e-1};
484 const Double_t kBunchLenght = 25e-9; // LHC clock
485 TObjArray *hits = mod->GetHits();
486 Int_t nhits = hits->GetEntriesFast();
487 Int_t h,ix,iz,i;
488 Int_t idtrack;
489 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0;
490 Double_t x,y,z,t,st,el,sig,sigx,sigz,fda;
491 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
492 AliITSSimuParam* simpar = fDetType->GetSimuParam();
493 Double_t thick = 0.5*kmictocm*seg->Dy(); // Half thickness
494 simpar->GetSPDSigmaDiffusionAsymmetry(fda);
5bfe44ce 495// cout << "Half Thickness " << thick << endl; // dom
496// cout << "Diffusion asymm " << fda << endl; // dom
497
20eef074 498 AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling));
499 if(nhits<=0) return;
500 for(h=0;h<nhits;h++){
501 if(AliDebugLevel()>0) {
502 AliDebug(1,Form("Hits, %d", h));
503 cout << *(mod->GetHit(h)) << endl;
504 } // end if GetDebug
505 // Check if the hit is inside readout window
506 if (fStrobe)
507 if ((mod->GetHit(h)->GetTOF() < fStrobePhase) ||
508 (mod->GetHit(h)->GetTOF() > (fStrobePhase+(Double_t)fStrobeLenght*kBunchLenght))) continue;
509 if(!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
510 st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
511 if(st>0.0) for(i=0;i<kn10;i++){ // Integrate over t
512 t = kti[i];
513 x = x0+x1*t;
514 y = y0+y1*t;
515 z = z0+z1*t;
516 if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
517 el = kwi[i]*de/simpar->GetGeVToCharge();
518 if(GetDebug(1)){
519 if(el<=0.0) cout<<"el="<<el<<" kwi["<<i<<"]="<<kwi[i]
520 <<" de="<<de<<endl;
521 } // end if GetDebug
522 sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
523 sigx=sig;
524 sigz=sig*fda;
525 if (fLorentz) ld=(y+thick)*fTanLorAng;
526 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
8e2f60b2 527// cout << "sigx sigz " << sigx << " " << sigz << endl; // dom
20eef074 528 } // end for i // End Integrate over t
529 else { // st == 0.0 deposit it at this point
530 x = x0;
531 y = y0;
532 z = z0;
533 if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
534 el = de / simpar->GetGeVToCharge();
535 sig = simpar->SigmaDiffusion1D(TMath::Abs(thick + y));
536 sigx=sig;
537 sigz=sig*fda;
538 if (fLorentz) ld=(y+thick)*fTanLorAng;
539 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
540 } // end if st>0.0
541
542 } // Loop over all hits h
543
544 // Coupling
545 switch (fCoupling) {
546 default:
547 break;
548 case 1: // case 3:
549 for(i=0;i<GetMap()->GetEntries();i++)
550 if(GetMap()->GetpListItem(i)==0) continue;
551 else{
552 GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
2cad3dbf 553 SetCoupling(iz,ix);
20eef074 554 } // end for i
555 break;
556 case 2: // case 4:
557 for(i=0;i<GetMap()->GetEntries();i++)
558 if(GetMap()->GetpListItem(i)==0) continue;
559 else{
560 GetMap()->GetMapIndex(GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
2cad3dbf 561 SetCouplingOld(iz,ix);
20eef074 562 } // end for i
563 break;
564 } // end switch
565 if(GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling);
21b825a4 566}
c7a4dac0 567//______________________________________________________________________
5bfe44ce 568void AliITSsimulationSPD::SpreadCharge(Double_t x0,Double_t z0,
20eef074 569 Int_t ix0,Int_t iz0,
a4005be7 570 Double_t el,Double_t sig,Double_t ld,
571 Int_t t,Int_t hi){
20eef074 572 // Spreads the charge over neighboring cells. Assume charge is distributed
573 // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
574 // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
575 // if fLorentz=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
576 // Defined this way, the integral over all x and z is el.
577 // Inputs:
578 // Double_t x0 x position of point where charge is liberated
579 // Double_t z0 z position of point where charge is liberated
580 // Int_t ix0 row of cell corresponding to point x0
581 // Int_t iz0 columb of cell corresponding to point z0
582 // Double_t el number of electrons liberated in this step
583 // Double_t sig Sigma difusion for this step (y0 dependent)
584 // Double_t ld lorentz drift in x for this step (y0 dependent)
585 // Int_t t track number
586 // Int_t ti hit track index number
587 // Int_t hi hit "hit" index number
588 // Outputs:
589 // none.
590 // Return:
591 // none.
592 const Int_t knx = 3,knz = 2;
593 const Double_t kRoot2 = 1.414213562; // Sqrt(2).
594 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
595 Int_t ix,iz,ixs,ixe,izs,ize;
596 Float_t x,z;
597 Double_t x1,x2,z1,z2,s,sp;
598 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
599
600
601 if(GetDebug(4)) Info("SpreadCharge","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,"
602 "sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi);
603 if(sig<=0.0) { // if sig<=0 No diffusion to simulate.
604 GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
605 if(GetDebug(2)){
606 cout << "sig<=0.0=" << sig << endl;
607 } // end if GetDebug
608 return;
609 } // end if
610 sp = 1.0/(sig*kRoot2);
611 if(GetDebug(2)){
612 cout << "sig=" << sig << " sp=" << sp << endl;
613 } // end if GetDebug
614 ixs = TMath::Max(-knx+ix0,0);
615 ixe = TMath::Min(knx+ix0,seg->Npx()-1);
616 izs = TMath::Max(-knz+iz0,0);
617 ize = TMath::Min(knz+iz0,seg->Npz()-1);
618 for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
619 seg->DetToLocal(ix,iz,x,z); // pixel center
620 x1 = x;
621 z1 = z;
622 x2 = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper
623 x1 -= 0.5*kmictocm*seg->Dpx(ix); // Lower
624 z2 = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper
625 z1 -= 0.5*kmictocm*seg->Dpz(iz); // Lower
626 x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
627 x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
628 z1 -= z0; // Distance from where track traveled
629 z2 -= z0; // Distance from where track traveled
630 s = 0.25; // Correction based on definision of Erfc
bb7e41dd 631 s *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
20eef074 632 if(GetDebug(3)){
633 cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
634 " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<<
635 " sp*x1="<<sp*x1<<" sp*x2="<<sp*x2<<" s="<<s;
636 } // end if GetDebug
bb7e41dd 637 s *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
20eef074 638 if(GetDebug(3)){
639 cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl;
640 } // end if GetDebug
641 GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
642 } // end for ix, iz
c7a4dac0 643}
644//______________________________________________________________________
5bfe44ce 645void AliITSsimulationSPD::SpreadChargeAsym(Double_t x0,Double_t z0,
20eef074 646 Int_t ix0,Int_t iz0,
647 Double_t el,Double_t sigx,Double_t sigz,
648 Double_t ld,Int_t t,Int_t hi){
649 // Spreads the charge over neighboring cells. Assume charge is distributed
650 // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
651 // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
652 // if fLorentz=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
653 // Defined this way, the integral over all x and z is el.
654 // Inputs:
655 // Double_t x0 x position of point where charge is liberated
656 // Double_t z0 z position of point where charge is liberated
657 // Int_t ix0 row of cell corresponding to point x0
658 // Int_t iz0 columb of cell corresponding to point z0
659 // Double_t el number of electrons liberated in this step
660 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
661 // Double_t sigz Sigma difusion along z for this step (y0 dependent)
662 // Double_t ld lorentz drift in x for this stip (y0 dependent)
663 // Int_t t track number
664 // Int_t ti hit track index number
665 // Int_t hi hit "hit" index number
666 // Outputs:
667 // none.
668 // Return:
669 // none.
670 const Int_t knx = 3,knz = 2;
671 const Double_t kRoot2 = 1.414213562; // Sqrt(2).
672 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
673 Int_t ix,iz,ixs,ixe,izs,ize;
674 Float_t x,z;
675 Double_t x1,x2,z1,z2,s,spx,spz;
676 AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
677
678
679 if(GetDebug(4)) Info("SpreadChargeAsym","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,"
a5a317a9 680 "sigx=%e, sigz=%e, t=%d,i=%d)",x0,z0,ix0,iz0,el,sigx,sigz,t,hi);
20eef074 681 if(sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
682 GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
683 if(GetDebug(2)){
684 cout << "sigx<=0.0=" << sigx << endl;
685 cout << "sigz<=0.0=" << sigz << endl;
686 } // end if GetDebug
687 return;
688 } // end if
689 spx = 1.0/(sigx*kRoot2); spz = 1.0/(sigz*kRoot2);
690 if(GetDebug(2)){
691 cout << "sigx=" << sigx << " spx=" << spx << endl;
692 cout << "sigz=" << sigz << " spz=" << spz << endl;
693 } // end if GetDebug
694 ixs = TMath::Max(-knx+ix0,0);
695 ixe = TMath::Min(knx+ix0,seg->Npx()-1);
696 izs = TMath::Max(-knz+iz0,0);
697 ize = TMath::Min(knz+iz0,seg->Npz()-1);
698 for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
699 seg->DetToLocal(ix,iz,x,z); // pixel center
700 x1 = x;
701 z1 = z;
702 x2 = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper
703 x1 -= 0.5*kmictocm*seg->Dpx(ix); // Lower
704 z2 = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper
705 z1 -= 0.5*kmictocm*seg->Dpz(iz); // Lower
706 x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
707 x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
708 z1 -= z0; // Distance from where track traveled
709 z2 -= z0; // Distance from where track traveled
710 s = 0.25; // Correction based on definision of Erfc
bb7e41dd 711 s *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
20eef074 712 if(GetDebug(3)){
713 cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
714 " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<<
715 " spx*x1="<<spx*x1<<" spx*x2="<<spx*x2<<" s="<<s;
716 } // end if GetDebug
bb7e41dd 717 s *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
20eef074 718 if(GetDebug(3)){
719 cout<<" spz*z1="<<spz*z1<<" spz*z2="<<spz*z2<<" s="<<s<< endl;
720 } // end if GetDebug
721 GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
722 } // end for ix, iz
c7a4dac0 723}
724//______________________________________________________________________
ad7f2bfa 725void AliITSsimulationSPD::RemoveDeadPixels(){
726 // Removes dead pixels on each module (ladder)
727 // This should be called before going from sdigits to digits (FrompListToDigits)
728 Int_t mod = GetModuleNumber();
729 AliITSCalibrationSPD* calObj = (AliITSCalibrationSPD*) fDetType->GetCalibrationModel(mod);
730
731 Int_t nrDead = calObj->GetNrBad();
732 for (Int_t i=0; i<nrDead; i++) {
733 GetMap()->DeleteHit(calObj->GetBadColAt(i), calObj->GetBadRowAt(i));
734 }
735}
736//______________________________________________________________________
737void AliITSsimulationSPD::AddNoisyPixels() {
738 // Adds noisy pixels on each module (ladder)
739 // This should be called before going from sdigits to digits (FrompListToDigits)
740 Int_t mod = GetModuleNumber();
741 AliITSCalibrationSPD* calObj = (AliITSCalibrationSPD*) fDetType->GetSPDNoisyModel(mod);
742
743 Int_t nrNoisy = calObj->GetNrBad();
744 for (Int_t i=0; i<nrNoisy; i++) {
745 // adding 10 times the threshold will for sure make this pixel fire...
746 GetMap()->AddNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), mod, 10*GetThreshold());
747 }
5bfe44ce 748}
749//______________________________________________________________________
590d15ee 750void AliITSsimulationSPD::FrompListToDigits(){
20eef074 751 // add noise and electronics, perform the zero suppression and add the
752 // digit to the list
753 // Inputs:
754 // none.
755 // Outputs:
756 // none.
757 // Return:
758 // none.
759 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
760 Int_t j,ix,iz;
761 Double_t electronics;
762 Double_t sig;
763 const Int_t knmaxtrk=AliITSdigit::GetNTracks();
764 static AliITSdigitSPD dig;
765 AliITSSimuParam *simpar = fDetType->GetSimuParam();
766 if(GetDebug(1)) Info("FrompListToDigits","()");
767 for(iz=0; iz<GetNPixelsZ(); iz++) for(ix=0; ix<GetNPixelsX(); ix++){
5bfe44ce 768// NEW (for the moment plugged by hand, in the future possibly read from Data Base)
769// here parametrize the efficiency of the pixel along the row for the test columns (1,9,17,25)
770// if(iz==1 || iz == 9 || iz == 17 || iz == 25) {
771// Double_t eff,p1=0.,p2=0.;
772// Double_t x=ix;
773// switch (iz) {
774// case 1: p1=0.63460;p2=0.42438E-01;break;
775// case 9: p1=0.41090;p2=0.75914E-01;break;
776// case 17: p1=0.31883;p2=0.91502E-01;break;
777// case 25: p1=0.48828;p2=0.57975E-01;break;
778// } // end switch
779// eff=1.-p1*exp(-p2*x);
780// if (gRandom->Rndm() >= eff) continue;
781// } // end if
2ae37d58 782// END parametrize the efficiency
783//
20eef074 784 electronics = simpar->ApplySPDBaselineAndNoise();
785 UpdateMapNoise(ix,iz,electronics);
786 //
787 // Apply Threshold and write Digits.
788 sig = GetMap()->GetSignalOnly(iz,ix);
789 FillHistograms(ix,iz,sig+electronics);
790 if(GetDebug(3)){
791 cout<<sig<<"+"<<electronics<<">threshold("<<ix<<","<<iz
792 <<")="<<GetThreshold() <<endl;
793 } // end if GetDebug
ad7f2bfa 794 // if (sig+electronics <= GetThreshold()) continue;
795 if (GetMap()->GetSignal(iz,ix) <= GetThreshold()) continue;
20eef074 796 dig.SetCoord1(iz);
797 dig.SetCoord2(ix);
798 dig.SetSignal(1);
590d15ee 799
800// dig.SetSignalSPD((Int_t) GetMap()->GetSignal(iz,ix));
20eef074 801 Double_t aSignal = GetMap()->GetSignal(iz,ix);
802 if (TMath::Abs(aSignal)>2147483647.0) {
803 //PH 2147483647 is the max. integer
804 //PH This apparently is a problem which needs investigation
805 AliWarning(Form("Too big or too small signal value %f",aSignal));
806 aSignal = TMath::Sign((Double_t)2147483647,aSignal);
807 }
808 dig.SetSignalSPD((Int_t)aSignal);
809
810 for(j=0;j<knmaxtrk;j++){
811 if (j<GetMap()->GetNEntries()) {
812 dig.SetTrack(j,GetMap()->GetTrack(iz,ix,j));
813 dig.SetHit(j,GetMap()->GetHit(iz,ix,j));
814 }else { // Default values
815 dig.SetTrack(j,-3);
816 dig.SetHit(j,-1);
817 } // end if GetMap()
818 } // end for j
819 if(GetDebug(3)){
820 cout<<iz<<","<<ix<<","<<*(GetMap()->GetpListItem(iz,ix))<<endl;
821 } // end if GetDebug
822 aliITS->AddSimDigit(0,&dig);
ad7f2bfa 823 // simulate fo signal response for this pixel hit:
824 fDetType->ProcessSPDDigitForFastOr(fModule, dig.GetCoord1(), dig.GetCoord2());
20eef074 825 } // for ix/iz
c7a4dac0 826}
827//______________________________________________________________________
5bfe44ce 828void AliITSsimulationSPD::CreateHistograms(){
20eef074 829 // create 1D histograms for tests
830 // Inputs:
831 // none.
832 // Outputs:
833 // none.
834 // Return:
835 // none.
836
837 if(GetDebug(1)) Info("CreateHistograms","create histograms");
838
839 fHis = new TObjArray(GetNPixelsZ());
840 fSPDname="spd_";
841 for(Int_t i=0;i<GetNPixelsZ();i++) {
842 Char_t pixelz[4];
ce189d7c 843 snprintf(pixelz,3,"%d",i);
20eef074 844 fSPDname.Append(pixelz);
845 fHis->AddAt(new TH1F(fSPDname.Data(),"SPD maps",
846 GetNPixelsX(),0.,(Double_t)GetNPixelsX()),i);
847 } // end for i
c7a4dac0 848}
849//______________________________________________________________________
5bfe44ce 850void AliITSsimulationSPD::FillHistograms(Int_t ix,Int_t iz,Double_t v){
20eef074 851 // Fill the histogram
852 // Inputs:
853 // none.
854 // Outputs:
855 // none.
856 // Return:
857 // none.
858
859 if(!GetHistArray()) return; // Only fill if setup.
860 if(GetDebug(2)) Info("FillHistograms","fill histograms");
861 GetHistogram(iz)->Fill(ix,v);
c7a4dac0 862}
863//______________________________________________________________________
5bfe44ce 864void AliITSsimulationSPD::ResetHistograms(){
20eef074 865 // Reset histograms for this detector
866 // Inputs:
867 // none.
868 // Outputs:
869 // none.
870 // Return:
871 // none.
872
873 if(!GetHistArray()) return; // Only fill if setup.
874 if(GetDebug(2)) Info("FillHistograms","fill histograms");
875 for ( int i=0;i<GetNPixelsZ();i++ ) {
876 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();
877 } // end for i
21b825a4 878}
5bfe44ce 879
c7a4dac0 880//______________________________________________________________________
2cad3dbf 881void AliITSsimulationSPD::SetCoupling(Int_t col, Int_t row) {
20eef074 882 // Take into account the coupling between adiacent pixels.
883 // The parameters probcol and probrow are the probability of the
884 // signal in one pixel shared in the two adjacent pixels along
885 // the column and row direction, respectively.
886 // Note pList is goten via GetMap() and module is not need any more.
887 // Otherwise it is identical to that coded by Tiziano Virgili (BSN).
888 //Begin_Html
889 /*
890 <img src="picts/ITS/barimodel_3.gif">
891 </pre>
892 <br clear=left>
893 <font size=+2 color=red>
894 <a href="mailto:tiziano.virgili@cern.ch"></a>.
895 </font>
896 <pre>
897 */
898 //End_Html
899 // Inputs:
900 // Int_t col z cell index
901 // Int_t row x cell index
20eef074 902 // Outputs:
903 // none.
904 // Return:
905 // none.
906 Int_t j1,j2,flag=0;
907 Double_t pulse1,pulse2;
908 Double_t couplR=0.0,couplC=0.0;
909 Double_t xr=0.;
910
911 GetCouplings(couplC,couplR);
2cad3dbf 912 if(GetDebug(3)) Info("SetCoupling","(col=%d,row=%d) "
20eef074 913 "Calling SetCoupling couplC=%e couplR=%e",
2cad3dbf 914 col,row,couplC,couplR);
20eef074 915 j1 = col;
916 j2 = row;
917 pulse1 = GetMap()->GetSignalOnly(col,row);
918 pulse2 = pulse1;
919 for (Int_t isign=-1;isign<=1;isign+=2){// loop in col direction
920 do{
921 j1 += isign;
922 xr = gRandom->Rndm();
923 if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplC)){
924 j1 = col;
925 flag = 1;
926 }else{
2cad3dbf 927 UpdateMapNoise(row,j1,pulse1);
20eef074 928 // flag = 0;
929 flag = 1; // only first next!!
930 } // end if
931 } while(flag == 0);
932 // loop in row direction
933 do{
934 j2 += isign;
935 xr = gRandom->Rndm();
936 if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplR)){
937 j2 = row;
938 flag = 1;
939 }else{
2cad3dbf 940 UpdateMapNoise(j2,col,pulse2);
20eef074 941 // flag = 0;
942 flag = 1; // only first next!!
943 } // end if
944 } while(flag == 0);
945 } // for isign
f8d9a5b8 946}
947//______________________________________________________________________
2cad3dbf 948void AliITSsimulationSPD::SetCouplingOld(Int_t col, Int_t row) {
20eef074 949 // Take into account the coupling between adiacent pixels.
950 // The parameters probcol and probrow are the fractions of the
951 // signal in one pixel shared in the two adjacent pixels along
952 // the column and row direction, respectively.
953 //Begin_Html
954 /*
955 <img src="picts/ITS/barimodel_3.gif">
956 </pre>
957 <br clear=left>
958 <font size=+2 color=red>
959 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
960 </font>
961 <pre>
962 */
963 //End_Html
964 // Inputs:
965 // Int_t col z cell index
966 // Int_t row x cell index
20eef074 967 // Int_t module module number
968 // Outputs:
969 // none.
970 // Return:
971 // none.
972 Int_t j1,j2,flag=0;
973 Double_t pulse1,pulse2;
974 Double_t couplR=0.0,couplC=0.0;
975
976 GetCouplings(couplC,couplR);
977
978 // Debugging ...
5bfe44ce 979// cout << "Threshold --> " << GetThreshold() << endl; // dom
8e2f60b2 980// cout << "Couplings --> " << couplC << " " << couplR << endl; // dom
5bfe44ce 981
2cad3dbf 982 if(GetDebug(3)) Info("SetCouplingOld","(col=%d,row=%d) "
20eef074 983 "Calling SetCoupling couplC=%e couplR=%e",
2cad3dbf 984 col,row,couplC,couplR);
20eef074 985 for (Int_t isign=-1;isign<=1;isign+=2){// loop in col direction
986 pulse1 = GetMap()->GetSignalOnly(col,row);
987 pulse2 = pulse1;
988 j1 = col;
989 j2 = row;
990 do{
991 j1 += isign;
992 pulse1 *= couplC;
993 if ((j1<0)||(j1>GetNPixelsZ()-1)||(pulse1<GetThreshold())){
994 pulse1 = GetMap()->GetSignalOnly(col,row);
995 j1 = col;
996 flag = 1;
997 }else{
2cad3dbf 998 UpdateMapNoise(row,j1,pulse1);
20eef074 999 // flag = 0;
1000 flag = 1; // only first next !!
1001 } // end if
1002 } while(flag == 0);
1003 // loop in row direction
1004 do{
1005 j2 += isign;
1006 pulse2 *= couplR;
1007 if((j2<0)||(j2>(GetNPixelsX()-1))||(pulse2<GetThreshold())){
1008 pulse2 = GetMap()->GetSignalOnly(col,row);
1009 j2 = row;
1010 flag = 1;
1011 }else{
2cad3dbf 1012 UpdateMapNoise(j2,col,pulse2);
20eef074 1013 // flag = 0;
1014 flag = 1; // only first next!!
1015 } // end if
1016 } while(flag == 0);
1017 } // for isign
21b825a4 1018}
20eef074 1019//______________________________________________________________________
1020void AliITSsimulationSPD::GenerateStrobePhase()
1021{
1022 // Generate randomly the strobe
1023 // phase w.r.t to the LHC clock
1024 // Done once per event
1025 const Double_t kBunchLenght = 25e-9; // LHC clock
1026 fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
1027 (Double_t)fStrobeLenght*kBunchLenght+
1028 kBunchLenght/2;
1029}
1030