ITS/AliITSsimulationSPD.cxx

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