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