]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Fixing Coding violations
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index f006ca43ecade53a0b7cea94787c831121e7f747..bda9cdd77f0ebc7ca4c7a0c66eeab31f09f74bb1 100644 (file)
@@ -51,6 +51,8 @@
 #include <TTree.h>
 #include <TVector.h>
 #include <TVirtualMC.h>
+#include <TString.h>
+#include <TF2.h>
 
 #include "AliArrayBranch.h"
 #include "AliClusters.h"
 #include "AliTPCTrackHits.h"
 #include "AliTPCTrackHitsV2.h"
 #include "AliTPCcluster.h"
-#include "AliTPCclusterer.h"
-#include "AliTPCtracker.h"
 #include "AliTrackReference.h"
+#include "AliMC.h"
+#include "AliTPCDigitizer.h"
+#include "AliTPCclustererMI.h"
+#include "AliTPCtrackerMI.h"
+#include "AliTPCpidESD.h"
 
 
 ClassImp(AliTPC) 
@@ -86,13 +91,18 @@ ClassImp(AliTPC)
 
 class AliTPCFastMatrix : public TMatrix {
 public :
-  AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb);
-  inline Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
-  inline Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+  AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb);
+#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
+  Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
+  Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
+#else
+  Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+  Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+#endif
 };
 
-AliTPCFastMatrix::AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb):
-  TMatrix(row_lwb, row_upb,col_lwb,col_upb)
+AliTPCFastMatrix::AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb):
+  TMatrix(rowlwb, rowupb,collwb,colupb)
    {
    };
 //_____________________________________________________________________________
@@ -100,7 +110,7 @@ class AliTPCFastVector : public TVector {
 public :
   AliTPCFastVector(Int_t size);
   virtual ~AliTPCFastVector(){;}
-  inline Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces  
+  Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces  
   
 };
 
@@ -140,7 +150,7 @@ AliTPC::AliTPC(const char *name, const char *title)
   //
   // Initialise arrays of hits and digits 
   fHits     = new TClonesArray("AliTPChit",  176);
-  gAlice->AddHitList(fHits); 
+  gAlice->GetMCApp()->AddHitList(fHits); 
   fDigitsArray = 0;
   fClustersArray= 0;
   fDefaults = 0;
@@ -184,6 +194,11 @@ AliTPC::AliTPC(const char *name, const char *title)
 }
 
 //_____________________________________________________________________________
+AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
+  //
+  // dummy copy constructor
+  //
+}
 AliTPC::~AliTPC()
 {
   //
@@ -302,7 +317,7 @@ void AliTPC::BuildGeometry()
 }    
 
 //_____________________________________________________________________________
-Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
+Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const
 {
   //
   // Calculate distance from TPC to mouse on the display
@@ -311,15 +326,16 @@ Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
   return 9999;
 }
 
-void AliTPC::Clusters2Tracks() 
+void AliTPC::Clusters2Tracks() const 
  {
   //-----------------------------------------------------------------
   // This is a track finder.
   //-----------------------------------------------------------------
-  AliTPCtracker tracker(fTPCParam,0,fLoader->GetEventFolder()->GetName());
-  tracker.Clusters2Tracks();
+  Error("Clusters2Tracks",
+  "Dummy function !  Call AliTPCtracker::Clusters2Tracks(...) instead !");
  }
 
+
 //_____________________________________________________________________________
 void AliTPC::CreateMaterials()
 {
@@ -432,11 +448,11 @@ void AliTPC::CreateMaterials()
   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
   //----------------------------------------------------------------
 
-  char namate[21]; 
+  char namate[21]=""
   density = 0.;
   Float_t am=0;
   Int_t nc;
-  Float_t rho,absl,X0,buf[1];
+  Float_t rho,absl,x0,buf[1];
   Int_t nbuf;
   Float_t a,z;
 
@@ -445,7 +461,7 @@ void AliTPC::CreateMaterials()
     
       // retrive material constants
       
-      gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
+      gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
 
       amat[nc] = a;
       zmat[nc] = z;
@@ -675,7 +691,7 @@ void AliTPC::CreateMaterials()
 
   // get epoxy
 
-  gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf);
+  gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf);
 
   // Carbon fiber
 
@@ -688,7 +704,7 @@ void AliTPC::CreateMaterials()
 
   // get SiO2
 
-  gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); 
+  gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf); 
 
   wmat[0]=0.725; // by weight!
   wmat[1]=0.275;
@@ -758,7 +774,7 @@ Float_t AliTPC::GetNoise()
 }
 
 
-Bool_t  AliTPC::IsSectorActive(Int_t sec)
+Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 {
   //
   // check if the sector is active
@@ -802,7 +818,7 @@ void    AliTPC::SetActiveSectors(Int_t flag)
   else branch = TreeH()->GetBranch("TPC");
   Stat_t ntracks = TreeH()->GetEntries();
   // loop over all hits
-  cout<<"\nAliTPC::SetActiveSectors():  Got "<<ntracks<<" tracks\n";
+  if (GetDebug()) cout<<"\nAliTPC::SetActiveSectors():  Got "<<ntracks<<" tracks\n";
   
   for(Int_t track=0;track<ntracks;track++)
    {
@@ -828,24 +844,24 @@ void    AliTPC::SetActiveSectors(Int_t flag)
       } 
     }    
   }
-  
 }  
 
 
 
 
-void AliTPC::Digits2Clusters(Int_t eventnumber)
+void AliTPC::Digits2Clusters(Int_t /*eventnumber*/) const
 {
   //-----------------------------------------------------------------
   // This is a simple cluster finder.
   //-----------------------------------------------------------------
-  AliTPCclusterer::Digits2Clusters(fTPCParam, fLoader,eventnumber);
+  Error("Digits2Clusters",
+  "Dummy function !  Call AliTPCclusterer::Digits2Clusters(...) instead !");
 }
 
 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
 extern Double_t SigmaZ2(Double_t, Double_t);
 //_____________________________________________________________________________
-void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
+void AliTPC::Hits2Clusters(Int_t /*eventn*/)
 {
   //--------------------------------------------------------
   // TPC simple cluster generator from hits
@@ -915,9 +931,9 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
   
   cout<<"fTPCParam->GetTitle() = "<<fTPCParam->GetTitle()<<endl;
   
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
+  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->CdGAFile();
-  fTPCParam->Write(fTPCParam->GetTitle());
+  //fTPCParam->Write(fTPCParam->GetTitle());
 
   AliTPCClustersArray carray;
   carray.Setup(fTPCParam);
@@ -966,7 +982,7 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
         continue; 
        }
        ipart=tpcHit->Track();
-       particle=gAlice->Particle(ipart);
+       particle=gAlice->GetMCApp()->Particle(ipart);
        pl=particle->Pz();
        pt=particle->Pt();
        if(pt < 1.e-9) pt=1.e-9;
@@ -1036,7 +1052,7 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
 
   } // end of loop over sectors  
 
-  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
+  //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
   fLoader->WriteRecPoints("OVERWRITE");
   
   
@@ -1083,7 +1099,7 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
   SetTreeAddress();
 
   Stat_t ntracks = tH->GetEntries();
-  Int_t npart = gAlice->GetNtrack();
+  Int_t npart = gAlice->GetMCApp()->GetNtrack();
   //MI change
   TBranch * branch=0;
   if (fHitType>1) branch = tH->GetBranch("TPC2");
@@ -1120,7 +1136,7 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
       }
 
       ipart=tpcHit->Track();
-      if (ipart<npart) particle=gAlice->Particle(ipart);
+      if (ipart<npart) particle=gAlice->GetMCApp()->Particle(ipart);
       
       //find row number
 
@@ -1227,8 +1243,13 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
 
 
 
+//______________________________________________________________________
+AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
+{
+  return new AliTPCDigitizer(manager);
+}
 //__
-void AliTPC::SDigits2Digits2(Int_t eventnumber)  
+void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)  
 {
   //create digits from summable digits
   GenerNoise(500000); //create teble with noise
@@ -1351,14 +1372,16 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
 }
 //__________________________________________________________________
 void AliTPC::SetDefaults(){
-
+  //
+  // setting the defaults
+  //
    
-   cerr<<"Setting default parameters...\n";
+  //   cerr<<"Setting default parameters...\n";
 
   // Set response functions
 
   //
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
+  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->CdGAFile();
   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
   if(param){
@@ -1389,9 +1412,24 @@ void AliTPC::SetDefaults(){
     cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
      exit(3);
   }
+
+  TString s;
   prfinner->Read("prf_07504_Gati_056068_d02");
+  //PH Set different names
+  s=prfinner->GetGRF()->GetName();
+  s+="in";
+  prfinner->GetGRF()->SetName(s.Data());
+
   prfouter1->Read("prf_10006_Gati_047051_d03");
+  s=prfouter1->GetGRF()->GetName();
+  s+="out1";
+  prfouter1->GetGRF()->SetName(s.Data());
+
   prfouter2->Read("prf_15006_Gati_047051_d03");  
+  s=prfouter2->GetGRF()->GetName();
+  s+="out2";
+  prfouter2->GetGRF()->SetName(s.Data());
+
   f->Close();
   savedir->cd();
 
@@ -1409,12 +1447,32 @@ void AliTPC::SetDefaults(){
 
 }
 //__________________________________________________________________  
+void AliTPC::Hits2Digits()  
+{
+  //
+  // creates digits from hits
+  //
+
+  fLoader->LoadHits("read");
+  fLoader->LoadDigits("recreate");
+  AliRunLoader* runLoader = fLoader->GetRunLoader(); 
+
+  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+    runLoader->GetEvent(iEvent);
+    SetActiveSectors();   
+    Hits2Digits(iEvent);
+  }
+
+  fLoader->UnloadHits();
+  fLoader->UnloadDigits();
+} 
+//__________________________________________________________________  
 void AliTPC::Hits2Digits(Int_t eventnumber)  
 { 
  //----------------------------------------------------
  // Loop over all sectors for a single event
  //----------------------------------------------------
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
+  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->GetEvent(eventnumber);
   if (fLoader->TreeH() == 0x0)
    {
@@ -1451,24 +1509,30 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
 
   fDigitsSwitch=0; // standard digits
 
-  cerr<<"Digitizing TPC -- normal digits...\n";
+  //  cerr<<"Digitizing TPC -- normal digits...\n";
 
  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) 
   if (IsSectorActive(isec)) 
    {
-    cout<<"Sector "<<isec<<"is active\n";
+    if (fDebug) Info("Hits2Digits","Sector %d is active.",isec);
     Hits2DigitsSector(isec);
    }
   else
    {
-    cout<<"Sector "<<isec<<"is NOT active\n";
+    if (fDebug) Info("Hits2Digits","Sector %d is NOT active.",isec);
    }
 
   fLoader->WriteDigits("OVERWRITE"); 
+  
+//this line prevents the crash in the similar one
+//on the beginning of this method
+//destructor attempts to reset the tree, which is deleted by the loader
+//need to be redesign
+ if(GetDigitsArray()) delete GetDigitsArray();
+ SetDigitsArray(0x0);
+  
 }
 
-
-
 //__________________________________________________________________
 void AliTPC::Hits2SDigits2(Int_t eventnumber)  
 { 
@@ -1516,7 +1580,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
 
   SetDigitsArray(arr);
 
-  cerr<<"Digitizing TPC -- summable digits...\n"; 
+  //  cerr<<"Digitizing TPC -- summable digits...\n"; 
 
   fDigitsSwitch=1; // summable digits
   
@@ -1549,53 +1613,19 @@ void AliTPC::Hits2SDigits()
   //   summable digits - 16 bit "ADC", no noise, no saturation
   //-----------------------------------------------------------
 
- //----------------------------------------------------
- // Loop over all sectors for a single event
- //----------------------------------------------------
-  //MI change - for pp run
-//  Int_t eventnumber = gAlice->GetEvNumber();
-//  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
-//  rl->GetEvent(eventnumber);
-// 12/05/2003 This method is obsolete and not used. It should be redesingned
-// M.Kowalski
-
-  if (fLoader->TreeH() == 0x0)
-   {
-     if(fLoader->LoadHits())
-      {
-        Error("Hits2Digits","Can not load hits.");
-      }
-   }
-  SetTreeAddress();
-
-  if(fDefaults == 0) SetDefaults();
-  GenerNoise(500000); //create table with noise
-
-  //setup TPCDigitsArray 
+  fLoader->LoadHits("read");
+  fLoader->LoadSDigits("recreate");
+  AliRunLoader* runLoader = fLoader->GetRunLoader(); 
 
-  if(GetDigitsArray()) delete GetDigitsArray();
-
-  if (fLoader->TreeS() == 0x0 ) 
-   {
-     fLoader->MakeTree("S");
-   }
-  
-  AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
-  arr->SetClass("AliSimDigits");
-  arr->Setup(fTPCParam);
-  arr->MakeTree(fLoader->TreeS());
-  SetDigitsArray(arr);
-
-  cerr<<"Digitizing TPC -- summable digits...\n"; 
-
-  //  fDigitsSwitch=1; // summable digits  -for the moment direct
-
-  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
+  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+    runLoader->GetEvent(iEvent);
+    SetTreeAddress();
+    SetActiveSectors();
+    Hits2SDigits2(iEvent);
+  }
 
-  // write results
-  //
-  cout<<"Why method TPC::Hits2SDigits writes digits and not sdigits? skowron\n";
-  fLoader->WriteDigits("OVERWRITE");
+  fLoader->UnloadHits();
+  fLoader->UnloadSDigits();
 }
 //_____________________________________________________________________________
 
@@ -2427,7 +2457,7 @@ AliHit(shunt,track)
 //________________________________________________________________________
 // Additional code because of the AliTPCTrackHitsV2
 
-void AliTPC::MakeBranch2(Option_t *option,const char *file)
+void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
 {
   //
   // Create a new branch in the current Root Tree
@@ -2485,18 +2515,21 @@ void AliTPC::SetTreeAddress2()
     if (branch) 
      {
        branch->SetAddress(&fTrackHits);
-       cout<<"AliTPC::SetTreeAddress2 fHitType&4 Setting"<<endl;       
+       if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Setting");
      }
-    else cout<<"AliTPC::SetTreeAddress2 fHitType&4  Failed"<<endl;
+    else 
+    if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Failed (can not find branch)");
+    
   }
   if ((treeH)&&(fHitType&2)) {
     branch = treeH->GetBranch(branchname);
     if (branch) 
      {
        branch->SetAddress(&fTrackHitsOld);
-       cout<<"AliTPC::SetTreeAddress2 fHitType&2 Setting"<<endl;       
+       if (GetDebug()) Info("SetTreeAddress2","fHitType&2 Setting");
      }
-    else cout<<"AliTPC::SetTreeAddress2 fHitType&2  Failed"<<endl;
+    else if (GetDebug()) 
+      Info("SetTreeAddress2","fHitType&2 Failed (can not find branch)");
   }
   //set address to TREETR
 
@@ -2521,12 +2554,12 @@ void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
   // add hit to the list  
   Int_t rtrack;
   if (fIshunt) {
-    int primary = gAlice->GetPrimary(track);
-    gAlice->Particle(primary)->SetBit(kKeepBit);
+    int primary = gAlice->GetMCApp()->GetPrimary(track);
+    gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
     rtrack=primary;
   } else {
     rtrack=track;
-    gAlice->FlagTrack(track);
+    gAlice->GetMCApp()->FlagTrack(track);
   }  
   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
   //if (hit->fTrack!=rtrack)
@@ -2562,6 +2595,9 @@ AliHit* AliTPC::FirstHit(Int_t track)
 }
 AliHit* AliTPC::NextHit()
 {
+  //
+  // gets next hit
+  //
   if (fHitType>1) return NextHit2();
   
   return AliDetector::NextHit();
@@ -2628,6 +2664,9 @@ void AliTPC::LoadPoints(Int_t)
 
 void AliTPC::RemapTrackHitIDs(Int_t *map)
 {
+  //
+  // remapping
+  //
   if (!fTrackHits) return;
   
   if (fTrackHitsOld && fHitType&2){
@@ -2701,7 +2740,7 @@ void AliTPC::LoadPoints2(Int_t)
   if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
   
   if (nhits == 0) return;
-  Int_t tracks = gAlice->GetNtrack();
+  Int_t tracks = gAlice->GetMCApp()->GetNtrack();
   if (fPoints == 0) fPoints = new TObjArray(tracks);
   AliHit *ahit;
   //
@@ -2776,7 +2815,7 @@ void AliTPC::LoadPoints3(Int_t)
   //
   Int_t nhits = fTrackHits->GetEntriesFast();
   if (nhits == 0) return;
-  Int_t tracks = gAlice->GetNtrack();
+  Int_t tracks = gAlice->GetMCApp()->GetNtrack();
   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
   fPoints->Expand(2*tracks);
   AliHit *ahit;
@@ -2854,7 +2893,7 @@ void AliTPC::LoadPoints3(Int_t)
 
 
 
-void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
+void AliTPC::FindTrackHitsIntersection(TClonesArray * /*arr*/)
 {
 
   //
@@ -3027,19 +3066,8 @@ void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
 
 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
 {
- cout<<"AliTPC::MakeLoader ";
+//Makes TPC loader
  fLoader = new AliTPCLoader(GetName(),topfoldername);
- if (fLoader)
-  {
-   cout<<"Success"<<endl;
-  }
- else
-  {
-   cout<<"Failure"<<endl;
-  }
-
  return fLoader;
 }
 
@@ -3089,3 +3117,47 @@ AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
 
 }
 
+
+//____________________________________________________________________________
+Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
+{
+  //
+  // Parametrised error of the cluster reconstruction (pad direction)
+  //
+  // Sigma rphi
+  const Float_t kArphi=0.41818e-2;
+  const Float_t kBrphi=0.17460e-4;
+  const Float_t kCrphi=0.30993e-2;
+  const Float_t kDrphi=0.41061e-3;
+
+  pt=TMath::Abs(pt)*1000.;
+  Double_t x=r/pt;
+  tgl=TMath::Abs(tgl);
+  Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
+  if (s<0.4e-3) s=0.4e-3;
+  s*=1.3; //Iouri Belikov
+
+  return s;
+}
+
+
+//____________________________________________________________________________
+Double_t SigmaZ2(Double_t r, Double_t tgl)
+{
+  //
+  // Parametrised error of the cluster reconstruction (drift direction)
+  //
+  // Sigma z
+  const Float_t kAz=0.39614e-2;
+  const Float_t kBz=0.22443e-4;
+  const Float_t kCz=0.51504e-1;
+
+
+  tgl=TMath::Abs(tgl);
+  Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
+  if (s<0.4e-3) s=0.4e-3;
+  s*=1.3; //Iouri Belikov
+
+  return s;
+}
+