Transition to NewIO
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPDdubna.cxx
CommitLineData
f74211b0 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$Log$
88cb7938 18Revision 1.3.6.2 2003/05/19 16:24:40 hristov
19Merging with v3-09-09 (partial, AliITCclustererV2 is excluded from compilation and has to be adapted to NewIO, AliITStrackerV2 needs additional changes)
20
21Revision 1.10 2003/01/22 01:58:38 nilsen
22Small patch to allow new Dubna modle to reproduce BS model.
23
99df08fd 24Revision 1.9 2002/10/22 14:45:45 alibrary
25Introducing Riostream.h
26
4ae5bbc4 27Revision 1.8 2002/10/14 14:57:08 hristov
28Merging the VirtualMC branch to the main development branch (HEAD)
29
b9d0a01d 30Revision 1.3.8.2 2002/10/14 13:14:08 hristov
31Updating VirtualMC to v3-09-02
32
33Revision 1.7 2002/09/09 17:23:28 nilsen
34Minor changes in support of changes to AliITSdigitS?D class'.
35
ee86d557 36Revision 1.6 2002/08/21 22:09:58 nilsen
37Updated SPD simulation with difusion effects. ReWritten Hit to SDigits
38code.
39
f74211b0 40*/
4ae5bbc4 41#include <Riostream.h>
409f8c84 42#include <TRandom.h>
43#include <TH1.h>
44#include <TMath.h>
45#include <TString.h>
46#include <TParticle.h>
47
409f8c84 48#include "AliRun.h"
49#include "AliITS.h"
50#include "AliITShit.h"
51#include "AliITSdigit.h"
52#include "AliITSmodule.h"
53#include "AliITSMapA2.h"
2cc6b29a 54#include "AliITSpList.h"
409f8c84 55#include "AliITSsimulationSPDdubna.h"
f74211b0 56#include "AliITSsegmentationSPD.h"
57#include "AliITSresponseSPDdubna.h"
409f8c84 58
f74211b0 59//#define DEBUG
409f8c84 60
61ClassImp(AliITSsimulationSPDdubna)
62////////////////////////////////////////////////////////////////////////
63// Version: 0
64// Written by Boris Batyunya
65// December 20 1999
66//
67// AliITSsimulationSPDdubna is the simulation of SPDs
2cc6b29a 68//______________________________________________________________________
69
70
71AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(){
72 // constructor
73
74 fResponse = 0;
75 fSegmentation = 0;
76 fMapA2 = 0;
77 fpList = 0;
78 fModule = 0;
79 fEvent = 0;
80 fHis = 0;
81 fNoise = 0.;
82 fBaseline = 0.;
83 fNPixelsZ = 0;
84 fNPixelsX = 0;
409f8c84 85}
2cc6b29a 86//______________________________________________________________________
87AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg,
88 AliITSresponse *resp){
89 // standard constructor
f74211b0 90 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
409f8c84 91
2cc6b29a 92 fHis = 0;
93 fResponse = resp;
94 fSegmentation = seg;
95 fModule = 0;
96 fEvent = 0;
409f8c84 97
99df08fd 98 fNPixelsZ=GetSeg()->Npz();
99 fNPixelsX=GetSeg()->Npx();
409f8c84 100
99df08fd 101 GetResp()->GetNoiseParam(fNoise,fBaseline);
102 GetResp()->SetDistanceOverVoltage(kmictocm*GetSeg()->Dy(),50.0);
409f8c84 103
99df08fd 104// fMapA2 = new AliITSMapA2(GetSeg());
f74211b0 105 fMapA2 = 0;
409f8c84 106
2cc6b29a 107 fpList = new AliITSpList(fNPixelsZ+1,fNPixelsX+1);
409f8c84 108
2cc6b29a 109}
110//______________________________________________________________________
111AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna(){
112 // destructor
409f8c84 113
f74211b0 114 if(fMapA2) delete fMapA2;
409f8c84 115
2cc6b29a 116 if (fHis) {
117 fHis->Delete();
118 delete fHis;
119 } // end if fHis
409f8c84 120}
2cc6b29a 121//______________________________________________________________________
122AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const
123 AliITSsimulationSPDdubna
124 &source){
125 // Copy Constructor
126 if(&source == this) return;
127 this->fMapA2 = source.fMapA2;
128 this->fNoise = source.fNoise;
129 this->fBaseline = source.fBaseline;
130 this->fNPixelsX = source.fNPixelsX;
131 this->fNPixelsZ = source.fNPixelsZ;
132 this->fHis = source.fHis;
133 return;
134}
135//______________________________________________________________________
136AliITSsimulationSPDdubna& AliITSsimulationSPDdubna::operator=(const
137 AliITSsimulationSPDdubna &source){
138 // Assignment operator
139 if(&source == this) return *this;
140 this->fMapA2 = source.fMapA2;
141 this->fNoise = source.fNoise;
142 this->fBaseline = source.fBaseline;
143 this->fNPixelsX = source.fNPixelsX;
144 this->fNPixelsZ = source.fNPixelsZ;
145 this->fHis = source.fHis;
146 return *this;
147}
148//______________________________________________________________________
149void AliITSsimulationSPDdubna::InitSimulationModule(Int_t module, Int_t event){
150 // This function creates maps to build the list of tracks for each
151 // summable digit.
152 //
153 // Inputs:
154 // Int_t module // Module number to be simulated
155 // Int_t event // Event number to be simulated
156 //
157 // Outputs:
158 // none
159 //
160 // Returns:
161 // none
409f8c84 162
2cc6b29a 163 fModule = module;
164 fEvent = event;
f74211b0 165// fMapA2->ClearMap();
2cc6b29a 166 fpList->ClearMap();
167}
168//_____________________________________________________________________
169void AliITSsimulationSPDdubna::SDigitiseModule(AliITSmodule *mod, Int_t mask,
170 Int_t event){
171 // This function begins the work of creating S-Digits
172 //
173 // Inputs:
174 // AliITSmodule *mod // module
175 // Int_t mask // mask to be applied to the module
176 //
177 // Outputs:
178 // none
179 //
180 // Return:
181 // test // test returns kTRUE if the module contained hits
182 // // test returns kFALSE if it did not contain hits
409f8c84 183
2cc6b29a 184 Int_t module = 0;
409f8c84 185
2cc6b29a 186 if(!(mod->GetNhits())) return;// if module has no hits don't create Sdigits
187 fModule = mod->GetIndex();
188 HitToSDigit(mod, module, mask, fpList);
189 WriteSDigits(fpList);
f74211b0 190// fMapA2->ClearMap();
2cc6b29a 191 fpList->ClearMap();
409f8c84 192}
2cc6b29a 193//______________________________________________________________________
194void AliITSsimulationSPDdubna::WriteSDigits(AliITSpList *pList){
195 // This function adds each S-Digit to pList
196 //
197 // Inputs:
198 // AliITSpList *pList
199 //
200 // Outputs:
201 // none
202 //
203 // Return:
204 // none
7d50ea21 205 Int_t ix, nix, iz, niz;
2cc6b29a 206 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
207
7d50ea21 208 pList->GetMaxMapIndex(niz, nix);
99df08fd 209 for(iz=0; iz<niz-1; iz++)for(ix=0; ix<nix-1; ix++){
210 if(pList->GetSignalOnly(iz+1,ix+1)>0.0){
211 aliITS->AddSumDigit(*(pList->GetpListItem(iz+1,ix+1)));
f74211b0 212#ifdef DEBUG
213 cout <<"SDigits " << iz << "," << ix << "," <<
99df08fd 214 *(pList->GetpListItem(iz+1,ix+1)) << endl;
f74211b0 215#endif
2cc6b29a 216 } // end if pList
7d50ea21 217 } // end for iz,ix
2cc6b29a 218 return;
219}
220//______________________________________________________________________
221void AliITSsimulationSPDdubna::FinishSDigitiseModule(){
222 // This function calls SDigitsToDigits which creates Digits from SDigits
223 //
224 // Inputs:
225 // none
226 //
227 // Outputs:
228 // none
229 // Return
230 // none
409f8c84 231
2cc6b29a 232 SDigitsToDigits(fModule, fpList);
233 return;
234}
235//______________________________________________________________________
236void AliITSsimulationSPDdubna::SDigitsToDigits(Int_t module,
237 AliITSpList *pList){
238 // This function adds electronic noise to the S-Digits and then adds them
239 // to a new pList
240 //
241 // Inputs:
242 // Int_t module // module number
243 // AliITSpList *pList // pList
244 //
245 // Outputs:
246 // pList is passed along to the functions ChargeToSignal and GetList
247 //
248 // Return:
249 // none
409f8c84 250
2cc6b29a 251 fModule = module;
252 ChargeToSignal(pList); // Charge To Signal both adds noise and
f74211b0 253// fMapA2->ClearMap();
2cc6b29a 254 pList->ClearMap();
409f8c84 255}
2cc6b29a 256//______________________________________________________________________
257void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module,
258 Int_t dummy){
259 // This function creates Digits straight from the hits and then adds
260 // electronic noise to the digits before adding them to pList
261 //
262 // Inputs:
263 // AliITSmodule *mod // module
264 // Int_t module // module number Dummy.
265 // Int_t dummy
266 //
267 // Outputs:
268 // Each of the input variables is passed along to HitToSDigit
269 //
270 // Return:
271 // none
409f8c84 272
2cc6b29a 273 fModule = mod->GetIndex(); //This calls the module for HitToSDigit
274 HitToSDigit(mod,fModule, dummy, fpList);
275 ChargeToSignal(fpList);
f74211b0 276// fMapA2->ClearMap();
2cc6b29a 277 fpList->ClearMap();
278}
279//______________________________________________________________________
7d50ea21 280void AliITSsimulationSPDdubna::UpdateMapSignal(Int_t iz, Int_t ix, Int_t trk,
2cc6b29a 281 Int_t ht, Int_t module,
282 Double_t signal,
283 AliITSpList *pList){
284 // This function adds a signal to the pList from the pList class
285 //
286 // Inputs:
7d50ea21 287 // Int_t iz // row number
288 // Int_t ix // column number
2cc6b29a 289 // Int_t trk // track number
290 // Int_t ht // hit number
291 // Double_t signal // signal strength
292 // AliITSpList *pList // pList
293 //
294 // Outputs:
295 // All of the inputs are passed to AliITSpList::AddSignal
296 // Int_t ix // row number
297 // Int_t iz // column number
298 // Double_t sig // signal strength
299 // // These three variables are defined to preserve the
300 // // assignments used in the function AliITSMapA2::AddSignal
301 //
302 // Return:
303 // none
2cc6b29a 304
f74211b0 305// fMapA2->AddSignal(iz, ix, signal);
99df08fd 306 pList->AddSignal(iz+1,ix+1, trk, ht, fModule, signal);
2cc6b29a 307}
308//______________________________________________________________________
7d50ea21 309void AliITSsimulationSPDdubna::UpdateMapNoise(Int_t iz,
310 Int_t ix, Int_t fModule,
2cc6b29a 311 Double_t sig, Float_t noise,
312 AliITSpList *pList){
313 // This function adds noise to data in the MapA2 as well as the pList
314 //
315 // Inputs:
7d50ea21 316 // Int_t iz // row number
317 // Int_t ix // column number
2cc6b29a 318 // Int_t mod // module number
319 // Double_t sig // signal strength
320 // Double_t noise // electronic noise generated by ChargeToSignal
321 // AliITSpList *pList // pList
322 //
323 // Outputs:
324 // All of the inputs are passed to AliITSMapA2::AddSignal or
325 // AliITSpList::AddNoise
326 //
327 // Return:
328 // none
409f8c84 329
f74211b0 330// fMapA2->AddSignal(iz, ix, noise);
99df08fd 331 pList->AddNoise(iz+1,ix+1, fModule, noise);
2cc6b29a 332}
333//______________________________________________________________________
334void AliITSsimulationSPDdubna::HitToDigit(AliITSmodule *mod, Int_t module,
335 Int_t dummy){
336 DigitiseModule(mod, module, dummy);
337}
338//______________________________________________________________________
339void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod, Int_t module,
f74211b0 340 Int_t dummy,AliITSpList *pList){
341 // Does the charge distributions using Gaussian diffusion charge charing.
342 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
343 TObjArray *hits = mod->GetHits();
344 Int_t nhits = hits->GetEntriesFast();
345 Int_t h,ix,iz;
346 Int_t idtrack;
347 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
348 Double_t x,y,z,t,tp,st,dt=0.2,el,sig;
349 Double_t thick = kmictocm*GetSeg()->Dy();
350
351 if(nhits<=0) return;
352 for(h=0;h<nhits;h++){
353#ifdef DEBUG
354 cout << "Hits=" << h << "," << *(mod->GetHit(h)) << endl;
355#endif
356 if(mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)){
357 st =TMath::Sqrt(x1*x1+y1*y1+z1*z1);
99df08fd 358 if(st>0.0){
359 st = (Double_t)((Int_t)(1.0E+04*st)); // number of microns
360 if(st<=0.0) st = 1.0;
361 dt = 1.0/st;
362 for(t=0;t<1.0;t+=dt){ // Integrate over t
363 tp = t+0.5*dt;
364 el = GetResp()->GeVToCharge((Float_t)(dt*de));
f74211b0 365#ifdef DEBUG
99df08fd 366 if(el<=0.0) cout << "el="<<el<<" dt="<<dt<<" de="<<de<<endl;
f74211b0 367#endif
99df08fd 368 x = x0+x1*tp;
369 y = y0+y1*tp;
370 z = z0+z1*tp;
371 GetSeg()->LocalToDet(x,z,ix,iz);
372 sig = GetResp()->SigmaDiffusion1D(thick + y);
373 SpreadCharge(x,y,z,ix,iz,el,sig,idtrack,
374 mod->GetHitTrackIndex(h),h,mod->GetIndex());
375 } // end for t
376 } else { // st == 0.0 deposit it at this point
377 el = GetResp()->GeVToCharge((Float_t)de);
378 x = x0;
379 y = y0;
380 z = z0;
f74211b0 381 GetSeg()->LocalToDet(x,z,ix,iz);
382 sig = GetResp()->SigmaDiffusion1D(thick + y);
383 SpreadCharge(x,y,z,ix,iz,el,sig,
384 idtrack,mod->GetHitTrackIndex(h),h,mod->GetIndex());
99df08fd 385 } // end if st>0.0
386 }} // Loop over all hits h
387}/*
388//______________________________________________________________________
389void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod, Int_t module,
390 Int_t dummy,AliITSpList *pList){
391 // Does the charge distributions using Gaussian diffusion charge charing.
392 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
393 TObjArray *hits = mod->GetHits();
394 Int_t nhits = hits->GetEntriesFast();
395 Int_t h,ix,iz,i,n;
396 Int_t idtrack;
397 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
398 Double_t x,y,z,*ta,t,tp,st,dt=0.2,el,sig;
399 Double_t thick = kmictocm*GetSeg()->Dy();
400
401 if(nhits<=0) return;
402 for(h=0;h<nhits;h++){
403#ifdef DEBUG
404 cout << "Hits=" << h << "," << *(mod->GetHit(h)) << endl;
405#endif
406 if(mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)){
407 st =TMath::Sqrt(x1*x1+y1*y1+z1*z1);
408 if(st>0.0){
409 st =TMath::Sqrt(x1*x1+y1*y1+z1*z1)*(ta[i+1]-ta[i]);
410 ta = CreateFindCellEdges(x0,x1,z0,z1,n);
411 for(i=0;i<n-1;i++){
412 dt = TMath::Min((1.0E-4)/st,);
413 for(t=ta[i];t<ta[i+1];t+=dt){ // Integrate over t
414 tp = t+0.5*dt;
415 el = GetResp()->GeVToCharge((Float_t)(dt*de));
416#ifdef DEBUG
417 if(el<=0.0) cout << "el="<<el<<" dt="<<dt<<" de="<<de<<endl;
418#endif
419 x = x0+x1*tp;
420 y = y0+y1*tp;
421 z = z0+z1*tp;
422 GetSeg()->LocalToDet(x,z,ix,iz);
423 sig = GetResp()->SigmaDiffusion1D(thick + y);
424 SpreadCharge(x,y,z,ix,iz,el,sig,idtrack,
425 mod->GetHitTrackIndex(h),h,mod->GetIndex());
426 } // end for t[i]
427 delete[] t;
f74211b0 428 } else { // st == 0.0 deposit it at this point
429 el = GetResp()->GeVToCharge((Float_t)de);
430 x = x0;
431 y = y0;
432 z = z0;
433 GetSeg()->LocalToDet(x,z,ix,iz);
434 sig = GetResp()->SigmaDiffusion1D(thick + y);
435 SpreadCharge(x,y,z,ix,iz,el,sig,
436 idtrack,mod->GetHitTrackIndex(h),h,mod->GetIndex());
437 } // end if st>0.0
438 }} // Loop over all hits h
99df08fd 439 }*/
f74211b0 440//______________________________________________________________________
441void AliITSsimulationSPDdubna::SpreadCharge(Double_t x0,Double_t y0,
442 Double_t z0,Int_t ix0,Int_t iz0,
443 Double_t el,Double_t sig,Int_t t,
444 Int_t ti,Int_t hi,Int_t mod){
445 // Spreads the charge over neighboring cells. Assume charge is distributed
446 // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
447 // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
448 // Defined this way, the integral over all x and z is el.
449 const Int_t knx = 3,knz = 2;
450 const Double_t kRoot2 = 1.414213562; // Sqrt(2).
451 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
452 Int_t ix,iz,ixs,ixe,izs,ize;
453 Float_t x,z;
454 Double_t x1,x2,z1,z2,s,sp;
455
456 if(sig<=0.0) {
99df08fd 457 fpList->AddSignal(iz0+1,ix0+1,t,hi,mod,el);
f74211b0 458 return;
459 } // end if
460 sp = 1.0/(sig*kRoot2);
461#ifdef DEBUG
462 cout << "sig=" << sig << " sp=" << sp << endl;
463#endif
464 ixs = TMath::Max(-knx+ix0,0);
465 ixe = TMath::Min(knx+ix0,GetSeg()->Npx()-1);
466 izs = TMath::Max(-knz+iz0,0);
467 ize = TMath::Min(knz+iz0,GetSeg()->Npz()-1);
468 for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
469 GetSeg()->DetToLocal(ix,iz,x,z); // pixel center
470 x1 = x;
471 z1 = z;
472 x2 = x1 + 0.5*kmictocm*GetSeg()->Dpx(ix); // Upper
473 x1 -= 0.5*kmictocm*GetSeg()->Dpx(ix); // Lower
474 z2 = z1 + 0.5*kmictocm*GetSeg()->Dpz(iz); // Upper
475 z1 -= 0.5*kmictocm*GetSeg()->Dpz(iz); // Lower
476 x1 -= x0; // Distance from where track traveled
477 x2 -= x0; // Distance from where track traveled
478 z1 -= z0; // Distance from where track traveled
479 z2 -= z0; // Distance from where track traveled
480 s = 0.25; // Correction based on definision of Erfc
481 s *= TMath::Erfc(sp*x1) - TMath::Erfc(sp*x2);
482#ifdef DEBUG
483 cout << "el=" << el << " ix0=" << ix0 << " ix=" << ix << " x0="<< x <<
484 " iz0=" << iz0 << " iz=" << iz << " z0=" << z <<
485 " sp*x1=" << sp*x1 <<" sp*x2=" << sp*x2 << " s=" << s;
486#endif
487 s *= TMath::Erfc(sp*z1) - TMath::Erfc(sp*z2);
488#ifdef DEBUG
489 cout << " sp*z1=" << sp*z1 <<" sp*z2=" << sp*z2 << " s=" << s << endl;
490#endif
99df08fd 491 fpList->AddSignal(iz+1,ix+1,t,hi,mod,s*el);
f74211b0 492 } // end for ix, iz
493}
494//______________________________________________________________________
99df08fd 495Double_t *AliITSsimulationSPDdubna::CreateFindCellEdges(Double_t x0,Double_t x1,
496 Double_t z0,Double_t z1,Int_t &n){
497 // Note: This function is a potensial source for a memory leak. The memory
498 // pointed to in its return, must be deleted.
499 // Inputs:
500 // Double_t x0 The starting location of the track step in x
501 // Double_t x1 The distance allong x for the track step
502 // Double_t z0 The starting location of the track step in z
503 // Double_t z1 The distance allong z for the track step
504 // Output:
505 // Int)t &n The size of the array returned. Minimal n=2.
506 // Return:
507 // The pointer to the array of track steps.
508 Int_t ix0,ix1,ix,iz0,iz1,iz,i;
509 Double_t x,z,lx,ux,lz,uz,a,b,c,d;
510 Double_t *t;
511
512 GetSeg()->LocalToDet(x0,z0,ix0,iz0);
513 GetSeg()->LocalToDet(x1,z1,ix1,iz1);
514 n = 2 + TMath::Abs(ix1-ix0) + TMath::Abs(iz1-iz0);
515 t = new Double_t[n];
516 t[0] = 0.0;
517 t[n-1] = 1.0;
518 x = x0;
519 z = z0;
520 for(i=1;i<n-1;i++){
521 GetSeg()->LocalToDet(x,z,ix,iz);
522 GetSeg()->CellBoundries(ix,iz,lx,ux,lz,uz);
523 a = (lx-x0)/x1;
524 if(a<=t[i-1]) a = 1.0;
525 b = (ux-x0)/x1;
526 if(b<=t[i-1]) b = 1.0;
527 c = (lz-z0)/z1;
528 if(c<=t[i-1]) c = 1.0;
529 d = (uz-z0)/z1;
530 if(d<=t[i-1]) d = 1.0;
531 t[i] = TMath::Min(TMath::Min(TMath::Min(a,b),c),d);
532 x = x0+x1*(t[i]*1.00000001);
533 z = z0+z1*(t[i]*1.00000001);
534 i++;
535 } // end for i
536 return t;
537}
538//______________________________________________________________________
f74211b0 539void AliITSsimulationSPDdubna::HitToSDigitOld(AliITSmodule *mod, Int_t module,
2cc6b29a 540 Int_t dummy, AliITSpList *pList){
541 // digitize module
409f8c84 542 const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
543 // for 3.6 eV/pair
544 const Float_t kconv = 10000.; // cm -> microns
545
99df08fd 546 Float_t spdLength = GetSeg()->Dz();
547 Float_t spdWidth = GetSeg()->Dx();
548 Float_t spdThickness = GetSeg()->Dy();
409f8c84 549 Float_t difCoef, dum;
99df08fd 550 GetResp()->DiffCoeff(difCoef,dum);
338e4f06 551 if(spdThickness > 290) difCoef = 0.00613;
409f8c84 552
553 Float_t zPix0 = 1e+6;
554 Float_t xPix0 = 1e+6;
555 Float_t yPrev = 1e+6;
556
99df08fd 557 Float_t zPitch = GetSeg()->Dpz(0);
558 Float_t xPitch = GetSeg()->Dpx(0);
409f8c84 559
560 TObjArray *fHits = mod->GetHits();
2cc6b29a 561 module = mod->GetIndex();
409f8c84 562 Int_t nhits = fHits->GetEntriesFast();
563 if (!nhits) return;
f74211b0 564#ifdef DEBUG
2cc6b29a 565 cout<<"len,wid,thickness,nx,nz,pitchx,pitchz,difcoef ="<<spdLength<<","
566 <<spdWidth<<","<<spdThickness<<","<<fNPixelsX<<","<<fNPixelsZ<<","
567 <<xPitch<<","<<zPitch<<","<<difCoef<<endl;
f74211b0 568#endif
2cc6b29a 569 // Array of pointers to the label-signal list
409f8c84 570 Int_t indexRange[4] = {0,0,0,0};
571
572 // Fill detector maps with GEANT hits
573 // loop over hits in the module
574 static Bool_t first;
575 Int_t lasttrack=-2;
576 Int_t hit, iZi, jz, jx;
338e4f06 577 Int_t idhit=-1; //!
f74211b0 578#ifdef DEBUG
338e4f06 579 cout<<"SPDdubna: module,nhits ="<<module<<","<<nhits<<endl;
f74211b0 580#endif
409f8c84 581 for (hit=0;hit<nhits;hit++) {
582 AliITShit *iHit = (AliITShit*) fHits->At(hit);
f74211b0 583#ifdef DEBUG
584 cout << "Hits=" << hit << "," << *iHit << endl;
585#endif
338e4f06 586 //Int_t layer = iHit->GetLayer();
587 Float_t yPix0 = -spdThickness/2;
588
589 // work with the idtrack=entry number in the TreeH
590 //Int_t idhit,idtrack; //!
591 //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit); //!
592 //Int_t idtrack=mod->GetHitTrackIndex(hit);
593 // or store straight away the particle position in the array
594 // of particles :
409f8c84 595 if(iHit->StatusEntering()) idhit=hit;
596 Int_t itrack = iHit->GetTrack();
597 Int_t dray = 0;
598
599 if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
600
338e4f06 601 //Int_t parent = iHit->GetParticle()->GetFirstMother();
409f8c84 602 Int_t partcode = iHit->GetParticle()->GetPdgCode();
603
2cc6b29a 604 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0,
605 // 211 - pi+, 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
409f8c84 606
607 Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
608 // vertex
609 pmod *= 1000;
610
409f8c84 611 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
612 // at p < 6 MeV/c
613
409f8c84 614 // Get hit z and x(r*phi) cordinates for each module (detector)
615 // in local system.
616
617 Float_t zPix = kconv*iHit->GetZL();
618 Float_t xPix = kconv*iHit->GetXL();
619 Float_t yPix = kconv*iHit->GetYL();
620
621 // Get track status
622 Int_t status = iHit->GetTrackStatus();
409f8c84 623
624 // Check boundaries
625 if(zPix > spdLength/2) {
f74211b0 626#ifdef DEBUG
627 cout<<"!!! SPD: z outside ="<<zPix<<endl;
628#endif
2cc6b29a 629 zPix = spdLength/2 - 10;
409f8c84 630 }
631 if(zPix < 0 && zPix < -spdLength/2) {
f74211b0 632#ifdef DEBUG
633 cout<<"!!! SPD: z outside ="<<zPix<<endl;
634#endif
2cc6b29a 635 zPix = -spdLength/2 + 10;
409f8c84 636 }
637 if(xPix > spdWidth/2) {
f74211b0 638#ifdef DEBUG
639 cout<<"!!! SPD: x outside ="<<xPix<<endl;
640#endif
2cc6b29a 641 xPix = spdWidth/2 - 10;
409f8c84 642 }
643 if(xPix < 0 && xPix < -spdWidth/2) {
f74211b0 644#ifdef DEBUG
645 cout<<"!!! SPD: x outside ="<<xPix<<endl;
646#endif
2cc6b29a 647 xPix = -spdWidth/2 + 10;
409f8c84 648 }
649 Int_t trdown = 0;
650
651 // enter Si or after event in Si
652 if (status == 66 ) {
2cc6b29a 653 zPix0 = zPix;
654 xPix0 = xPix;
655 yPrev = yPix;
656 } // end if status == 66
409f8c84 657
658 Float_t depEnergy = iHit->GetIonization();
659 // skip if the input point to Si
660
661 if(depEnergy <= 0.) continue;
662
663 // if track returns to the opposite direction:
664 if (yPix < yPrev) {
665 trdown = 1;
2cc6b29a 666 } // end if yPix < yPrev
409f8c84 667
668 // take into account the holes diffusion inside the Silicon
669 // the straight line between the entrance and exit points in Si is
670 // divided into the several steps; the diffusion is considered
671 // for each end point of step and charge
672 // is distributed between the pixels through the diffusion.
409f8c84 673
674 // ---------- the diffusion in Z (beam) direction -------
409f8c84 675 Float_t charge = depEnergy*kEnToEl; // charge in e-
676 Float_t drPath = 0.;
677 Float_t tang = 0.;
678 Float_t sigmaDif = 0.;
679 Float_t zdif = zPix - zPix0;
680 Float_t xdif = xPix - xPix0;
681 Float_t ydif = TMath::Abs(yPix - yPrev);
682 Float_t ydif0 = TMath::Abs(yPrev - yPix0);
683
684 if(ydif < 1) continue; // ydif is not zero
685
686 Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
687
688 Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
689 Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
690
691 // number of the steps along the track:
692 Int_t nsteps = ndZ;
693 if(ndX > ndZ) nsteps = ndX;
338e4f06 694 if(nsteps < 20) nsteps = 20; // minimum number of the steps
409f8c84 695
696 if (projDif < 5 ) {
2cc6b29a 697 drPath = (yPix-yPix0)*1.e-4;
698 drPath = TMath::Abs(drPath); // drift path in cm
699 sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
700 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
701 nsteps = 1;
702 } // end if projDif < 5
409f8c84 703
704 if(projDif > 5) tang = ydif/projDif;
705 Float_t dCharge = charge/nsteps; // charge in e- for one step
706 Float_t dZ = zdif/nsteps;
707 Float_t dX = xdif/nsteps;
708
2cc6b29a 709 for (iZi = 1; iZi <= nsteps;iZi++) {
710 Float_t dZn = iZi*dZ;
409f8c84 711 Float_t dXn = iZi*dX;
712 Float_t zPixn = zPix0 + dZn;
713 Float_t xPixn = xPix0 + dXn;
714
715 if(projDif >= 5) {
2cc6b29a 716 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
717 drPath = dProjn*tang*1.e-4; // drift path for iZi+1 step in cm
718 if(trdown == 0) {
719 drPath = TMath::Abs(drPath) + ydif0*1.e-4;
720 }// end if trdow ==0
721 if(trdown == 1) {
722 drPath = ydif0*1.e-4 - TMath::Abs(drPath);
723 drPath = TMath::Abs(drPath);
724 } // end if trdown == 1
725 sigmaDif = difCoef*sqrt(drPath);
726 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
727 } // end if projdif >= 5
409f8c84 728
729 zPixn = (zPixn + spdLength/2.);
730 xPixn = (xPixn + spdWidth/2.);
731 Int_t nZpix, nXpix;
99df08fd 732 GetSeg()->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
733 zPitch = GetSeg()->Dpz(nZpix);
734 GetSeg()->GetPadTxz(xPixn,zPixn);
409f8c84 735 // set the window for the integration
736 Int_t jzmin = 1;
737 Int_t jzmax = 3;
738 if(nZpix == 1) jzmin =2;
739 if(nZpix == fNPixelsZ) jzmax = 2;
740
741 Int_t jxmin = 1;
742 Int_t jxmax = 3;
743 if(nXpix == 1) jxmin =2;
744 if(nXpix == fNPixelsX) jxmax = 2;
745
746 Float_t zpix = nZpix;
747 Float_t dZright = zPitch*(zpix - zPixn);
748 Float_t dZleft = zPitch - dZright;
749
750 Float_t xpix = nXpix;
751 Float_t dXright = xPitch*(xpix - xPixn);
752 Float_t dXleft = xPitch - dXright;
753
754 Float_t dZprev = 0.;
755 Float_t dZnext = 0.;
756 Float_t dXprev = 0.;
757 Float_t dXnext = 0.;
758
759 for(jz=jzmin; jz <=jzmax; jz++) {
760 if(jz == 1) {
2cc6b29a 761 dZprev = -zPitch - dZleft;
762 dZnext = -dZleft;
763 } else if(jz == 2) {
764 dZprev = -dZleft;
765 dZnext = dZright;
766 } else if(jz == 3) {
767 dZprev = dZright;
768 dZnext = dZright + zPitch;
769 } // end if jz
409f8c84 770 // kz changes from 1 to the fNofPixels(270)
771 Int_t kz = nZpix + jz -2;
772
773 Float_t zArg1 = dZprev/sigmaDif;
774 Float_t zArg2 = dZnext/sigmaDif;
775 Float_t zProb1 = TMath::Erfc(zArg1);
776 Float_t zProb2 = TMath::Erfc(zArg2);
777 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
778
779
780 // ----------- holes diffusion in X(r*phi) direction --------
781
782 if(dZCharge > 1.) {
2cc6b29a 783 for(jx=jxmin; jx <=jxmax; jx++) {
784 if(jx == 1) {
785 dXprev = -xPitch - dXleft;
786 dXnext = -dXleft;
787 } else if(jx == 2) {
788 dXprev = -dXleft;
789 dXnext = dXright;
790 } else if(jx == 3) {
791 dXprev = dXright;
792 dXnext = dXright + xPitch;
793 } // end if jx
794 Int_t kx = nXpix + jx -2;
795 Float_t xArg1 = dXprev/sigmaDif;
796 Float_t xArg2 = dXnext/sigmaDif;
797 Float_t xProb1 = TMath::Erfc(xArg1);
798 Float_t xProb2 = TMath::Erfc(xArg2);
799 Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
800
801 if(dXCharge > 1.) {
2cc6b29a 802 if (first) {
f74211b0 803 indexRange[0]=indexRange[1]=kz-1;
2cc6b29a 804 indexRange[2]=indexRange[3]=kx-1;
805 first=kFALSE;
806 } // end if first
807 indexRange[0]=TMath::Min(indexRange[0],kz-1);
808 indexRange[1]=TMath::Max(indexRange[1],kz-1);
809 indexRange[2]=TMath::Min(indexRange[2],kx-1);
810 indexRange[3]=TMath::Max(indexRange[3],kx-1);
f74211b0 811/*
2cc6b29a 812 // build the list of digits for this module
f74211b0 813 Double_t signal = fMapA2->GetSignal(kz-1,kx-1);
2cc6b29a 814 signal+=dXCharge;
f74211b0 815 fMapA2->SetHit(kz-1,kx-1,(double)signal);
816*/
2cc6b29a 817 // The calling sequence for UpdateMapSignal was
818 // moved into the (dx > 1 e-) loop because it
819 // needs to call signal which is defined inside
820 // this loop
2cc6b29a 821 fModule = module;//Defined because functions
822 // called by UpdateMapSignal
823 // expect module to be an
824 // integer
7d50ea21 825 UpdateMapSignal(kz-1,kx-1,
ee86d557 826// mod->GetHitTrackIndex(hit),
827 ((AliITShit*)(mod->GetHit(hit)))->GetTrack(),
f74211b0 828 hit,fModule,dXCharge,pList);
2cc6b29a 829 } // dXCharge > 1 e-
830 } // jx loop
409f8c84 831 } // dZCharge > 1 e-
832 } // jz loop
833 } // iZi loop
409f8c84 834 if (status == 65) { // the step is inside of Si
2cc6b29a 835 zPix0 = zPix;
836 xPix0 = xPix;
837 } // end if status == 65
838 yPrev = yPix;
409f8c84 839 } // hit loop inside the module
409f8c84 840}
2cc6b29a 841//______________________________________________________________________
842void AliITSsimulationSPDdubna::ChargeToSignal(AliITSpList *pList){
843 // add noise and electronics, perform the zero suppression and add the
844 // digit to the list
f74211b0 845 static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
99df08fd 846 Float_t threshold = (float)GetResp()->MinVal();
f74211b0 847 Int_t j;
848// Int_t digits[3], tracks[3], hits[3];
849// Float_t charges[3];
2cc6b29a 850 Float_t electronics;
f74211b0 851// Float_t phys;
2cc6b29a 852 Double_t sig;
ee86d557 853 const Int_t nmaxtrk=AliITSdigitSPD::GetNTracks();
f74211b0 854 static AliITSdigitSPD dig;
855
2cc6b29a 856 for(Int_t iz=0; iz<fNPixelsZ; iz++){
857 for(Int_t ix=0; ix<fNPixelsX; ix++){
858 electronics = fBaseline + fNoise*gRandom->Gaus();
99df08fd 859 sig = pList->GetSignalOnly(iz+1,ix+1);
f74211b0 860 UpdateMapNoise(iz,ix,fModule,sig,electronics,pList);
861#ifdef DEBUG
862// cout << sig << "+" << electronics <<">threshold=" << threshold
863// << endl;
864#endif
865 if (sig+electronics > threshold) {
866 dig.fCoord1 = iz;
867 dig.fCoord2 = ix;
868 dig.fSignal = 1;
99df08fd 869 dig.fSignalSPD = (Int_t) pList->GetSignal(iz+1,ix+1);
f74211b0 870 /*
871 digits[0] = iz;
872 digits[1] = ix;
873 digits[2] = 1; */
874 for(j=0;j<nmaxtrk;j++){
875// charges[j] = 0.0;
ee86d557 876 if (j<pList->GetNEnteries()) {
99df08fd 877 dig.fTracks[j] = pList->GetTrack(iz+1,ix+1,j);
878 dig.fHits[j] = pList->GetHit(iz+1,ix+1,j);
f74211b0 879 /*
99df08fd 880 tracks[j] = pList->GetTrack(iz+1,ix+1,j);
881 hits[j] = pList->GetHit(iz+1,ix+1,j);
f74211b0 882 */
883 }else { // Default values
ee86d557 884 dig.fTracks[j] = -3;
885 dig.fHits[j] = -1;
f74211b0 886/* tracks[j] = -2; //noise
887 hits[j] = -1; */
2cc6b29a 888 } // end if pList
f74211b0 889 } // end for j
99df08fd 890// charges[0] = (Float_t) pList->GetSumSignal(iz+1,ix+1);
f74211b0 891/*
2cc6b29a 892 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
893 tracks[1] = -3;
894 hits[1] = -1;
895 tracks[2] = -3;
896 hits[2] = -1;
897 } else if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
898 tracks[1] = -3;
899 hits[1] = -1;
900 } else if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
901 tracks[2] = -3;
902 hits[2] = -1;
903 } else if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
904 tracks[2] = -3;
905 hits[2] = -1;
906 } // end if
f74211b0 907*/
908// phys = 0.0;
909#ifdef DEBUG
910 cout << iz << "," << ix << "," <<
99df08fd 911 *(pList->GetpListItem(iz+1,ix+1)) << endl;
f74211b0 912#endif
913// aliITS->AddSimDigit(0, phys, digits, tracks, hits, charges);
914 aliITS->AddSimDigit(0,&dig);
2cc6b29a 915 } //
916 } //
917 } //
409f8c84 918}
2cc6b29a 919//______________________________________________________________________
920void AliITSsimulationSPDdubna::CreateHistograms(){
921 // create 1D histograms for tests
922
923 printf("SPD - create histograms\n");
924
925 fHis=new TObjArray(fNPixelsZ);
926 TString spdName("spd_");
927 for (Int_t i=0;i<fNPixelsZ;i++) {
928 Char_t pixelz[4];
929 sprintf(pixelz,"%d",i+1);
930 spdName.Append(pixelz);
931 //PH (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
932 //PH fNPixelsX,0.,(Float_t) fNPixelsX);
933 fHis->AddAt(new TH1F(spdName.Data(),"SPD maps",
934 fNPixelsX,0.,(Float_t) fNPixelsX), i);
935 } // end for i
409f8c84 936}
2cc6b29a 937//______________________________________________________________________
938void AliITSsimulationSPDdubna::ResetHistograms(){
409f8c84 939 //
940 // Reset histograms for this detector
941 //
942
943 for ( int i=0;i<fNPixelsZ;i++ ) {
2cc6b29a 944 //PH if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
2682e810 945 if (fHis->At(i)) ((TH1F*)fHis->At(i))->Reset();
2cc6b29a 946 } // end for i
409f8c84 947}