| | 1 | /* |
| | 2 | questions to experts: why RemoveDeadPixels should be called before FrompListToDigits ? |
| | 3 | |
| | 4 | |
| | 5 | */ |
| | 6 | |
| | 7 | /************************************************************************** |
| | 8 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
| | 9 | * * |
| | 10 | * Author: The ALICE Off-line Project. * |
| | 11 | * Contributors are mentioned in the code where appropriate. * |
| | 12 | * * |
| | 13 | * Permission to use, copy, modify and distribute this software and its * |
| | 14 | * documentation strictly for non-commercial purposes is hereby granted * |
| | 15 | * without fee, provided that the above copyright notice appears in all * |
| | 16 | * copies and that both the copyright notice and this permission notice * |
| | 17 | * appear in the supporting documentation. The authors make no claims * |
| | 18 | * about the suitability of this software for any purpose. It is * |
| | 19 | * provided "as is" without express or implied warranty. * |
| | 20 | **************************************************************************/ |
| | 21 | |
| | 22 | #include <TGeoGlobalMagField.h> |
| | 23 | #include <TH1.h> |
| | 24 | #include <TString.h> |
| | 25 | #include "AliITSU.h" |
| | 26 | #include "AliITSUDigitPix.h" |
| | 27 | #include "AliITSUHit.h" |
| | 28 | #include "AliITSUModule.h" |
| | 29 | #include "AliITSUSensMap.h" |
| | 30 | #include "AliITSUCalibrationPix.h" |
| | 31 | #include "AliITSUSegmentationPix.h" |
| | 32 | #include "AliITSUSimulationPix.h" |
| | 33 | #include "AliLog.h" |
| | 34 | #include "AliRun.h" |
| | 35 | #include "AliMagF.h" |
| | 36 | #include "AliMathBase.h" |
| | 37 | #include "AliITSUSimuParam.h" |
| | 38 | #include "AliITSUSDigit.h" |
| | 39 | |
| | 40 | ClassImp(AliITSUSimulationPix) |
| | 41 | //////////////////////////////////////////////////////////////////////// |
| | 42 | // Version: 1 |
| | 43 | // Modified by D. Elia, G.E. Bruno, H. Tydesjo |
| | 44 | // Fast diffusion code by Bjorn S. Nilsen |
| | 45 | // March-April 2006 |
| | 46 | // October 2007: GetCalibrationObjects() removed |
| | 47 | // |
| | 48 | // Version: 0 |
| | 49 | // Written by Boris Batyunya |
| | 50 | // December 20 1999 |
| | 51 | // |
| | 52 | // Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch |
| | 53 | // |
| | 54 | // AliITSUSimulationPix is to do the simulation of pixels |
| | 55 | // |
| | 56 | //////////////////////////////////////////////////////////////////////// |
| | 57 | |
| | 58 | //______________________________________________________________________ |
| | 59 | AliITSUSimulationPix::AliITSUSimulationPix() |
| | 60 | : fTanLorAng(0) |
| | 61 | ,fStrobe(kTRUE) |
| | 62 | ,fStrobeLenght(4) |
| | 63 | ,fStrobePhase(-12.5e-9) |
| | 64 | { |
| | 65 | // Default constructor. |
| | 66 | } |
| | 67 | |
| | 68 | //______________________________________________________________________ |
| | 69 | AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map) |
| | 70 | :AliITSUSimulation(sim,map) |
| | 71 | ,fTanLorAng(0) |
| | 72 | ,fStrobe(kTRUE) |
| | 73 | ,fStrobeLenght(4) |
| | 74 | ,fStrobePhase(-12.5e-9) |
| | 75 | { |
| | 76 | // standard constructor |
| | 77 | Init(); |
| | 78 | } |
| | 79 | |
| | 80 | //______________________________________________________________________ |
| | 81 | AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s) |
| | 82 | :AliITSUSimulation(s) |
| | 83 | ,fTanLorAng(s.fTanLorAng) |
| | 84 | ,fStrobe(s.fStrobe) |
| | 85 | ,fStrobeLenght(s.fStrobeLenght) |
| | 86 | ,fStrobePhase(s.fStrobePhase) |
| | 87 | { |
| | 88 | // Copy Constructor |
| | 89 | } |
| | 90 | |
| | 91 | |
| | 92 | //______________________________________________________________________ |
| | 93 | AliITSUSimulationPix::~AliITSUSimulationPix() |
| | 94 | { |
| | 95 | // destructor |
| | 96 | // only the sens map is owned and it is deleted by ~AliITSUSimulation |
| | 97 | } |
| | 98 | |
| | 99 | //______________________________________________________________________ |
| | 100 | AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s) |
| | 101 | { |
| | 102 | // Assignment operator |
| | 103 | if (&s == this) return *this; |
| | 104 | AliITSUSimulation::operator=(s); |
| | 105 | fStrobe = s.fStrobe; |
| | 106 | fStrobeLenght = s.fStrobeLenght; |
| | 107 | fStrobePhase = s.fStrobePhase; |
| | 108 | return *this; |
| | 109 | } |
| | 110 | |
| | 111 | //______________________________________________________________________ |
| | 112 | void AliITSUSimulationPix::Init() |
| | 113 | { |
| | 114 | // Initilization |
| | 115 | if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight()); |
| | 116 | } |
| | 117 | |
| | 118 | //______________________________________________________________________ |
| | 119 | Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole) |
| | 120 | { |
| | 121 | // This function set the Tangent of the Lorentz angle. |
| | 122 | // A weighted average is used for electrons and holes |
| | 123 | // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1] |
| | 124 | // output: Bool_t : kTRUE in case of success |
| | 125 | // |
| | 126 | if (weightHole<0) { |
| | 127 | weightHole=0.; |
| | 128 | AliWarning("You have asked for negative Hole weight"); |
| | 129 | AliWarning("I'm going to use only electrons"); |
| | 130 | } |
| | 131 | if (weightHole>1) { |
| | 132 | weightHole=1.; |
| | 133 | AliWarning("You have asked for weight > 1"); |
| | 134 | AliWarning("I'm going to use only holes"); |
| | 135 | } |
| | 136 | Double_t weightEle=1.-weightHole; |
| | 137 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); |
| | 138 | if (!fld) AliFatal("The field is not initialized"); |
| | 139 | Double_t bz = fld->SolenoidField(); |
| | 140 | fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) + |
| | 141 | weightEle*fSimuParam->LorentzAngleElectron(bz)); |
| | 142 | return kTRUE; |
| | 143 | } |
| | 144 | |
| | 145 | //_____________________________________________________________________ |
| | 146 | void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg) |
| | 147 | { |
| | 148 | // This function begins the work of creating S-Digits. |
| | 149 | if (!(mod->GetNHits())) { |
| | 150 | AliDebug(1,Form("In event %d module %d there are %d hits returning.", |
| | 151 | event, mod->GetIndex(),mod->GetNHits())); |
| | 152 | return; |
| | 153 | } |
| | 154 | // |
| | 155 | if (fStrobe) if (event != fEvent) GenerateStrobePhase(); |
| | 156 | InitSimulationModule(mod->GetIndex(), event, seg ); |
| | 157 | // Hits2SDigits(mod); |
| | 158 | Hits2SDigitsFast(mod); |
| | 159 | if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels(); |
| | 160 | if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels(); |
| | 161 | WriteSDigits(); |
| | 162 | ClearMap(); |
| | 163 | } |
| | 164 | |
| | 165 | //______________________________________________________________________ |
| | 166 | void AliITSUSimulationPix::WriteSDigits() |
| | 167 | { |
| | 168 | // This function adds each S-Digit to pList |
| | 169 | static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS"); |
| | 170 | int nsd = fSensMap->GetEntries(); |
| | 171 | for (int i=0;i<nsd;i++) { |
| | 172 | AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index |
| | 173 | if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue; |
| | 174 | aliITS->AddSumDigit(*sd); |
| | 175 | } |
| | 176 | return; |
| | 177 | } |
| | 178 | |
| | 179 | //______________________________________________________________________ |
| | 180 | void AliITSUSimulationPix::FinishSDigitiseModule() |
| | 181 | { |
| | 182 | // This function calls SDigitsToDigits which creates Digits from SDigits |
| | 183 | FrompListToDigits(); |
| | 184 | ClearMap(); |
| | 185 | return; |
| | 186 | } |
| | 187 | |
| | 188 | //______________________________________________________________________ |
| | 189 | void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg) |
| | 190 | { |
| | 191 | // This function creates Digits straight from the hits and then adds |
| | 192 | // electronic noise to the digits before adding them to pList |
| | 193 | // Each of the input variables is passed along to Hits2SDigits |
| | 194 | // |
| | 195 | if (fStrobe) if (event != fEvent) GenerateStrobePhase(); |
| | 196 | InitSimulationModule( mod->GetIndex(), event, seg ); |
| | 197 | // Hits2SDigits(mod); |
| | 198 | Hits2SDigitsFast(mod); |
| | 199 | // |
| | 200 | if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels(); |
| | 201 | if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels(); |
| | 202 | FrompListToDigits(); |
| | 203 | ClearMap(); |
| | 204 | } |
| | 205 | |
| | 206 | //______________________________________________________________________ |
| | 207 | void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod) |
| | 208 | { |
| | 209 | // Does the charge distributions using Gaussian diffusion charge charing. |
| | 210 | const Double_t kmictocm = 1.0e-4; // convert microns to cm. |
| | 211 | const Double_t kcmtomic = 1.e4; |
| | 212 | const Double_t kBunchLenght = 25e-9; // LHC clock |
| | 213 | Int_t nhits = mod->GetNHits(); |
| | 214 | if (!nhits) return; |
| | 215 | // |
| | 216 | Int_t h,ix,iz,i; |
| | 217 | Int_t idtrack; |
| | 218 | Float_t x,y,z; // keep coordinates float (required by AliSegmentation) |
| | 219 | 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; |
| | 220 | Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda; |
| | 221 | Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half Thickness |
| | 222 | fSimuParam->GetPixSigmaDiffusionAsymmetry(fda); |
| | 223 | // |
| | 224 | for (h=0;h<nhits;h++) { |
| | 225 | if (fStrobe && |
| | 226 | ((mod->GetHit(h)->GetTOF()<fStrobePhase) || |
| | 227 | (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght))) |
| | 228 | ) continue; |
| | 229 | // |
| | 230 | if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; |
| | 231 | st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); |
| | 232 | if (st>0.0) { |
| | 233 | st = (Double_t)((Int_t)(st*kcmtomic)); // number of microns |
| | 234 | if (st<=1.0) st = 1.0; |
| | 235 | dt = 1.0/st; |
| | 236 | for (t=0.0;t<1.0;t+=dt) { // Integrate over t |
| | 237 | tp = t+0.5*dt; |
| | 238 | x = x0+x1*tp; |
| | 239 | y = y0+y1*tp; |
| | 240 | z = z0+z1*tp; |
| | 241 | if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside |
| | 242 | el = dt * de / fSimuParam->GetGeVToCharge(); |
| | 243 | // |
| | 244 | sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); |
| | 245 | sigx=sig; |
| | 246 | sigz=sig*fda; |
| | 247 | if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng; |
| | 248 | SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); |
| | 249 | } // end for t |
| | 250 | } else { // st == 0.0 deposit it at this point |
| | 251 | x = x0; |
| | 252 | y = y0; |
| | 253 | z = z0; |
| | 254 | if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside |
| | 255 | el = de / fSimuParam->GetGeVToCharge(); |
| | 256 | sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); |
| | 257 | sigx = sig; |
| | 258 | sigz = sig*fda; |
| | 259 | if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng; |
| | 260 | SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); |
| | 261 | } // end if st>0.0 |
| | 262 | } // Loop over all hits h |
| | 263 | // |
| | 264 | // Coupling |
| | 265 | int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster |
| | 266 | AliITSUSDigit* dg = 0; |
| | 267 | switch (fSimuParam->GetPixCouplingOption()) { |
| | 268 | case AliITSUSimuParam::kNewCouplingPix : |
| | 269 | for (i=nd;i--;) { |
| | 270 | dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); |
| | 271 | if (fSensMap->IsDisabled(dg)) continue; |
| | 272 | SetCoupling(dg,idtrack,h); |
| | 273 | } |
| | 274 | break; |
| | 275 | case AliITSUSimuParam::kOldCouplingPix: |
| | 276 | for (i=nd;i--;) { |
| | 277 | dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); |
| | 278 | if (fSensMap->IsDisabled(dg)) continue; |
| | 279 | SetCouplingOld(dg,idtrack,h); |
| | 280 | } |
| | 281 | break; |
| | 282 | default: |
| | 283 | break; |
| | 284 | |
| | 285 | } // end switch |
| | 286 | if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption())); |
| | 287 | } |
| | 288 | |
| | 289 | //______________________________________________________________________ |
| | 290 | void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod) |
| | 291 | { |
| | 292 | // Does the charge distributions using Gaussian diffusion charge charing. // Inputs: |
| | 293 | // AliITSUModule *mod Pointer to this module |
| | 294 | // Output: |
| | 295 | // none. |
| | 296 | // Return: |
| | 297 | // none. |
| | 298 | const Double_t kmictocm = 1.0e-4; // convert microns to cm. |
| | 299 | const Int_t kn10=10; |
| | 300 | const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1, |
| | 301 | 4.325316833e-1,4.869532643e-1,5.130467358e-1, |
| | 302 | 5.674683167e-1,6.602952159e-1,7.833023029e-1, |
| | 303 | 9.255628306e-1}; |
| | 304 | const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1, |
| | 305 | 7.472567455e-2,3.333567215e-2,3.333567215e-2, |
| | 306 | 7.472567455e-2,1.095431813e-1,1.346333597e-1, |
| | 307 | 1.477621124e-1}; |
| | 308 | const Double_t kBunchLenght = 25e-9; // LHC clock |
| | 309 | TObjArray *hits = mod->GetHits(); |
| | 310 | Int_t nhits = hits->GetEntriesFast(); |
| | 311 | if (nhits<=0) return; |
| | 312 | // |
| | 313 | Int_t h,ix,iz,i; |
| | 314 | Int_t idtrack; |
| | 315 | Float_t x,y,z; // keep coordinates float (required by AliSegmentation) |
| | 316 | Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; |
| | 317 | Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0; |
| | 318 | Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half thickness |
| | 319 | fSimuParam->GetPixSigmaDiffusionAsymmetry(fda); |
| | 320 | // |
| | 321 | for (h=0;h<nhits;h++) { |
| | 322 | // Check if the hit is inside readout window |
| | 323 | if (fStrobe && |
| | 324 | ((mod->GetHit(h)->GetTOF()<fStrobePhase) || |
| | 325 | (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght))) |
| | 326 | ) continue; |
| | 327 | // |
| | 328 | if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue; |
| | 329 | st = TMath::Sqrt(x1*x1+y1*y1+z1*z1); |
| | 330 | if (st>0.0) |
| | 331 | for (i=0;i<kn10;i++) { // Integrate over t |
| | 332 | t = kti[i]; |
| | 333 | x = x0+x1*t; |
| | 334 | y = y0+y1*t; |
| | 335 | z = z0+z1*t; |
| | 336 | if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside |
| | 337 | el = kwi[i]*de/fSimuParam->GetGeVToCharge(); |
| | 338 | sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); |
| | 339 | sigx=sig; |
| | 340 | sigz=sig*fda; |
| | 341 | if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng; |
| | 342 | SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); |
| | 343 | } // end for i // End Integrate over t |
| | 344 | else { // st == 0.0 deposit it at this point |
| | 345 | x = x0; |
| | 346 | y = y0; |
| | 347 | z = z0; |
| | 348 | if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside |
| | 349 | el = de / fSimuParam->GetGeVToCharge(); |
| | 350 | sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); |
| | 351 | sigx=sig; |
| | 352 | sigz=sig*fda; |
| | 353 | if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng; |
| | 354 | SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h); |
| | 355 | } // end if st>0.0 |
| | 356 | |
| | 357 | } // Loop over all hits h |
| | 358 | |
| | 359 | // Coupling |
| | 360 | int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster |
| | 361 | AliITSUSDigit* dg = 0; |
| | 362 | switch (fSimuParam->GetPixCouplingOption()) { |
| | 363 | case AliITSUSimuParam::kNewCouplingPix : |
| | 364 | for (i=nd;i--;) { |
| | 365 | dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); |
| | 366 | if (fSensMap->IsDisabled(dg)) continue; |
| | 367 | SetCoupling(dg,idtrack,h); |
| | 368 | } |
| | 369 | case AliITSUSimuParam::kOldCouplingPix: |
| | 370 | for (i=nd;i--;) { |
| | 371 | dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i); |
| | 372 | if (fSensMap->IsDisabled(dg)) continue; |
| | 373 | SetCouplingOld(dg,idtrack,h); |
| | 374 | } |
| | 375 | break; |
| | 376 | default: |
| | 377 | break; |
| | 378 | } // end switch |
| | 379 | if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption())); |
| | 380 | } |
| | 381 | |
| | 382 | //______________________________________________________________________ |
| | 383 | void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0, |
| | 384 | Int_t ix0,Int_t iz0, |
| | 385 | Double_t el,Double_t sig,Double_t ld, |
| | 386 | Int_t t,Int_t hi) |
| | 387 | { |
| | 388 | // Spreads the charge over neighboring cells. Assume charge is distributed |
| | 389 | // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg) |
| | 390 | // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig) |
| | 391 | // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account) |
| | 392 | // Defined this way, the integral over all x and z is el. |
| | 393 | // Inputs: |
| | 394 | // Double_t x0 x position of point where charge is liberated |
| | 395 | // Double_t z0 z position of point where charge is liberated |
| | 396 | // Int_t ix0 row of cell corresponding to point x0 |
| | 397 | // Int_t iz0 columb of cell corresponding to point z0 |
| | 398 | // Double_t el number of electrons liberated in this step |
| | 399 | // Double_t sig Sigma difusion for this step (y0 dependent) |
| | 400 | // Double_t ld lorentz drift in x for this step (y0 dependent) |
| | 401 | // Int_t t track number |
| | 402 | // Int_t ti hit track index number |
| | 403 | // Int_t hi hit "hit" index number |
| | 404 | // Outputs: |
| | 405 | // none. |
| | 406 | // Return: |
| | 407 | // none. |
| | 408 | const Int_t knx = 3,knz = 2; |
| | 409 | const Double_t kRoot2 = 1.414213562; // Sqrt(2). |
| | 410 | const Double_t kmictocm = 1.0e-4; // convert microns to cm. |
| | 411 | Int_t ix,iz,ixs,ixe,izs,ize; |
| | 412 | Float_t x,z; // keep coordinates float (required by AliSegmentation) |
| | 413 | Double_t s,sp,x1,x2,z1,z2; |
| | 414 | // |
| | 415 | if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi)); |
| | 416 | if (sig<=0.0) { // if sig<=0 No diffusion to simulate. |
| | 417 | UpdateMapSignal(iz0,ix0,t,hi,el); |
| | 418 | return; |
| | 419 | } // end if |
| | 420 | sp = 1.0/(sig*kRoot2); |
| | 421 | ixs = TMath::Max(-knx+ix0,0); |
| | 422 | ixe = TMath::Min(knx+ix0,fSeg->Npx()-1); |
| | 423 | izs = TMath::Max(-knz+iz0,0); |
| | 424 | ize = TMath::Min(knz+iz0,fSeg->Npz()-1); |
| | 425 | for (ix=ixs;ix<=ixe;ix++) |
| | 426 | for (iz=izs;iz<=ize;iz++) { |
| | 427 | fSeg->DetToLocal(ix,iz,x,z); // pixel center |
| | 428 | double dxi = 0.5*kmictocm*fSeg->Dpx(ix); |
| | 429 | double dzi = 0.5*kmictocm*fSeg->Dpz(iz); |
| | 430 | x1 = x; |
| | 431 | z1 = z; |
| | 432 | x2 = x1 + dxi; // Upper |
| | 433 | x1 -= dxi; // Lower |
| | 434 | z2 = z1 + dzi; // Upper |
| | 435 | z1 -= dzi; // Lower |
| | 436 | x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) |
| | 437 | x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) |
| | 438 | z1 -= z0; // Distance from where track traveled |
| | 439 | z2 -= z0; // Distance from where track traveled |
| | 440 | s = el*0.25; // Correction based on definision of Erfc |
| | 441 | s *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2); |
| | 442 | s *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2); |
| | 443 | if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s); |
| | 444 | } // end for ix, iz |
| | 445 | // |
| | 446 | } |
| | 447 | |
| | 448 | //______________________________________________________________________ |
| | 449 | void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0, |
| | 450 | Int_t ix0,Int_t iz0, |
| | 451 | Double_t el,Double_t sigx,Double_t sigz, |
| | 452 | Double_t ld,Int_t t,Int_t hi) |
| | 453 | { |
| | 454 | // Spreads the charge over neighboring cells. Assume charge is distributed |
| | 455 | // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg) |
| | 456 | // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz) |
| | 457 | // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account) |
| | 458 | // Defined this way, the integral over all x and z is el. |
| | 459 | // Inputs: |
| | 460 | // Double_t x0 x position of point where charge is liberated |
| | 461 | // Double_t z0 z position of point where charge is liberated |
| | 462 | // Int_t ix0 row of cell corresponding to point x0 |
| | 463 | // Int_t iz0 columb of cell corresponding to point z0 |
| | 464 | // Double_t el number of electrons liberated in this step |
| | 465 | // Double_t sigx Sigma difusion along x for this step (y0 dependent) |
| | 466 | // Double_t sigz Sigma difusion along z for this step (z0 dependent) |
| | 467 | // Double_t ld lorentz drift in x for this stip (y0 dependent) |
| | 468 | // Int_t t track number |
| | 469 | // Int_t ti hit track index number |
| | 470 | // Int_t hi hit "hit" index number |
| | 471 | // Outputs: |
| | 472 | // none. |
| | 473 | // Return: |
| | 474 | // none. |
| | 475 | const Int_t knx = 3,knz = 3; // RS: TO TUNE |
| | 476 | const Double_t kRoot2 = 1.414213562; // Sqrt(2). |
| | 477 | const Double_t kmictocm = 1.0e-4; // convert microns to cm. |
| | 478 | Int_t ix,iz,ixs,ixe,izs,ize; |
| | 479 | Float_t x,z; // keep coordinates float (required by AliSegmentation) |
| | 480 | Double_t s,spx,spz,x1,x2,z1,z2; |
| | 481 | // |
| | 482 | if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d,ld=%e)", |
| | 483 | x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld)); |
| | 484 | if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate. |
| | 485 | UpdateMapSignal(iz0,ix0,t,hi,el); |
| | 486 | return; |
| | 487 | } // end if |
| | 488 | spx = 1.0/(sigx*kRoot2); |
| | 489 | spz = 1.0/(sigz*kRoot2); |
| | 490 | ixs = TMath::Max(-knx+ix0,0); |
| | 491 | ixe = TMath::Min(knx+ix0,fSeg->Npx()-1); |
| | 492 | izs = TMath::Max(-knz+iz0,0); |
| | 493 | ize = TMath::Min(knz+iz0,fSeg->Npz()-1); |
| | 494 | for (ix=ixs;ix<=ixe;ix++) |
| | 495 | for (iz=izs;iz<=ize;iz++) { |
| | 496 | fSeg->DetToLocal(ix,iz,x,z); // pixel center |
| | 497 | double dxi = 0.5*kmictocm*fSeg->Dpx(ix); |
| | 498 | double dzi = 0.5*kmictocm*fSeg->Dpz(iz); |
| | 499 | x1 = x; |
| | 500 | z1 = z; |
| | 501 | x2 = x1 + dxi; // Upper |
| | 502 | x1 -= dxi; // Lower |
| | 503 | z2 = z1 + dzi; // Upper |
| | 504 | z1 -= dzi; // Lower |
| | 505 | x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) |
| | 506 | x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift) |
| | 507 | z1 -= z0; // Distance from where track traveled |
| | 508 | z2 -= z0; // Distance from where track traveled |
| | 509 | s = el*0.25; // Correction based on definision of Erfc |
| | 510 | s *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2); |
| | 511 | s *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2); |
| | 512 | if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s); |
| | 513 | } // end for ix, iz |
| | 514 | // |
| | 515 | } |
| | 516 | |
| | 517 | //______________________________________________________________________ |
| | 518 | void AliITSUSimulationPix::RemoveDeadPixels() |
| | 519 | { |
| | 520 | // Removes dead pixels on each module (ladder) |
| | 521 | // This should be called before going from sdigits to digits (FrompListToDigits) |
| | 522 | |
| | 523 | AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead(); |
| | 524 | if (!calObj) return; |
| | 525 | // |
| | 526 | if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked |
| | 527 | // |
| | 528 | // remove single bad pixels one by one |
| | 529 | int nsingle = calObj->GetNrBadSingle(); |
| | 530 | UInt_t col,row; |
| | 531 | for (int i=nsingle;i--;) { |
| | 532 | calObj->GetBadPixelSingle(i,row,col); |
| | 533 | fSensMap->DeleteItem(col,row); |
| | 534 | } |
| | 535 | int nsd = fSensMap->GetEntriesUnsorted(); |
| | 536 | for (int isd=nsd;isd--;) { |
| | 537 | AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd); |
| | 538 | if (fSensMap->IsDisabled(sd)) continue; |
| | 539 | fSensMap->GetMapIndex(sd->GetUniqueID(),col,row); |
| | 540 | int chip = fSeg->GetChipFromChannel(0,col); |
| | 541 | // if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad |
| | 542 | if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list |
| | 543 | } |
| | 544 | // |
| | 545 | } |
| | 546 | |
| | 547 | //______________________________________________________________________ |
| | 548 | void AliITSUSimulationPix::AddNoisyPixels() |
| | 549 | { |
| | 550 | // Adds noisy pixels on each module (ladder) |
| | 551 | // This should be called before going from sdigits to digits (FrompListToDigits) |
| | 552 | AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy(); |
| | 553 | if (!calObj) return; |
| | 554 | for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), |
| | 555 | 10*fSimuParam->GetPixThreshold(fModule)); |
| | 556 | // |
| | 557 | } |
| | 558 | |
| | 559 | //______________________________________________________________________ |
| | 560 | void AliITSUSimulationPix::FrompListToDigits() |
| | 561 | { |
| | 562 | // add noise and electronics, perform the zero suppression and add the |
| | 563 | // digit to the list |
| | 564 | static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS"); |
| | 565 | UInt_t ix,iz; |
| | 566 | Double_t sig; |
| | 567 | const Int_t knmaxtrk=AliITSdigit::GetNTracks(); |
| | 568 | static AliITSUDigitPix dig; |
| | 569 | // RS: in principle: |
| | 570 | // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold. |
| | 571 | // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold |
| | 572 | // With many channels this will be too time consuming, hence I do the following |
| | 573 | // 1) Precalculate the probability that the nois alone will exceed the threshold. |
| | 574 | // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold. |
| | 575 | // 3) For pixels having a hits apply the usual noise and compare with threshold |
| | 576 | // |
| | 577 | // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf |
| | 578 | // |
| | 579 | int maxInd = fSensMap->GetMaxIndex(); |
| | 580 | double minProb = 0.1/maxInd; |
| | 581 | // |
| | 582 | int nsd = fSensMap->GetEntries(); |
| | 583 | Int_t prevID=0,curID=0; |
| | 584 | TArrayI ordSampleInd(100),ordSample(100); |
| | 585 | // |
| | 586 | double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule); |
| | 587 | fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean); |
| | 588 | probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold |
| | 589 | // |
| | 590 | for (int i=0;i<nsd;i++) { |
| | 591 | AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index |
| | 592 | if (fSensMap->IsDisabled(sd)) continue; |
| | 593 | curID = (int)sd->GetUniqueID(); |
| | 594 | // |
| | 595 | if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current |
| | 596 | CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean); |
| | 597 | prevID = curID+1; |
| | 598 | } |
| | 599 | // |
| | 600 | if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue; |
| | 601 | if (TMath::Abs(sig)>2147483647.0) { //RS? |
| | 602 | //PH 2147483647 is the max. integer |
| | 603 | //PH This apparently is a problem which needs investigation |
| | 604 | AliWarning(Form("Too big or too small signal value %f",sig)); |
| | 605 | } |
| | 606 | fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix); |
| | 607 | dig.SetCoord1(iz); |
| | 608 | dig.SetCoord2(ix); |
| | 609 | dig.SetSignal(1); |
| | 610 | dig.SetSignalPix((Int_t)sig); |
| | 611 | int ntr = sd->GetNTracks(); |
| | 612 | for (int j=0;j<ntr;j++) { |
| | 613 | dig.SetTrack(j,sd->GetTrack(j)); |
| | 614 | dig.SetHit(j,sd->GetHit(j)); |
| | 615 | } |
| | 616 | for (int j=ntr;j<knmaxtrk;j++) { |
| | 617 | dig.SetTrack(j,-3); |
| | 618 | dig.SetHit(j,-1); |
| | 619 | } |
| | 620 | aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig); |
| | 621 | } |
| | 622 | // if needed, add noisy pixels with id from last real hit to maxID |
| | 623 | if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean); |
| | 624 | // |
| | 625 | } |
| | 626 | |
| | 627 | //______________________________________________________________________ |
| | 628 | Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base) |
| | 629 | { |
| | 630 | // create random noisy digits above threshold within id range [minID,maxID[ |
| | 631 | // see FrompListToDigits for details |
| | 632 | // |
| | 633 | static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS"); |
| | 634 | UInt_t ix,iz; |
| | 635 | static AliITSUDigitPix dig; |
| | 636 | static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!! |
| | 637 | const Int_t knmaxtrk=AliITSdigit::GetNTracks(); |
| | 638 | // |
| | 639 | Int_t ncand = 0; |
| | 640 | int npix = maxID-minID; |
| | 641 | if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added |
| | 642 | ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd); |
| | 643 | int* ordV = ordSample.GetArray(); |
| | 644 | int* ordI = ordSampleInd.GetArray(); |
| | 645 | for (int j=0;j<ncand;j++) { |
| | 646 | fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix); // create noisy digit |
| | 647 | dig.SetCoord1(iz); |
| | 648 | dig.SetCoord2(ix); |
| | 649 | dig.SetSignal(1); |
| | 650 | dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise)); |
| | 651 | for (int k=knmaxtrk;k--;) { |
| | 652 | dig.SetTrack(k,-3); |
| | 653 | dig.SetHit(k,-1); |
| | 654 | } |
| | 655 | aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig); |
| | 656 | if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix())); |
| | 657 | } |
| | 658 | return ncand; |
| | 659 | } |
| | 660 | |
| | 661 | //______________________________________________________________________ |
| | 662 | void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit) |
| | 663 | { |
| | 664 | // Take into account the coupling between adiacent pixels. |
| | 665 | // The parameters probcol and probrow are the probability of the |
| | 666 | // signal in one pixel shared in the two adjacent pixels along |
| | 667 | // the column and row direction, respectively. |
| | 668 | // Note pList is goten via GetMap() and module is not need any more. |
| | 669 | // Otherwise it is identical to that coded by Tiziano Virgili (BSN). |
| | 670 | //Begin_Html |
| | 671 | /* |
| | 672 | <img src="picts/ITS/barimodel_3.gif"> |
| | 673 | </pre> |
| | 674 | <br clear=left> |
| | 675 | <font size=+2 color=red> |
| | 676 | <a href="mailto:tiziano.virgili@cern.ch"></a>. |
| | 677 | </font> |
| | 678 | <pre> |
| | 679 | */ |
| | 680 | //End_Html |
| | 681 | // Inputs: |
| | 682 | // old existing AliITSUSDigit |
| | 683 | // Int_t ntrack track incex number |
| | 684 | // Int_t idhit hit index number |
| | 685 | UInt_t col,row; |
| | 686 | Double_t pulse1,pulse2; |
| | 687 | Double_t couplR=0.0,couplC=0.0; |
| | 688 | Double_t xr=0.; |
| | 689 | // |
| | 690 | fSensMap->GetMapIndex(old->GetUniqueID(),col,row); |
| | 691 | fSimuParam->GetPixCouplingParam(couplC,couplR); |
| | 692 | if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d) couplC=%e couplR=%e", |
| | 693 | col,row,ntrack,idhit,couplC,couplR)); |
| | 694 | pulse2 = pulse1 = old->GetSignal(); |
| | 695 | if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal |
| | 696 | for (Int_t isign=-1;isign<=1;isign+=2) { |
| | 697 | // |
| | 698 | // loop in col direction |
| | 699 | int j1 = int(col) + isign; |
| | 700 | xr = gRandom->Rndm(); |
| | 701 | if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1); |
| | 702 | // |
| | 703 | // loop in row direction |
| | 704 | int j2 = int(row) + isign; |
| | 705 | xr = gRandom->Rndm(); |
| | 706 | if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2); |
| | 707 | } |
| | 708 | // |
| | 709 | } |
| | 710 | |
| | 711 | //______________________________________________________________________ |
| | 712 | void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit) |
| | 713 | { |
| | 714 | // Take into account the coupling between adiacent pixels. |
| | 715 | // The parameters probcol and probrow are the fractions of the |
| | 716 | // signal in one pixel shared in the two adjacent pixels along |
| | 717 | // the column and row direction, respectively. |
| | 718 | //Begin_Html |
| | 719 | /* |
| | 720 | <img src="picts/ITS/barimodel_3.gif"> |
| | 721 | </pre> |
| | 722 | <br clear=left> |
| | 723 | <font size=+2 color=red> |
| | 724 | <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>. |
| | 725 | </font> |
| | 726 | <pre> |
| | 727 | */ |
| | 728 | //End_Html |
| | 729 | // Inputs: |
| | 730 | // old existing AliITSUSDigit |
| | 731 | // ntrack track incex number |
| | 732 | // idhit hit index number |
| | 733 | // module module number |
| | 734 | // |
| | 735 | UInt_t col,row; |
| | 736 | Double_t pulse1,pulse2; |
| | 737 | Double_t couplR=0.0,couplC=0.0; |
| | 738 | // |
| | 739 | fSensMap->GetMapIndex(old->GetUniqueID(),col,row); |
| | 740 | fSimuParam->GetPixCouplingParam(couplC,couplR); |
| | 741 | if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d) couplC=%e couplR=%e", |
| | 742 | col,row,ntrack,idhit,couplC,couplR)); |
| | 743 | // |
| | 744 | if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal |
| | 745 | for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction |
| | 746 | pulse2 = pulse1 = old->GetSignal(); |
| | 747 | // |
| | 748 | int j1 = int(col)+isign; |
| | 749 | pulse1 *= couplC; |
| | 750 | if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal(); |
| | 751 | else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1); |
| | 752 | |
| | 753 | // loop in row direction |
| | 754 | int j2 = int(row) + isign; |
| | 755 | pulse2 *= couplR; |
| | 756 | if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal(); |
| | 757 | else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2); |
| | 758 | } // for isign |
| | 759 | } |
| | 760 | |
| | 761 | //______________________________________________________________________ |
| | 762 | void AliITSUSimulationPix::GenerateStrobePhase() |
| | 763 | { |
| | 764 | // Generate randomly the strobe |
| | 765 | // phase w.r.t to the LHC clock |
| | 766 | // Done once per event |
| | 767 | const Double_t kBunchLenght = 25e-9; // LHC clock |
| | 768 | fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght- |
| | 769 | (Double_t)fStrobeLenght*kBunchLenght+ |
| | 770 | kBunchLenght/2; |
| | 771 | } |
| | 772 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental