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