Simulation response function treatment.
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 5 Feb 2013 20:15:10 +0000 (20:15 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 5 Feb 2013 20:15:10 +0000 (20:15 +0000)
13 files changed:
ITS/UPGRADE/AliITSU.cxx
ITS/UPGRADE/AliITSU.h
ITS/UPGRADE/AliITSUDigitizer.cxx
ITS/UPGRADE/AliITSURecoDet.cxx
ITS/UPGRADE/AliITSUSimuParam.cxx
ITS/UPGRADE/AliITSUSimuParam.h
ITS/UPGRADE/AliITSUSimulation.cxx
ITS/UPGRADE/AliITSUSimulation.h
ITS/UPGRADE/AliITSUSimulationPix.cxx
ITS/UPGRADE/AliITSUSimulationPix.h
ITS/UPGRADE/AliITSUv11.cxx
ITS/UPGRADE/testITSU/MakeITSUSimuParam.C [new file with mode: 0644]
ITS/UPGRADE/testITSU/sim.C

index 51fc5be..c315c23 100644 (file)
@@ -97,6 +97,9 @@ the AliITS class.
 #include "AliITSUSegmentationPix.h"
 #include "AliITSUSimuParam.h"
 #include "AliITSFOSignalsSPD.h"
+#include "AliParamList.h"
+#include "AliCDBManager.h" // tmp! Later the simuparam should be loaded centrally
+#include "AliCDBEntry.h"
 
 ClassImp(AliITSU)
 
@@ -117,8 +120,9 @@ AliDetector()
   ,fModuleHits(0)
   ,fDetDigits(0)
   ,fSensMap(0)
-  ,fSimulation(0)
-  ,fSegmentation(0)
+  ,fSimModelLr(0)
+  ,fSegModelLr(0)
+  ,fResponseLr(0)
   ,fCalibration(0)
   ,fRunNumber(0)
   ,fSimInitDone(kFALSE)
@@ -143,8 +147,9 @@ AliITSU::AliITSU(const Char_t *title, Int_t nlay) :
   ,fModuleHits(0)
   ,fDetDigits(0)
   ,fSensMap(0)
-  ,fSimulation(0)
-  ,fSegmentation(0)
+  ,fSimModelLr(0)
+  ,fSegModelLr(0)
+  ,fResponseLr(0)
   ,fCalibration(0)
   ,fRunNumber(0)
   ,fSimInitDone(kFALSE)
@@ -164,11 +169,25 @@ AliITSU::~AliITSU()
   // Default destructor for ITS.
   //  
   delete fHits;
-  delete fSimuParam;
+  //  delete fSimuParam; // provided by the CDBManager
   delete fSensMap;
-  if (fSimulation) fSimulation->Delete();
-  delete fSimulation;
-  delete fSegmentation;
+  if (fSimModelLr) {
+    for (int i=fNLayers;i--;) { // different layers may use the same simulation model
+      for (int j=i;j--;) if (fSimModelLr[j]==fSimModelLr[i]) fSimModelLr[j] = 0;
+      delete fSimModelLr[i];
+    }
+    delete[] fSimModelLr;
+  }
+  if (fSegModelLr) {
+    for (int i=fNLayers;i--;) { // different layers may use the same simulation model     
+      for (int j=i;j--;) if (fSegModelLr[j]==fSegModelLr[i]) fSegModelLr[j] = 0;
+      delete fSegModelLr[i];
+    }
+    delete[] fSegModelLr;
+  }
+  //
+  delete fResponseLr; // note: the response data is owned by the CDBManager, we don't delete them
+  //
   delete[] fLayerName;  // Array of TStrings
   delete[] fIdSens;
   //
@@ -486,11 +505,12 @@ void AliITSU::Hits2SDigits(Int_t evNumber,Int_t bgrev,Option_t *option,const cha
   FillModules(bgrev,option,filename);
   //
   Int_t nmodules = fGeomTGeo->GetNModules();
+  
   for(int module=0;module<nmodules;module++) {
-    int id = fGeomTGeo->GetModuleDetTypeID(module);
-    AliITSUSimulation* sim = GetSimulationModel(id);
-    if (!sim) AliFatal(Form("The sim.class for module %d of DetTypeID %d is missing",module,id));
-    sim->SDigitiseModule( GetModule(module), module, evNumber, GetSegmentation(id));
+    int lr = fGeomTGeo->GetLayer(module);
+    AliITSUSimulation* sim = GetSimulationModel(lr);
+    sim->InitSimulationModule(GetModule(module),evNumber/*,gAlice->GetEvNumber()*/,GetSegmentation(lr),GetResponseParam(lr));
+    sim->SDigitiseModule();
     fLoader->TreeS()->Fill();      // fills all branches - wasted disk space
     ResetSDigits();
   } 
@@ -544,10 +564,10 @@ void AliITSU::Hits2Digits(Int_t evNumber,Int_t bgrev,Option_t *option,const char
   // 
   Int_t nmodules = fGeomTGeo->GetNModules();
   for (Int_t module=0;module<nmodules;module++) {
-    int id = fGeomTGeo->GetModuleDetTypeID(module);
-    AliITSUSimulation* sim = GetSimulationModel(id);
-    if (!sim) AliFatal(Form("The sim.class for module %d of DetTypeID %d is missing",module,id));
-    sim->DigitiseModule( GetModule(module) ,module, evNumber, GetSegmentation(id));
+    int lr = fGeomTGeo->GetLayer(module);
+    AliITSUSimulation* sim = GetSimulationModel(lr);
+    sim->InitSimulationModule(GetModule(module),evNumber/*gAlice->GetEvNumber()*/,GetSegmentation(lr),GetResponseParam(lr));
+    sim->DigitiseModule();
     // fills all branches - wasted disk space
     fLoader->TreeD()->Fill(); 
     ResetDigits();
@@ -786,13 +806,9 @@ void AliITSU::SDigits2Digits()
   //
   int nmodules = fGeomTGeo->GetNModules();
   for (int module=0;module<nmodules;module++) {
-    int id = fGeomTGeo->GetModuleDetTypeID(module);
-    AliITSUSimulation *sim = GetSimulationModel(id);
-    if(!sim){
-      AliFatal(Form("The sim.class for module %d of DetTypeID %d is missing",module,id));
-      exit(1);
-    }
-    sim->InitSimulationModule(module,gAlice->GetEvNumber(),GetSegmentation(id));
+    int lr = fGeomTGeo->GetLayer(module);
+    AliITSUSimulation* sim = GetSimulationModel(lr);
+    sim->InitSimulationModule(GetModule(module),gAlice->GetEvNumber(),GetSegmentation(lr),GetResponseParam(lr));
     fSDigits->Clear();
     brchSDigits->GetEvent(module);
     sim->AddSDigitsToModule(fSDigits,0);
@@ -814,27 +830,51 @@ void AliITSU::InitSimulation()
   //
   if (fSimInitDone) {AliInfo("Already done"); return;}
   //
-  fSimuParam    = new AliITSUSimuParam();
+  AliCDBEntry* cdbEnt = AliCDBManager::Instance()->Get("ITS/Calib/SimuParam"); // tmp: load it centrally
+  if (!cdbEnt) {AliFatal("Failed to find ITS/Calib/SimuParam on CDB"); exit(1);}
+  fSimuParam    = (AliITSUSimuParam*)cdbEnt->GetObject();
+  //
   fSensMap      = new AliITSUSensMap("AliITSUSDigit",0,0);
-  fSimulation   = new TObjArray(kNDetTypes);
-  fSegmentation = new TObjArray();
-  AliITSUSegmentationPix::LoadSegmentations(fSegmentation, AliITSUGeomTGeo::GetITSsegmentationFileName());
-  fSegmentation->SetOwner(kTRUE);
+  fSimModelLr   = new AliITSUSimulation*[fNLayers];
+  fSegModelLr   = new AliITSsegmentation*[fNLayers];
+  fResponseLr   = new AliParamList*[fNLayers];
+  //
+  TObjArray arrSeg;
+  AliITSUSegmentationPix::LoadSegmentations(&arrSeg, AliITSUGeomTGeo::GetITSsegmentationFileName());
   //
   // add known simulation types used in the setup
   for (int i=fNLayers;i--;) {
-    int sType = fGeomTGeo->GetLayerDetTypeID(i)/AliITSUGeomTGeo::kMaxSegmPerDetType;
-    if (fSimulation->At(sType)) continue;
+    int dType = fGeomTGeo->GetLayerDetTypeID(i);           // fine detector type: class + segmentation
+    int sType = dType/AliITSUGeomTGeo::kMaxSegmPerDetType; // detector simulation class
     //
+    // check if the simulation of this sType was already created for preceeding layers
     AliITSUSimulation* simUpg = 0;
-    switch (sType) {
-    case AliITSUGeomTGeo::kDetTypePix : 
-      simUpg = new AliITSUSimulationPix(fSimuParam,fSensMap); 
-      break;
-    default: AliFatal(Form("No %d detector type is defined",sType));
-    };
-    fSimulation->AddAtAndExpand(simUpg,sType);
+    for (int j=fNLayers;j>i;j--) {
+      simUpg = GetSimulationModel(j);
+      if (int(simUpg->GetUniqueID())==sType) break;
+      else simUpg = 0;
+    }
+    //
+    if (!simUpg) { // need to create simulation for detector class sType
+      switch (sType) 
+       {
+       case AliITSUGeomTGeo::kDetTypePix : 
+         simUpg = new AliITSUSimulationPix(fSimuParam,fSensMap); 
+         break;
+       default: AliFatal(Form("No %d detector type is defined",sType));
+       }
+    }
+    fSimModelLr[i] = simUpg;
+    //
+    // add segmentations used in the setup
+    if (!(fSegModelLr[i]=(AliITSsegmentation*)arrSeg[dType])) {AliFatal(Form("Segmentation for DetType#%d is not found",dType)); exit(1);}
+    //
+    // add response function for the detectors of this layer
+    if ( !(fResponseLr[i]=(AliParamList*)fSimuParam->FindRespFunParams(dType)) ) {AliFatal(Form("Response for DetType#%d is not found in SimuParams",dType)); exit(1);}
   }
+  // delete non needed segmentations
+  for (int i=fNLayers;i--;) arrSeg.Remove(fSegModelLr[i]);
+  arrSeg.Delete();
   //
   InitArrays();
   //
index 867e0fb..5be2487 100644 (file)
@@ -28,7 +28,7 @@ class AliITSdigit;
 class AliDigitizationInput;
 class AliITSUSensMap;
 class AliITSUSimuParam;
-
+class AliParamList;
 
 class AliITSU : public AliDetector {
 
@@ -69,8 +69,9 @@ class AliITSU : public AliDetector {
   virtual void MakeBranchD(const char* file);
   virtual void MakeBranchInTreeD(TTree* treeD, const char* file=0);
   virtual void SetTreeAddress();
-  virtual AliITSUSimulation* GetSimulationModel(Int_t id) {return (AliITSUSimulation*)fSimulation->At(id/AliITSUGeomTGeo::kMaxSegmPerDetType);}
-  virtual AliITSsegmentation*  GetSegmentation(Int_t id)    {return (AliITSsegmentation*) fSegmentation->At(id);}
+  virtual AliITSUSimulation*   GetSimulationModel(Int_t lr)   {return (AliITSUSimulation*)fSimModelLr[lr];}
+  virtual AliITSsegmentation*  GetSegmentation(Int_t lr)      {return (AliITSsegmentation*)fSegModelLr[lr];}
+  virtual AliParamList*        GetResponseParam(Int_t lr)     {return (AliParamList*)fResponseLr[lr];}
   //=================== Hits =========================================
   virtual void StepManager() {} // See Step Manager for specific geometry.
   //------------ sort hits by module for Digitisation ----------------
@@ -139,8 +140,9 @@ class AliITSU : public AliDetector {
   TObjArray*            fDetDigits;      //! AliDetector has TClonesArray fDigits, avoid same name
   AliITSUSensMap*       fSensMap;        //! sensor map for digitization
   //
-  TObjArray            *fSimulation;     //! simulation objects per det.type
-  TObjArray            *fSegmentation;   //! segmentation objects per det.type (and segmentation)
+  AliITSUSimulation    **fSimModelLr;     //! simulation objects per layer
+  AliITSsegmentation   **fSegModelLr;     //! segmentation objects per layar
+  AliParamList         **fResponseLr;     //! response parameters for each layer
   TObjArray            *fCalibration;    //! calibration objects
   Int_t                 fRunNumber;      //! run number
   Bool_t                fSimInitDone;    //! flag initialized simulation
index e92ca03..2b30a0a 100644 (file)
@@ -154,12 +154,12 @@ void AliITSUDigitizer::Digitize(Option_t* /*opt*/)
   for (int module=0; module<nModules; module++ ) {
     //
     if (!fRoif && !fModActive[module]) continue;
-    int id = geom->GetModuleDetTypeID(module);
-    AliITSUSimulation *sim = fITS->GetSimulationModel(id);
-    if (!sim) AliFatal(Form("The simulation model %d is not available",id));
+    int lr = geom->GetLayer(module);
+    AliITSUSimulation *sim = fITS->GetSimulationModel(lr);
+    if (!sim) AliFatal(Form("The simulation model for layer %d is not available",lr));
     //
     // Fill the module with the sum of SDigits
-    sim->InitSimulationModule(module, event, fITS->GetSegmentation(id));
+    sim->InitSimulationModule(fITS->GetModule(module), event, fITS->GetSegmentation(lr), fITS->GetResponseParam(lr));
     //
     for (int ifiles=0; ifiles<nfiles; ifiles++ ) {
       //
index 05bb208..4be281a 100644 (file)
@@ -105,15 +105,21 @@ Bool_t AliITSURecoDet::Build()
   }
   //
   // TPC-ITS wall
-  lrp = new AliITSURecoLayer("TPC-ITSwall");
-  lrp->SetRMin(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMin());
-  lrp->SetRMax(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMax());
-  lrp->SetR(0.5*(lrp->GetRMin()+lrp->GetRMax()));
-  lrp->SetZMin(-AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
-  lrp->SetZMax( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
-  lrp->SetMaxStep( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallMaxStep());
-  lrp->SetPassive(kTRUE);
-  AddLayer(lrp);
+  const AliITSURecoParam* recopar = AliITSUReconstructor::GetRecoParam();
+  if (recopar) {
+    lrp = new AliITSURecoLayer("TPC-ITSwall");
+    lrp->SetRMin(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMin());
+    lrp->SetRMax(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMax());
+    lrp->SetR(0.5*(lrp->GetRMin()+lrp->GetRMax()));
+    lrp->SetZMin(-AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
+    lrp->SetZMax( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
+    lrp->SetMaxStep( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallMaxStep());
+    lrp->SetPassive(kTRUE);
+    AddLayer(lrp);
+  }
+  else {
+    AliWarning("RecoParam is not available, TPC-ITS wall is not set");
+  }
   //
   IndexLayers();
   Print("lr");
index 37443be..63de282 100644 (file)
@@ -24,6 +24,8 @@
 ///////////////////////////////////////////////////////////////////
 #include "AliITSUSimuParam.h"
 #include "AliLog.h"
+#include "AliParamList.h"
+
 using namespace TMath;
 
 
@@ -74,8 +76,10 @@ AliITSUSimuParam::AliITSUSimuParam()
   ,fPixSigma(0)
   ,fPixNoise(0)
   ,fPixBaseline(0)
+  ,fRespFunParam(0)
 {  
   // default constructor
+  fRespFunParam.SetOwner(kTRUE);
 }
 
 //______________________________________________________________________
@@ -107,6 +111,7 @@ AliITSUSimuParam::AliITSUSimuParam(UInt_t nPix)
   ,fPixSigma(0)
   ,fPixNoise(0)
   ,fPixBaseline(0)
+  ,fRespFunParam(0)
 {  
   // regular constructor
   if (fNPix>0) {
@@ -119,6 +124,7 @@ AliITSUSimuParam::AliITSUSimuParam(UInt_t nPix)
   SetPixThreshold(fgkPixThreshDefault,fgkPixThrSigmaDefault);
   SetPixNoise(0.,0.);
   SetPixBiasVoltage(fgkPixBiasVoltageDefault);
+  fRespFunParam.SetOwner(kTRUE);
   //
 }
 
@@ -152,6 +158,7 @@ AliITSUSimuParam::AliITSUSimuParam(const AliITSUSimuParam &simpar)
   ,fPixSigma(0)
   ,fPixNoise(0)
   ,fPixBaseline(0)   
+  ,fRespFunParam(0)
    //
 {
   // copy constructor
@@ -170,6 +177,11 @@ AliITSUSimuParam::AliITSUSimuParam(const AliITSUSimuParam &simpar)
     fPixBaseline[i]    = simpar.fPixBaseline[i];
   }
   //
+  for (int i=0;i<simpar.fRespFunParam.GetEntriesFast();i++) {
+    AliParamList* pr = (AliParamList*)simpar.fRespFunParam[0];
+    if (pr) fRespFunParam.AddLast(new AliParamList(*pr));
+  }
+  fRespFunParam.SetOwner(kTRUE);
 }
 
 //______________________________________________________________________
@@ -420,71 +432,20 @@ Double_t AliITSUSimuParam::LorentzAngleElectron(Double_t B) const
   return angle;
 }
 
-//______________________________________________________________________
-Double_t AliITSUSimuParam::SigmaDiffusion3D(Double_t l) const 
+//___________________________________________________________
+const AliParamList* AliITSUSimuParam::FindRespFunParams(Int_t detId) const
 {
-  // Returns the Gaussian sigma^2 == <x^2+y^2+z^2> [cm^2] due to the
-  // defusion of electrons or holes through a distance l [cm] caused
-  // by an applied voltage v [volt] through a distance d [cm] in any
-  //  material at a temperature T [degree K]. The sigma diffusion when
-  //  expressed in terms of the distance over which the diffusion
-  // occures, l=time/speed, is independent of the mobility and therefore
-  //  the properties of the material. The charge distributions is given by
-  // n = exp(-r^2/4Dt)/(4piDt)^1.5. From this <r^2> = 6Dt where D=mkT/e
-  // (m==mobility, k==Boltzman's constant, T==temparature, e==electric
-  // charge. and vel=m*v/d. consiquently sigma^2=6kTdl/ev.
-  // Inputs:
-  //    Double_t l   Distance the charge has to travel.
-  // Output:
-  //    none.
-  // Return:
-  //    The Sigma due to the diffution of electrons. [cm]
-  const Double_t kcon = 5.17040258E-04; // == 6k/e [J/col or volts]  
-  return Sqrt(kcon*fT*fDOverV*l);  // [cm]
-}
-
-//______________________________________________________________________
-Double_t AliITSUSimuParam::SigmaDiffusion2D(Double_t l) const 
-{
-  // Returns the Gaussian sigma^2 == <x^2+z^2> [cm^2] due to the defusion
-  // of electrons or holes through a distance l [cm] caused by an applied
-  // voltage v [volt] through a distance d [cm] in any material at a
-  // temperature T [degree K]. The sigma diffusion when expressed in terms
-  // of the distance over which the diffusion occures, l=time/speed, is
-  // independent of the mobility and therefore the properties of the
-  // material. The charge distributions is given by
-  // n = exp(-r^2/4Dt)/(4piDt)^1.5. From this <x^2+z^2> = 4Dt where D=mkT/e
-  // (m==mobility, k==Boltzman's constant, T==temparature, e==electric
-  // charge. and vel=m*v/d. consiquently sigma^2=4kTdl/ev.
-  // Inputs:
-  //    Double_t l   Distance the charge has to travel.
-  // Output:
-  //    none.
-  // Return:
-  //    The Sigma due to the diffution of electrons. [cm]
-  const Double_t kcon = 3.446935053E-04; // == 4k/e [J/col or volts]
-  return Sqrt(kcon*fT*fDOverV*l);  // [cm]
+  // find parameters list for detID
+  for (int i=fRespFunParam.GetEntriesFast();i--;) {
+    const AliParamList* pr = GetRespFunParams(i);
+    if (int(pr->GetUniqueID())==detId) return pr;
+  }
+  return 0;
 }
 
-//______________________________________________________________________
-Double_t AliITSUSimuParam::SigmaDiffusion1D(Double_t l) const 
+//___________________________________________________________
+void AliITSUSimuParam::AddRespFunParam(AliParamList* pr)
 {
-  // Returns the Gaussian sigma^2 == <x^2> [cm^2] due to the defusion
-  // of electrons or holes through a distance l [cm] caused by an applied
-  // voltage v [volt] through a distance d [cm] in any material at a
-  // temperature T [degree K]. The sigma diffusion when expressed in terms
-  // of the distance over which the diffusion occures, l=time/speed, is
-  // independent of the mobility and therefore the properties of the
-  // material. The charge distributions is given by
-  // n = exp(-r^2/4Dt)/(4piDt)^1.5. From this <r^2> = 2Dt where D=mkT/e
-  // (m==mobility, k==Boltzman's constant, T==temparature, e==electric
-  // charge. and vel=m*v/d. consiquently sigma^2=2kTdl/ev.
-  // Inputs:
-  //    Double_t l   Distance the charge has to travel.
-  // Output:
-  //    none.
-  // Return:
-  //    The Sigma due to the diffution of electrons. [cm]
-  const Double_t kcon = 1.723467527E-04; // == 2k/e [J/col or volts]
-  return Sqrt(kcon*fT*fDOverV*l);  // [cm]
+  // add spread parameterization data
+  fRespFunParam.AddLast(pr);
 }
index 53b7ca8..83617be 100644 (file)
 ///////////////////////////////////////////////////////////////////
 #include <TRandom.h>
 #include <TObject.h>
+#include <TObjArray.h>
 #include <TMath.h>
 #include "AliMathBase.h"
+class AliParamList;
 
 class AliITSUSimuParam : public TObject {
 
@@ -46,19 +48,12 @@ class AliITSUSimuParam : public TObject {
   Double_t GetGeVToCharge()                                                  const  {return fGeVcharge;}
   Double_t GeVToCharge(Double_t gev)                                         const  {return gev/fGeVcharge;}
   //
-  void     SetDistanceOverVoltage(Double_t d,Double_t v)                            {fDOverV = d/v;}
-  void     SetDistanceOverVoltage(Double_t dv=fgkDOverVDefault)                     {fDOverV = dv;}
-  Double_t GetDistanceOverVoltage()                                          const  {return fDOverV;}
-  //
   void     SetPixCouplingOption(UInt_t opt);
   UInt_t   GetPixCouplingOption()                                         const  {return fPixCouplOpt;}
 
   void     SetPixCouplingParam(Double_t col, Double_t row)                       {fPixCouplCol = col; fPixCouplRow = row;}
   void     GetPixCouplingParam(Double_t &col, Double_t &row)              const  {col = fPixCouplCol; row = fPixCouplRow;}
 
-  void     SetPixSigmaDiffusionAsymmetry(Double_t ecc)                           {fPixEccDiff=ecc;}   
-  void     GetPixSigmaDiffusionAsymmetry(Double_t &ecc)                   const  {ecc=fPixEccDiff;}
-
   void     SetPixLorentzDrift(Bool_t ison)                                       {fPixLorentzDrift=ison;}
   Bool_t   GetPixLorentzDrift()                                           const  {return fPixLorentzDrift;}
   void     SetPixLorentzHoleWeight(Double_t weight)                               {fPixLorentzHoleWeight=weight;}
@@ -76,6 +71,11 @@ class AliITSUSimuParam : public TObject {
   Double_t SigmaDiffusion2D(Double_t l)                                      const;
   Double_t SigmaDiffusion1D(Double_t l)                                      const;
   //
+  Int_t    GetNRespFunParams()                                               const {return fRespFunParam.GetEntriesFast();}
+  const AliParamList* GetRespFunParams(Int_t i)                              const {return (const AliParamList*)fRespFunParam[i];}
+  const AliParamList* FindRespFunParams(Int_t detId)                         const;
+  void     AddRespFunParam(AliParamList* pr);
+  //
   virtual void Print(Option_t *opt = "")                                     const; 
   //
   static Double_t CalcProbNoiseOverThreshold(double base, double noise, double thresh);
@@ -130,6 +130,7 @@ class AliITSUSimuParam : public TObject {
   Double_t*  fPixNoise;       //[fNPix] Pix electronic noise: sigma
   Double_t*  fPixBaseline;    //[fNPix] Pix electronic noise: baseline
   //
+  TObjArray  fRespFunParam;   // set of parameterizations for response function (AliParamList)
 
   ClassDef(AliITSUSimuParam,1);  // ITSU simulataion params
 };
index 5a6a271..889ebac 100644 (file)
@@ -22,6 +22,8 @@
 #include "TSeqCollection.h"
 #include "AliITSUSimulation.h"
 #include "AliITSUSDigit.h"
+#include "AliITSUModule.h"
+#include "AliParamList.h"
 using namespace TMath;
 
 ClassImp(AliITSUSimulation)
@@ -33,6 +35,7 @@ AliITSUSimulation::AliITSUSimulation()
   ,fCalibNoisy(0)
   ,fSensMap(0)
   ,fSimuParam(0)
+  ,fResponseParam(0)
   ,fModule(0)
   ,fEvent(0)
   ,fDebug(0)
@@ -46,6 +49,7 @@ AliITSUSimulation::AliITSUSimulation(AliITSUSimuParam* sim,AliITSUSensMap* map)
   ,fCalibNoisy(0)
   ,fSensMap(map)
   ,fSimuParam(sim)
+  ,fResponseParam(0)
   ,fModule(0)
   ,fEvent(0)
   ,fDebug(0)
@@ -61,6 +65,7 @@ AliITSUSimulation::AliITSUSimulation(const AliITSUSimulation &s)
   ,fCalibNoisy(s.fCalibNoisy)
   ,fSensMap(s.fSensMap)
   ,fSimuParam(s.fSimuParam)   
+  ,fResponseParam(s.fResponseParam)
   ,fModule(s.fModule)
   ,fEvent(s.fEvent)
   ,fDebug(s.fDebug)
@@ -78,21 +83,26 @@ AliITSUSimulation&  AliITSUSimulation::operator=(const AliITSUSimulation &s)
   fCalibNoisy= s.fCalibNoisy;
   fSensMap   = s.fSensMap;
   fSimuParam = s.fSimuParam;
+  fResponseParam = s.fResponseParam;
   fModule    = s.fModule;
   fEvent     = s.fEvent;
   return *this;
 }
 
 //______________________________________________________________________
-void AliITSUSimulation::InitSimulationModule(Int_t module, Int_t event, AliITSsegmentation* seg)
+void AliITSUSimulation::InitSimulationModule(AliITSUModule* mod, Int_t event, AliITSsegmentation* seg, AliParamList* resp)
 {
   //  This function creates maps to build the list of tracks for each
   //  summable digit. Inputs defined by base class.
   //
-  SetModule(module);
-  SetEvent(event);
+  SetModule(mod);
   SetSegmentation(seg);
+  SetResponseParam(resp);
   ClearMap();
+  //
+  if (event != fEvent) GenerateStrobePhase(); 
+  SetEvent(event);
+  
 }
 
 //______________________________________________________________________
@@ -108,7 +118,7 @@ Bool_t AliITSUSimulation::AddSDigitsToModule(TSeqCollection *pItemArr,Int_t mask
   // 
   for( Int_t i=0; i<nItems; i++ ) {
     AliITSUSDigit * pItem = (AliITSUSDigit *)(pItemArr->At( i ));
-    if( pItem->GetModule() != fModule ) AliFatal(Form("SDigits module %d != current module %d: exit", pItem->GetModule(), fModule ));
+    if(pItem->GetModule() != int(fModule->GetIndex()) ) AliFatal(Form("SDigits module %d != current module %d: exit", pItem->GetModule(),fModule->GetIndex()));
     if(pItem->GetSumSignal()>0.0 ) sig = kTRUE;
     AliITSUSDigit* oldItem = (AliITSUSDigit*)fSensMap->GetItem(pItem);
     if (!oldItem) {
@@ -126,7 +136,7 @@ void AliITSUSimulation::UpdateMapSignal(UInt_t dim0,UInt_t dim1,Int_t trk,Int_t
   // update map with new hit
   UInt_t ind = fSensMap->GetIndex(dim0,dim1);
   AliITSUSDigit* oldItem = (AliITSUSDigit*)fSensMap->GetItem(ind);
-  if (!oldItem) fSensMap->RegisterItem( new(fSensMap->GetFree()) AliITSUSDigit(trk,ht,fModule,ind,signal) );
+  if (!oldItem) fSensMap->RegisterItem( new(fSensMap->GetFree()) AliITSUSDigit(trk,ht,fModule->GetIndex(),ind,signal) );
   else oldItem->AddSignal(trk,ht,signal);
 }
 
@@ -136,7 +146,7 @@ void AliITSUSimulation::UpdateMapNoise(UInt_t dim0,UInt_t dim1,Double_t noise)
   // update map with new hit
   UInt_t ind = fSensMap->GetIndex(dim0,dim1);
   AliITSUSDigit* oldItem = (AliITSUSDigit*)fSensMap->GetItem(ind);
-  if (!oldItem) fSensMap->RegisterItem( new(fSensMap->GetFree()) AliITSUSDigit(fModule,ind,noise) );
+  if (!oldItem) fSensMap->RegisterItem( new(fSensMap->GetFree()) AliITSUSDigit(fModule->GetIndex(),ind,noise) );
   else oldItem->AddNoise(noise);
 }
 
index d0df46c..5eb31fc 100644 (file)
 #include <TObject.h>
 #include "AliITSUSensMap.h"
 #include "AliITSsegmentation.h"
+#include "AliMathBase.h"
 
 class AliITSCalibration;
 class AliITSUSimuParam;
 class AliITSUModule;
 class TRandom;
 class TSegCollection;
+class AliParamList;
 
 // This is the base class for ITS detector signal simulations. Data members
 // include are a pointer to the detectors specific response and segmentation
@@ -45,33 +47,38 @@ class AliITSUSimulation : public TObject
   //
   void UpdateMapSignal(UInt_t dim0,UInt_t dim1, Int_t trk,Int_t ht,Double_t signal);
   void UpdateMapNoise(UInt_t dim0,UInt_t dim1, Double_t noise);
-  virtual void InitSimulationModule(Int_t, Int_t, AliITSsegmentation*);
+  virtual void InitSimulationModule(AliITSUModule* mod, Int_t ev, AliITSsegmentation* seg, AliParamList* resp);
   //
   // Hits -> SDigits
-  virtual void SDigitiseModule(AliITSUModule* mod, Int_t mask, Int_t event, AliITSsegmentation* seg) = 0;
+  virtual void SDigitiseModule() = 0;
   virtual void FinishSDigitiseModule() = 0;
   virtual Bool_t AddSDigitsToModule( TSeqCollection *pItemArray, Int_t mask );
   //
   // Hits -> Digits
-  virtual void DigitiseModule(AliITSUModule* mod, Int_t mask, Int_t event, AliITSsegmentation* seg) = 0;
+  virtual void DigitiseModule() = 0;
   virtual void CreateFastRecPoints(AliITSUModule *,Int_t,TRandom *,TClonesArray* /*recp*/) {}
   //
+  // readout phase (strobe, timing etc) generation
+  virtual void GenerateStrobePhase() {}
+  //
   AliITSCalibration*  GetCalibDead()                   const {return fCalibDead;}
   AliITSCalibration*  GetCalibNoisy()                  const {return fCalibNoisy;}
   AliITSsegmentation* GetSegmentation()                const {return fSeg;}
-  AliITSUSimuParam* GetSimuParam()                   const {return fSimuParam;}
-  AliITSUSensMap*   GetMap()                         const {return fSensMap;}
-  Int_t               GetModule()                      const {return fModule;}
+  AliITSUSimuParam*   GetSimuParam()                   const {return fSimuParam;}
+  AliITSUSensMap*     GetMap()                         const {return fSensMap;}
+  AliITSUModule*      GetModule()                      const {return fModule;}
+  AliParamList*       GetResponseParam()               const {return fResponseParam;}
   Int_t               GetEvent()                       const {return fEvent;}
   Bool_t              GetDebug(Int_t level=1)          const {return fDebug>=level;}
-
+  
   //
   void SetCalibDead(AliITSCalibration *calib)              {fCalibDead = calib;}
   void SetCalibNoisy(AliITSCalibration *calib)             {fCalibNoisy = calib;}
   void SetSegmentation(AliITSsegmentation *seg)            {fSeg = seg; if (seg&&fSensMap) fSensMap->SetDimensions(seg->Npz(),seg->Npx());}
-  void SetSimuParam(AliITSUSimuParam *sp)                {fSimuParam = sp;}
-  void SetMap(AliITSUSensMap *p)                         {fSensMap = p;}
-  void SetModule(Int_t mod)                                {fModule=mod;} 
+  void SetSimuParam(AliITSUSimuParam *sp)                  {fSimuParam = sp;}
+  virtual void SetResponseParam(AliParamList* resp)        {fResponseParam = resp;}
+  void SetMap(AliITSUSensMap *p)                           {fSensMap = p;}
+  void SetModule(AliITSUModule* mod)                       {fModule=mod;} 
   void SetEvent(Int_t evnt)                                {fEvent=evnt;} 
   void SetDebug(Int_t level=5)                             {fDebug=level;}
   void SetNoDebug()                                        {fDebug=0;}
@@ -79,13 +86,18 @@ class AliITSUSimulation : public TObject
   //
   static  Int_t GenOrderedSample(UInt_t nmax,UInt_t ngen,TArrayI &vals,TArrayI &ind);
   //
+  static  Double_t GausInt1D(Double_t sig,Double_t a,Double_t b);
+  static  Double_t GausInt2D(Double_t sig0,Double_t a0,Double_t b0,
+                            Double_t sig1,Double_t a1,Double_t b1);
+  //
  protected:
   AliITSsegmentation  *fSeg;            //! segmentation
   AliITSCalibration   *fCalibDead;      //! dead channels
   AliITSCalibration   *fCalibNoisy;     //! noisy channels
-  AliITSUSensMap    *fSensMap;        //! sensor map for hits manipulations
-  AliITSUSimuParam  *fSimuParam;      //! simulation parameters
-  Int_t                fModule;         //! module number being processed
+  AliITSUSensMap      *fSensMap;        //! sensor map for hits manipulations
+  AliITSUSimuParam    *fSimuParam;      //! simulation parameters
+  AliParamList        *fResponseParam;  //! response parameterization data
+  AliITSUModule       *fModule;         //! module being processed
   Int_t                fEvent;          //! event number being processed
   Int_t                fDebug;          //!  debug flag
 
@@ -93,4 +105,25 @@ class AliITSUSimulation : public TObject
     
 };
 
+//_____________________________________________________________________________
+inline Double_t AliITSUSimulation::GausInt1D(Double_t sig,Double_t a,Double_t b)
+{
+  // calculate gaussian integral from a to b (with respecto to mean)
+  const Double_t kRoot2 = 1.414213562; // Sqrt(2).
+  double sp = 1.0/(sig*kRoot2);
+  return 0.5*TMath::Abs( AliMathBase::ErfcFast(sp*a) - AliMathBase::ErfcFast(sp*b) );
+}
+
+//_____________________________________________________________________________
+inline Double_t AliITSUSimulation::GausInt2D(Double_t sig0,Double_t a0,Double_t b0,
+                                            Double_t sig1,Double_t a1,Double_t b1)
+{
+  // calculate gaussian 2D integral from x0 to x1 (with respect to mean)
+  const Double_t kRoot2 = 1.414213562; // Sqrt(2).
+  double sp0 = 1.0/(sig0*kRoot2);
+  double sp1 = 1.0/(sig1*kRoot2);
+  return 0.25*TMath::Abs( (AliMathBase::ErfcFast(sp0*a0) - AliMathBase::ErfcFast(sp0*b0)) *
+                         (AliMathBase::ErfcFast(sp1*a1) - AliMathBase::ErfcFast(sp1*b1)));
+}
+
 #endif
index 1bff3da..74acb4e 100644 (file)
@@ -36,6 +36,7 @@
 #include "AliMathBase.h"
 #include "AliITSUSimuParam.h"
 #include "AliITSUSDigit.h"
+#include "AliParamList.h"
 
 using std::cout;
 using std::endl;
@@ -65,8 +66,10 @@ AliITSUSimulationPix::AliITSUSimulationPix()
   ,fStrobe(kTRUE)
   ,fStrobeLenght(4)
   ,fStrobePhase(-12.5e-9)
+  ,fSpreadFun(0)
 {
    // Default constructor.
+  SetUniqueID(AliITSUGeomTGeo::kDetTypePix);
 }
 
 //______________________________________________________________________
@@ -76,8 +79,10 @@ AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap*
   ,fStrobe(kTRUE)
   ,fStrobeLenght(4)
   ,fStrobePhase(-12.5e-9)
+  ,fSpreadFun(0)
 {
   // standard constructor
+  SetUniqueID(AliITSUGeomTGeo::kDetTypePix);
   Init();
 }
 
@@ -88,6 +93,7 @@ AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
   ,fStrobe(s.fStrobe)
   ,fStrobeLenght(s.fStrobeLenght)
   ,fStrobePhase(s.fStrobePhase)
+  ,fSpreadFun(s.fSpreadFun)
 {
   //     Copy Constructor
 }
@@ -109,6 +115,8 @@ AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix
   fStrobe       = s.fStrobe;
   fStrobeLenght = s.fStrobeLenght;
   fStrobePhase  = s.fStrobePhase;
+  fSpreadFun    = s.fSpreadFun;
+  //
   return *this;
 }
 
@@ -147,19 +155,17 @@ Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
 }
 
 //_____________________________________________________________________
-void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg)
+void AliITSUSimulationPix::SDigitiseModule()
 {
   //  This function begins the work of creating S-Digits.
-  if (!(mod->GetNHits())) {
+  if (!(fModule->GetNHits())) {
     AliDebug(1,Form("In event %d module %d there are %d hits returning.",
-                   event, mod->GetIndex(),mod->GetNHits()));
+                   fEvent, fModule->GetIndex(),fModule->GetNHits()));
     return;
   } 
   //
-  if (fStrobe) if (event != fEvent) GenerateStrobePhase(); 
-  InitSimulationModule(mod->GetIndex(), event, seg );
-  // Hits2SDigits(mod);
-  Hits2SDigitsFast(mod);
+  Hits2SDigitsFast();
+  //
   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();  
   WriteSDigits();
@@ -190,16 +196,14 @@ void AliITSUSimulationPix::FinishSDigitiseModule()
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg)
+void AliITSUSimulationPix::DigitiseModule()
 {
   //  This function creates Digits straight from the hits and then adds
   //  electronic noise to the digits before adding them to pList
   //  Each of the input variables is passed along to Hits2SDigits
   //
-  if (fStrobe) if (event != fEvent) GenerateStrobePhase();
-  InitSimulationModule( mod->GetIndex(), event, seg );
-  // Hits2SDigits(mod);
-  Hits2SDigitsFast(mod);
+  // pick charge spread function
+  Hits2SDigitsFast();
   //
   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
@@ -208,58 +212,54 @@ void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
+void AliITSUSimulationPix::Hits2SDigits()
 {
   // Does the charge distributions using Gaussian diffusion charge charing.
   const Double_t kBunchLenght = 25e-9; // LHC clock
-  Int_t nhits = mod->GetNHits();
+  Int_t nhits = fModule->GetNHits();
   if (!nhits) return;
   //
   Int_t h,ix,iz,i;
   Int_t idtrack;
   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
-  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;
-  Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
+  Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
+  Double_t t,tp,st,dt=0.2,el;
   Double_t thick = 0.5*fSeg->Dy();  // Half Thickness
-  fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
+
   //
   for (h=0;h<nhits;h++) {
     if (fStrobe && 
-       ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
-        (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
+       ((fModule->GetHit(h)->GetTOF()<fStrobePhase) || 
+        (fModule->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
        ) continue;
     //
-    if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
+    if (!fModule->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
     st = Sqrt(x1*x1+y1*y1+z1*z1);
     if (st>0.0) {
       st = (Double_t)((Int_t)(st*1e4)); // number of microns
       if (st<=1.0) st = 1.0;
       dt = 1.0/st;               // RS TODO: do we need 1 micron steps?
+      double dy = dt*thick;
+      y = -0.5*dy;
       for (t=0.0;t<1.0;t+=dt) { // Integrate over t
        tp  = t+0.5*dt;
        x   = x0+x1*tp;
-       y   = y0+y1*tp;
+       y  += dy;
        z   = z0+z1*tp;
        if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
        el  = dt * de / fSimuParam->GetGeVToCharge();
        //
-       sig = fSimuParam->SigmaDiffusion1D(Abs(thick + y)); 
-       sigx=sig;
-       sigz=sig*fda;
-       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+       SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
       } // end for t
     } else { // st == 0.0 deposit it at this point
       x   = x0;
-      y   = y0;
+      y   = y0 + 0.5*thick;
       z   = z0;
       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
       el  = de / fSimuParam->GetGeVToCharge();
-      sig = fSimuParam->SigmaDiffusion1D(Abs(thick + y));
-      sigx = sig;
-      sigz = sig*fda;
-      if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-      SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+      if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+      SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
     } // end if st>0.0    
   } // Loop over all hits h
   //
@@ -289,28 +289,13 @@ void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
+void AliITSUSimulationPix::Hits2SDigitsFast()
 {
   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
   //    AliITSUModule *mod  Pointer to this module
-  // Output:
-  //    none.
-  // Return:
-  //    none.
-  /* 
-  // RSTmp this injection points partitioning should be reimplemented
-  const Int_t kn10=10;
-  const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
-                           4.325316833e-1,4.869532643e-1,5.130467358e-1,
-                           5.674683167e-1,6.602952159e-1,7.833023029e-1,
-                           9.255628306e-1};
-  const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
-                           7.472567455e-2,3.333567215e-2,3.333567215e-2,
-                           7.472567455e-2,1.095431813e-1,1.346333597e-1,
-                           1.477621124e-1};
-  */
+  //
   const Double_t kBunchLenght = 25e-9; // LHC clock
-  TObjArray *hits = mod->GetHits();
+  TObjArray *hits = fModule->GetHits();
   Int_t nhits = hits->GetEntriesFast();
   if (nhits<=0) return;
   //
@@ -318,54 +303,46 @@ void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
   Int_t idtrack;
   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
   Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; 
-  Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
-  Double_t thick = 0.5*fSeg->Dy();  // Half thickness
+  Double_t t,st,el,de=0.0;
   Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
-  fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
+  Double_t thick = fSeg->Dy();
   //
   for (h=0;h<nhits;h++) {
     // Check if the hit is inside readout window
     if (fStrobe && 
-       ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
-        (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
+       ((fModule->GetHit(h)->GetTOF()<fStrobePhase) || 
+        (fModule->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
        ) continue;
     //
-    if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
+    if (!fModule->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
     st = Sqrt(x1*x1+y1*y1+z1*z1); 
     if (st>0.0) {
       int np = int(1.5*st/minDim);  //RStmp: at the moment neglect kti,kwi: inject the points in such a way that there is ~1.5 point per cell
-      if (np<5) np = 5;             //RStmp
-      double dt = 1./(np+1);        //RStmp
-      double dw = 1./np;
-      //      printf("Dst: %f md:%f np=%d Tr#%d\n",st,minDim,np,idtrack);
+      if (np<3) np = 3;             //RStmp
+      double dt = 1./np;
+      double dy = dt*thick;
+      y = -0.5*dy;
       t = -0.5*dt;
       for (i=0;i<np;i++) {          //RStmp Integrate over t
        //      for (i=0;i<kn10;i++) { // Integrate over t
        t  += dt;  // RStmp kti[i];
        x   = x0+x1*t;
-       y   = y0+y1*t;
+       y  += dy;
        z   = z0+z1*t;
        if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
-       //      el  = kwi[i]*de/fSimuParam->GetGeVToCharge();
-       el  = dw*de/fSimuParam->GetGeVToCharge();
-       sig = fSimuParam->SigmaDiffusion1D(Abs(thick + y));
-       sigx=sig;
-       sigz=sig*fda;
-       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+       el  = dt*de/fSimuParam->GetGeVToCharge();
+       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+       SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
       } // end for i // End Integrate over t
     }
     else { // st == 0.0 deposit it at this point
       x   = x0;
-      y   = y0;
+      y   = y0+0.5*thick;
       z   = z0;
       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
       el  = de / fSimuParam->GetGeVToCharge();
-      sig = fSimuParam->SigmaDiffusion1D(Abs(thick + y));
-      sigx=sig;
-      sigz=sig*fda;
-      if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-      SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+      if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+      SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
     } // end if st>0.0
     
   } // Loop over all hits h
@@ -394,139 +371,109 @@ void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
-                                         Int_t ix0,Int_t iz0,
-                                         Double_t el,Double_t sig,Double_t ld,
-                                         Int_t t,Int_t hi)
-{
-   // Spreads the charge over neighboring cells. Assume charge is distributed
-   // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
-   // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
-   // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
-   // Defined this way, the integral over all x and z is el.
-   // Inputs:
-   //    Double_t x0   x position of point where charge is liberated
-   //    Double_t z0   z position of point where charge is liberated
-   //    Int_t    ix0  row of cell corresponding to point x0
-   //    Int_t    iz0  columb of cell corresponding to point z0
-   //    Double_t el   number of electrons liberated in this step
-   //    Double_t sig  Sigma difusion for this step (y0 dependent)
-   //    Double_t ld   lorentz drift in x for this step (y0 dependent)
-   //    Int_t    t    track number
-   //    Int_t    ti   hit track index number
-   //    Int_t    hi   hit "hit" index number
-   // Outputs:
-   //     none.
-   // Return:
-   //     none.
-   const Int_t knx = 3,knz = 2;
-   const Double_t kRoot2 = 1.414213562; // Sqrt(2).
-   Int_t ix,iz,ixs,ixe,izs,ize;
-   Float_t x,z;  // keep coordinates float (required by AliSegmentation)
-   Double_t s,sp,x1,x2,z1,z2; 
-   //
-   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));
-   if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
-     UpdateMapSignal(iz0,ix0,t,hi,el);
-     return;
-   } // end if
-   sp = 1.0/(sig*kRoot2);
-   ixs = Max(-knx+ix0,0);
-   ixe = Min(knx+ix0,fSeg->Npx()-1);
-   izs = Max(-knz+iz0,0);
-   ize = Min(knz+iz0,fSeg->Npz()-1);
-   for (ix=ixs;ix<=ixe;ix++) 
-     for (iz=izs;iz<=ize;iz++) {
-       fSeg->DetToLocal(ix,iz,x,z); // pixel center
-       double dxi = 0.5*fSeg->Dpx(ix);
-       double dzi = 0.5*fSeg->Dpz(iz);
-       x1  = x;
-       z1  = z;
-       x2  = x1 + dxi; // Upper
-       x1 -= dxi;  // Lower
-       z2  = z1 + dzi; // Upper
-       z1 -= dzi;  // Lower
-       x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-       x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-       z1 -= z0; // Distance from where track traveled
-       z2 -= z0; // Distance from where track traveled
-       s   = el*0.25; // Correction based on definision of Erfc
-       s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
-       s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
-       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
-     } // end for ix, iz
-   //
-}
-
-//______________________________________________________________________
-void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
-                                           Int_t ix0,Int_t iz0,
-                                           Double_t el,Double_t sigx,Double_t sigz,
-                                           Double_t ld,Int_t t,Int_t hi)
+void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
+                                         Double_t el, Int_t tID, Int_t hID)
 {
   // Spreads the charge over neighboring cells. Assume charge is distributed
   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
-  // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
   // Defined this way, the integral over all x and z is el.
   // Inputs:
-  //    Double_t x0   x position of point where charge is liberated
-  //    Double_t z0   z position of point where charge is liberated
+  //    Double_t x0   x position of point where charge is liberated (local)
+  //    Double_t z0   z position of point where charge is liberated (local)
+  //    Double_t dy   distance from the entrance surface (diffusion sigma may depend on it)
   //    Int_t    ix0  row of cell corresponding to point x0
   //    Int_t    iz0  columb of cell corresponding to point z0
   //    Double_t el   number of electrons liberated in this step
   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
-  //    Double_t ld   lorentz drift in x for this stip (y0 dependent)
-  //    Int_t    t    track number
-  //    Int_t    ti   hit track index number
-  //    Int_t    hi   hit "hit" index number
-  // Outputs:
-  //     none.
-  // Return:
-  //     none.
-  const Int_t knx = 3,knz = 3; // RS: TO TUNE
-  const Double_t kRoot2 = 1.414213562; // Sqrt(2).
+  //    Int_t    tID  track number
+  //    Int_t    hID  hit "hit" index number
+  //
   Int_t ix,iz,ixs,ixe,izs,ize;
   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
-  Double_t s,spx,spz,x1,x2,z1,z2; 
+  Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
   //
-  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)",
-                               x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
-  if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
-    UpdateMapSignal(iz0,ix0,t,hi,el);
-    return;
-  } // end if
-  spx = 1.0/(sigx*kRoot2);     
-  spz = 1.0/(sigz*kRoot2);
-  ixs = Max(-knx+ix0,0);
-  ixe = Min(knx+ix0,fSeg->Npx()-1);
-  izs = Max(-knz+iz0,0);
-  ize = Min(knz+iz0,fSeg->Npz()-1);
+  if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,dy=%e, ix0=%d,iz0=%d,el=%e,tID=%d,hID=%d)",
+                               x0,z0,dy,ix0,iz0,el,tID,hID));
+  //
+  Double_t &x1 = dtIn[kCellX1];
+  Double_t &x2 = dtIn[kCellX2];
+  Double_t &z1 = dtIn[kCellZ1];
+  Double_t &z2 = dtIn[kCellZ2];
+  //
+  int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
+  int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
+  //
+  dtIn[kCellYDepth]  = dy;
+  ixs = Max(-nx+ix0,0);
+  ixe = Min( nx+ix0,fSeg->Npx()-1);
+  izs = Max(-nz+iz0,0);
+  ize = Min( nz+iz0,fSeg->Npz()-1);
   for (ix=ixs;ix<=ixe;ix++) 
     for (iz=izs;iz<=ize;iz++) {
       fSeg->DetToLocal(ix,iz,x,z); // pixel center
       double dxi = 0.5*fSeg->Dpx(ix);
       double dzi = 0.5*fSeg->Dpz(iz);
-      x1  = x;
-      z1  = z;
+      x1  = x - x0;   // calculate distance of cell boundaries from injection center
+      z1  = z - z0;
       x2  = x1 + dxi; // Upper
       x1 -= dxi;      // Lower
       z2  = z1 + dzi; // Upper
       z1 -= dzi;      // Lower
-      x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-      x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-      z1 -= z0; // Distance from where track traveled
-      z2 -= z0; // Distance from where track traveled
-      s   = el*0.25; // Correction based on definision of Erfc
-      s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
-      s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
-      if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
+      s   = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
+      //(*((AliITSUSimulationPix*)this)->AliITSUSimulationPix::fSpreadFun)(dtIn)); // calculate charge deposited in the cell
+      if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s);
     } // end for ix, iz
   //
 }
 
 //______________________________________________________________________
+Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
+{
+  // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2] 
+  // and Z=dtIn[kCellZ1]:dtIn[kCellZ2] 
+  // The spread function is assumed to be double gaussian in 2D
+  // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
+  //
+  int ip = kParamStart;
+  // 1st gaussian
+  double intg1 = GausInt2D(fResponseParam->GetParameter(ip+1),  // sigX
+                          fResponseParam->GetParameter(ip+3),  // sigZ
+                          dtIn[kCellX1]-fResponseParam->GetParameter(ip),      // x1-xmean
+                          dtIn[kCellX2]-fResponseParam->GetParameter(ip),      // x2-xmean
+                          dtIn[kCellZ1]-fResponseParam->GetParameter(ip+2),    // z1-zmean
+                          dtIn[kCellZ2]-fResponseParam->GetParameter(ip+2));   // z2-zmean
+  // 2nd gaussian
+  double intg2 = GausInt2D(fResponseParam->GetParameter(ip+5),  // sigX
+                          fResponseParam->GetParameter(ip+7),  // sigZ
+                          dtIn[kCellX1]-fResponseParam->GetParameter(ip+4),    // x1-xmean
+                          dtIn[kCellX2]-fResponseParam->GetParameter(ip+4),    // x2-xmean
+                          dtIn[kCellZ1]-fResponseParam->GetParameter(ip+6),    // z1-zmean
+                          dtIn[kCellZ2]-fResponseParam->GetParameter(ip+6));   // z2-zmean
+  double scl = fResponseParam->GetParameter(ip+8);
+  return (intg1+intg2*scl)/(1+scl);
+  //
+} 
+
+//______________________________________________________________________
+Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
+{
+  // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2] 
+  // and Z=dtIn[kCellZ1]:dtIn[kCellZ2] 
+  // The spread function is assumed to be gaussian in 2D
+  // Parameters should be: mean0,sigma0
+  int ip = kParamStart;
+  return GausInt2D(fResponseParam->GetParameter(ip+1),  // sigX
+                  fResponseParam->GetParameter(ip+3),  // sigZ
+                  dtIn[kCellX1]-fResponseParam->GetParameter(ip),      // x1-xmean
+                  dtIn[kCellX2]-fResponseParam->GetParameter(ip),      // x2-xmean
+                  dtIn[kCellZ1]-fResponseParam->GetParameter(ip+2),    // z1-zmean
+                  dtIn[kCellZ2]-fResponseParam->GetParameter(ip+2));
+  //
+} 
+
+//______________________________________________________________________
 void AliITSUSimulationPix::RemoveDeadPixels() 
 {
   // Removes dead pixels on each module (ladder)
@@ -564,7 +511,7 @@ void AliITSUSimulationPix::AddNoisyPixels()
   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
   if (!calObj) return;
   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), 
-                                                       10*fSimuParam->GetPixThreshold(fModule));
+                                                       10*fSimuParam->GetPixThreshold(fModule->GetIndex()));
   //
 }
 
@@ -590,13 +537,14 @@ void AliITSUSimulationPix::FrompListToDigits()
   //
   int maxInd = fSensMap->GetMaxIndex();
   double minProb = 0.1/maxInd;
+  int modId = fModule->GetIndex();
   //
   int nsd = fSensMap->GetEntries();
   Int_t prevID=0,curID=0;
   TArrayI ordSampleInd(100),ordSample(100);
   //
-  double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
-  fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
+  double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(modId);
+  fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
   //
   for (int i=0;i<nsd;i++) {
@@ -609,7 +557,7 @@ void AliITSUSimulationPix::FrompListToDigits()
       prevID = curID+1;
     }
     //
-    if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
+    if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(modId)) continue;
     if (Abs(sig)>2147483647.0) { //RS?
       //PH 2147483647 is the max. integer
       //PH This apparently is a problem which needs investigation
@@ -745,6 +693,7 @@ void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t
   // module         module number
   //
   UInt_t col,row;
+  Int_t modId = fModule->GetIndex();
   Double_t pulse1,pulse2;
   Double_t couplR=0.0,couplC=0.0;
   //
@@ -759,13 +708,13 @@ void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t
    //
    int j1 = int(col)+isign;
    pulse1 *= couplC;    
-   if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
+   if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
    
    // loop in row direction
    int j2 = int(row) + isign;
    pulse2 *= couplR;
-   if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
+   if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
  } // for isign
 }
@@ -777,8 +726,23 @@ void AliITSUSimulationPix::GenerateStrobePhase()
   // phase w.r.t to the LHC clock
   // Done once per event
   const Double_t kBunchLenght = 25e-9; // LHC clock
+  if (!fStrobe) return;
   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
     (Double_t)fStrobeLenght*kBunchLenght+
     kBunchLenght/2;
 }
 
+//______________________________________________________________________
+void AliITSUSimulationPix::SetResponseParam(AliParamList* resp)
+{
+  // attach response parameterisation data
+  fResponseParam = resp;
+  switch (fResponseParam->GetID()) {
+  case kSpreadSingleGauss: fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D; 
+    break;
+  case kSpreadDoubleGauss: fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;       
+    break;
+  default: AliFatal(Form("Did not find requested spread function id=%d",fResponseParam->GetID()));
+  }
+  //
+}
index b03d72b..0462c28 100644 (file)
 class TH1F;
 class AliITSUModule;
 class AliITSUSimuParam;
+class AliParamList;
 
 //-------------------------------------------------------------------
 
 class AliITSUSimulationPix : public AliITSUSimulation {
 public:
+  enum {kCellX1,kCellX2,kCellZ1,kCellZ2,kCellYDepth,kNDtSpread}; // data used for ch. spread integral calc.
+  //
+  // charge spread functions defined
+  enum {kSpreadSingleGauss,                  // single gaussian in 2D, SpreadFunGauss2D
+       kSpreadDoubleGauss,      // double gaussian in 2D, SpreadFunDoubleGauss2D
+       kNSpreadFuns
+  };
+  // fist kParamStart entried of spread fun params are reserved for common parameters
+  enum {kSpreadFunParamNXoffs,               // number of pixels to consider +- from injection point (in X)
+       kSpreadFunParamNZoffs,               // number of pixels to consider +- from injection point (in Z)
+       kParamStart
+  };
+  //
   AliITSUSimulationPix();
   AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map);
   virtual ~AliITSUSimulationPix();
@@ -28,12 +42,12 @@ public:
   void Init();
   //
   void FinishSDigitiseModule();
-  void DigitiseModule(AliITSUModule *mod,Int_t mask, Int_t event, AliITSsegmentation* seg);
+  void DigitiseModule();
   //
-  void SDigitiseModule(AliITSUModule *mod, Int_t mask, Int_t event, AliITSsegmentation* seg);
+  void SDigitiseModule();
   void WriteSDigits();
-  void Hits2SDigits(AliITSUModule *mod);
-  void Hits2SDigitsFast(AliITSUModule *mod);
+  void Hits2SDigits();
+  void Hits2SDigitsFast();
   void AddNoisyPixels();   
   void RemoveDeadPixels();
   void FrompListToDigits();
@@ -46,12 +60,16 @@ public:
   
   // This sets fStrobe flag and allows generating the strobe and applying it to select hits 
   void SetStrobeGeneration(Bool_t b=kFALSE) {fStrobe=b;};
-  void GenerateStrobePhase();
+  virtual void GenerateStrobePhase();
+  //
+  Double_t SpreadFunDoubleGauss2D(const Double_t *dtIn);
+  Double_t SpreadFunGauss2D(const Double_t *dtIn);
+  //
+  virtual void SetResponseParam(AliParamList* resp);
   //
-  
  private:
-  void SpreadCharge(Double_t x0,Double_t z0,Int_t ix0,Int_t iz0,Double_t el,Double_t sig,Double_t ld,Int_t t,Int_t hi);
-  void SpreadChargeAsym(Double_t x0,Double_t z0,Int_t ix0,Int_t iz0,Double_t el,Double_t sigx,Double_t sigz,Double_t ld,Int_t t,Int_t hi);
+  void SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
+                     Double_t el, Int_t tID, Int_t hID);
   //
   void SetCoupling(AliITSUSDigit* old,Int_t ntrack,Int_t idhit);     // "New" coupling routine  Tiziano Virgili
   void SetCouplingOld(AliITSUSDigit* old,Int_t ntrack,Int_t idhit);  // "Old" coupling routine  Rocco Caliandro
@@ -61,6 +79,9 @@ public:
    Bool_t        fStrobe;       // kTRUE if readout strobe with proper phase applied to select hits
    Int_t         fStrobeLenght; // Strobe signal lenght in units of 25 ns
    Double_t      fStrobePhase;  // The phase of the strobe signal with respect to the trigger
+   //   
+   Double_t (AliITSUSimulationPix::*fSpreadFun)(const Double_t *dtIn); //! pointer on current spread function
+
    ClassDef(AliITSUSimulationPix,1)  // Simulation of pixel clusters
  };
 #endif 
index edf12a8..16597b3 100644 (file)
@@ -130,7 +130,7 @@ AliITSUv11::AliITSUv11(const char *title,const Int_t nlay)
       fModPerLadd[j] = 0;
       fLadWidth[j] = 0.;
       fDetThick[j] = 0.;
-      fDetTypeID[j] = 0.;
+      fDetTypeID[j] = 0;
       fUpGeom[j] = 0;
     }
   }
diff --git a/ITS/UPGRADE/testITSU/MakeITSUSimuParam.C b/ITS/UPGRADE/testITSU/MakeITSUSimuParam.C
new file mode 100644 (file)
index 0000000..1dd1395
--- /dev/null
@@ -0,0 +1,115 @@
+void MakeITSUSimuParam(const char* cdbURI="local://") {
+//========================================================================
+//
+// Steering macro for ITS simulation parameters
+//
+// Author: L.Molnar
+// Contact: levente.molnar@cern.ch
+//
+//========================================================================
+
+  const char* macroname = "MakeITSUSimuParam.C";
+  //
+  gSystem->Load("libITSUpgradeBase.so");
+  gSystem->Load("libITSUpgradeSim.so");
+  gSystem->Load("libITSUpgradeRec.so");
+  //
+  // Activate CDB storage and load geometry from CDB
+  AliCDBManager* cdb = AliCDBManager::Instance();
+  cdb->SetDefaultStorage(cdbURI);
+
+  AliITSUSimuParam* itsSimuParam = new AliITSUSimuParam();
+  //
+  // Add spread function parameterization data
+  AliParamList* parData = 0;
+  int offs = 0;
+  //
+  //------------------------ parameterization data for segmentation 0 ----------------------
+  parData = new AliParamList(AliITSUSimulationPix::kParamStart+11); // 2 common + 9 params for double gaussian
+  parData->SetUniqueID(0); // this is a function for detId=0
+  parData->SetID(AliITSUSimulationPix::kSpreadDoubleGauss); // and uses double gaussian
+  parData->SetNameTitle("Monopix_seg0","double gaussian for segmentation 0");
+  //
+  // obligatory params for all AliITSUSimulationPix functions: number of pixels in X,Z around
+  // injected one to consider
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNXoffs,3,"nPixX");
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNZoffs,3,"nPixZ"); 
+  // 
+  // now set the parameters according selected function
+  offs = AliITSUSimulationPix::kParamStart;
+  parData->SetParameter(offs++,-0.1e-4  , "G1 Mean_x");
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_x");
+  parData->SetParameter(offs++, 2.011e-4, "G1 Mean_z"); 
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_z"); 
+  parData->SetParameter(offs++,-0.069e-4, "G2 Mean_x");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_x");
+  parData->SetParameter(offs++,-8.425e-4, "G2 Mean_z");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_z"); 
+  parData->SetParameter(offs++,3.904037/59.468672, "G2 A2/A1");  // scaling of 2nd gaussian amplitude wrt 1st one
+  // 
+  itsSimuParam->AddRespFunParam(parData);
+  //
+  //
+  //------------------------ parameterization data for segmentation 1 ----------------------
+  parData = new AliParamList(AliITSUSimulationPix::kParamStart+9); // 2 common + 9 params for double gaussian
+  parData->SetUniqueID(1); // this is a function for detId=1
+  parData->SetID(AliITSUSimulationPix::kSpreadDoubleGauss); // and uses double gaussian
+  parData->SetNameTitle("Monopix_seg1","double gaussian for segmentation 1");
+  //
+  // obligatory params for all AliITSUSimulationPix functions: number of pixels in X,Z around
+  // injected one to consider
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNXoffs,3,"nPixX");
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNZoffs,3,"nPixZ"); 
+  // 
+  // now set the parameters according selected function
+  offs = AliITSUSimulationPix::kParamStart;
+  parData->SetParameter(offs++,-0.1e-4  , "G1 Mean_x");
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_x");
+  parData->SetParameter(offs++, 2.011e-4, "G1 Mean_z"); 
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_z"); 
+  parData->SetParameter(offs++,-0.069e-4, "G2 Mean_x");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_x");
+  parData->SetParameter(offs++,-8.425e-4, "G2 Mean_z");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_z"); 
+  parData->SetParameter(offs++,3.904037/59.468672, "G2 A2/A1");  // scaling of 2nd gaussian amplitude wrt 1st one
+  // 
+  itsSimuParam->AddRespFunParam(parData);
+  //
+  //
+  //------------------------ parameterization data for segmentation 2 ----------------------
+  parData = new AliParamList(AliITSUSimulationPix::kParamStart+9); // 2 common + 9 params for double gaussian
+  parData->SetUniqueID(2); // this is a function for detId=2
+  parData->SetID(AliITSUSimulationPix::kSpreadDoubleGauss); // and uses double gaussian
+  parData->SetNameTitle("Monopix_seg1","double gaussian for segmentation 1");
+  //
+  // obligatory params for all AliITSUSimulationPix functions: number of pixels in X,Z around
+  // injected one to consider
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNXoffs,3,"nPixX");
+  parData->SetParameter(AliITSUSimulationPix::kSpreadFunParamNZoffs,3,"nPixZ"); 
+  // 
+  // now set the parameters according selected function
+  offs = AliITSUSimulationPix::kParamStart;
+  parData->SetParameter(offs++,-0.1e-4  , "G1 Mean_x");
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_x");
+  parData->SetParameter(offs++, 2.011e-4, "G1 Mean_z"); 
+  parData->SetParameter(offs++, 8.125e-4, "G1 Sigma_z"); 
+  parData->SetParameter(offs++,-0.069e-4, "G2 Mean_x");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_x");
+  parData->SetParameter(offs++,-8.425e-4, "G2 Mean_z");
+  parData->SetParameter(offs++,15.050e-4, "G2 Sigma_z"); 
+  parData->SetParameter(offs++,3.904037/59.468672, "G2 A2/A1");  // scaling of 2nd gaussian amplitude wrt 1st one
+  // 
+  itsSimuParam->AddRespFunParam(parData);
+  //
+  // save in CDB storage
+  AliCDBMetaData *md= new AliCDBMetaData();
+  md->SetResponsible("ITS Upgrade Project");
+  md->SetComment("Simulation parameters for ITS Upgrade.");
+  md->SetAliRootVersion(gSystem->Getenv("ARVERSION"));
+  md->SetBeamPeriod(0);
+  AliCDBId id("ITS/Calib/SimuParam",0,AliCDBRunRange::Infinity());
+  cdb->GetDefaultStorage()->Put(itsSimuParam,id, md);
+  //
+  return;
+}
+
index 63f94c4..7f930ec 100644 (file)
@@ -11,6 +11,8 @@ void sim(Int_t nev=3) {
                               Form("local://%s",gSystem->pwd()));
   simulator.SetSpecificStorage("ITS/Align/Data",
                               Form("local://%s",gSystem->pwd()));
+  simulator.SetSpecificStorage("ITS/Calib/SimuParam",
+                              Form("local://%s",gSystem->pwd()));
   simulator.SetRunHLT("");
   simulator.SetRunQA(":");