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