latest changes from Marian
authorMikolaj Krzewicki <mikolaj.krzewicki@cern.ch>
Mon, 9 Dec 2013 11:18:30 +0000 (12:18 +0100)
committerMikolaj Krzewicki <mikolaj.krzewicki@cern.ch>
Mon, 9 Dec 2013 17:03:14 +0000 (18:03 +0100)
STAT/TStatToolkit.cxx
STAT/TStatToolkit.h
TPC/Calib/AliTPCPreprocessorOffline.cxx
TPC/Calib/AliTPCcalibCalib.cxx
TPC/Calib/AliTPCcalibGainMult.cxx
TPC/Calib/AliTPCcalibGainMult.h
TPC/Calib/AliTPCcalibSummary.cxx

index 0f2c984..1a856fb 100644 (file)
@@ -19,8 +19,8 @@
 // 
 // Subset of  matheamtical functions  not included in the TMath
 //
-
-///////////////////////////////////////////////////////////////////////////
+//
+/////////////////////////////////////////////////////////////////////////
 #include "TMath.h"
 #include "Riostream.h"
 #include "TH1F.h"
@@ -38,7 +38,6 @@
 #include "TCanvas.h"
 #include "TLatex.h"
 #include "TCut.h"
-
 //
 // includes neccessary for test functions
 //
@@ -298,6 +297,96 @@ void TStatToolkit::LTM(TH1 * his, TVectorD *param , Float_t fraction,  Bool_t ve
   }
 }
 
+
+void TStatToolkit::MedianFilter(TH1 * his1D, Int_t nmedian){
+  //
+  // Algorithm to filter  histogram
+  // author:  marian.ivanov@cern.ch
+  // Details of algorithm:
+  // http://en.wikipedia.org/w/index.php?title=Median_filter&oldid=582191524
+  // Input parameters:
+  //    his1D - input histogam - to be modiefied by Medianfilter
+  //    nmendian - number of bins in median filter
+  //
+  Int_t nbins    = his1D->GetNbinsX();
+  TVectorD vectorH(nbins);
+  for (Int_t ibin=0; ibin<nbins; ibin++) vectorH[ibin]=his1D->GetBinContent(ibin+1);
+  for (Int_t ibin=0; ibin<nbins; ibin++) {
+    Int_t index0=ibin-nmedian;
+    Int_t index1=ibin+nmedian;
+    if (index0<0) {index1+=-index0; index0=0;}
+    if (index1>=nbins) {index0-=index1-nbins+1; index1=nbins-1;}    
+    Double_t value= TMath::Median(index1-index0,&(vectorH.GetMatrixArray()[index0]));
+    his1D->SetBinContent(ibin+1, value);
+  }  
+}
+
+Bool_t TStatToolkit::LTMHisto(TH1 *his1D, TVectorD &params , Float_t fraction){
+  //
+  // LTM : Trimmed mean on histogram - Modified version for binned data
+  // Robust statistic to estimate properties of the distribution
+  // See http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
+  //
+  // Paramters:
+  //     his1D   - input histogram
+  //     params  - vector with parameters
+  //             - 0 - area
+  //             - 1 - mean
+  //             - 2 - rms
+  //             - 3 - error estimate of mean
+  //             - 4 - dummy
+  //             - 5 - first accepted bin position
+  //             - 6 - last accepted  bin position
+  //
+  Int_t nbins    = his1D->GetNbinsX();
+  Int_t nentries = (Int_t)his1D->GetEntries();
+  if (nentries<=0) return 0;
+  TVectorD vectorX(nbins);
+  TVectorD vectorMean(nbins);
+  TVectorD vectorRMS(nbins);
+  Double_t sumCont=0;
+  for (Int_t ibin0=1; ibin0<nbins; ibin0++) sumCont+=his1D->GetBinContent(ibin0);
+  //
+  Double_t minRMS=his1D->GetRMS()*10000;
+  Int_t maxBin=0;
+  //
+  for (Int_t ibin0=1; ibin0<nbins; ibin0++){
+    Double_t sum0=0, sum1=0, sum2=0;
+    Int_t ibin1=ibin0;
+    for ( ibin1=ibin0; ibin1<nbins; ibin1++){
+      Double_t cont=his1D->GetBinContent(ibin1);
+      Double_t x= his1D->GetBinCenter(ibin1);
+      sum0+=cont;
+      sum1+=cont*x;
+      sum2+=cont*x*x;
+      if (sum0>fraction*sumCont) break;
+    }
+    vectorX[ibin0]=his1D->GetBinCenter(ibin0);
+    if (sum0>fraction*sumCont){
+      vectorMean[ibin0]=sum1/sum0;
+      vectorRMS[ibin0]=TMath::Sqrt(sum2/sum0-vectorMean[ibin0]*vectorMean[ibin0]);
+      if (vectorRMS[ibin0]<minRMS){
+       minRMS=vectorRMS[ibin0];
+       params[0]=sum0;
+       params[1]=vectorMean[ibin0];
+       params[2]=vectorRMS[ibin0];
+       params[3]=vectorRMS[ibin0]/TMath::Sqrt(sumCont*fraction);
+       params[4]=0;
+       params[5]=ibin0;
+       params[6]=ibin1;
+       params[7]=his1D->GetBinCenter(ibin0);
+       params[8]=his1D->GetBinCenter(ibin1);
+       maxBin=ibin0;
+      }
+    }else{
+      break;
+    }
+  }
+  return kTRUE;
+
+}
+
+
 Double_t  TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
   //
   //  Fit histogram with gaussian function
@@ -691,7 +780,7 @@ TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t
       if (type==0) stat = projection->GetMean();
       if (type==1) stat = projection->GetRMS();
       if (type==2 || type==3){
-       TVectorD vec(3);
+       TVectorD vec(10);
        TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
        if (type==2) stat= vec[1];
        if (type==3) stat= vec[0];      
@@ -736,15 +825,16 @@ TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t frac
   TVectorD vecYErr(nbinx);
   //
   TF1 f1("f1","gaus");
-  TVectorD vecLTM(3);
+  TVectorD vecLTM(10);
 
-  for (Int_t ix=0; ix<nbinx;ix++){
-    Float_t xcenter = xaxis->GetBinCenter(ix); 
+  for (Int_t jx=1; jx<=nbinx;jx++){
+    Int_t ix=jx-1;
+    Float_t xcenter = xaxis->GetBinCenter(jx); 
     snprintf(name,1000,"%s_%d",his->GetName(), ix);
-    TH1 *projection = his->ProjectionY(name,TMath::Max(ix-deltaBin,1),TMath::Min(ix+deltaBin,nbinx));
+    TH1 *projection = his->ProjectionY(name,TMath::Max(jx-deltaBin,1),TMath::Min(jx+deltaBin,nbinx));
     Double_t stat= 0;
     Double_t err =0;
-    TStatToolkit::LTM((TH1F*)projection,&vecLTM,fraction);  
+    TStatToolkit::LTMHisto((TH1F*)projection,vecLTM,fraction);  
     //
     if (returnType==0) {
       stat = projection->GetMean();
@@ -759,7 +849,7 @@ TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t frac
       if (returnType==3) stat= vecLTM[2];      
     }
     if (returnType==4|| returnType==5){
-      projection->Fit(&f1,"QNR","QNR");
+      projection->Fit(&f1,"QN","QN", vecLTM[7], vecLTM[8]);
       if (returnType==4) {
        stat= f1.GetParameter(1);
        err=f1.GetParError(1);
index 6ec2816..a6c9205 100644 (file)
@@ -40,6 +40,9 @@ class TStatToolkit : public TObject
   // HISTOGRAMS TOOLS
   //
   static  void TruncatedMean(const TH1 * his, TVectorD *param, Float_t down=0, Float_t up=1.0, Bool_t verbose=kFALSE);
+  static void MedianFilter(TH1 * his1D, Int_t nmedian);
+  static Bool_t  LTMHisto(TH1 * his, TVectorD &param , Float_t fraction=1);
+  //
   static void LTM(TH1 * his, TVectorD *param=0 , Float_t fraction=1,  Bool_t verbose=kFALSE);
   static Double_t  FitGaus(TH1* his, TVectorD *param=0, TMatrixD *matrix=0, Float_t xmin=0, Float_t xmax=0,  Bool_t verbose=kFALSE);
   static Double_t  FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param=0, TMatrixD *matrix=0, Bool_t verbose=kFALSE);
index 1b8912e..063cf6b 100644 (file)
@@ -9,7 +9,7 @@
  * without fee, provided that the above copyright notice appears in all   *
  * copies and that both the copyright notice and this permission notice   *
  * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
+ * about the suitability of this software for any prpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
 #include "AliTPCPreprocessorOffline.h"
 #include "AliTPCCorrectionFit.h"
 
+#include "AliTPCClusterParam.h"
+#include "AliTPCRecoParam.h"
+
 using std::endl;
 using std::cout;
+
 ClassImp(AliTPCPreprocessorOffline)
 
 AliTPCPreprocessorOffline::AliTPCPreprocessorOffline():
@@ -253,8 +257,16 @@ void AliTPCPreprocessorOffline::CalibTimeVdrift(const Char_t* file, Int_t ustart
     if (alignmentTime) fVdriftArray->AddLast(alignmentTime);
   }
   //
+  // 5.) Add the RecoParam and ClusterParam - for compatibility checks -different sets of parameters can invalidate calibration 
+  //
+  AliTPCClusterParam *clParam =   AliTPCcalibDB::Instance()->GetClusterParam();
+  TObjArray *recoParams = new TObjArray(4) ;
+  for (Int_t i=0;i<4;i++) recoParams->AddAt(AliTPCcalibDB::Instance()->GetRecoParam(i),i);
+  fVdriftArray->AddLast(clParam);
+  fVdriftArray->AddLast(recoParams);
+  //
   //
-  // 5. update of OCDB
+  // 6. update of OCDB
   //
   //
   UpdateOCDBDrift(ustartRun,uendRun,fOCDBstorage);
@@ -969,7 +981,7 @@ void AliTPCPreprocessorOffline::ReadGainGlobal(const Char_t* fileName){
 Bool_t AliTPCPreprocessorOffline::AnalyzeGain(Int_t startRunNumber, Int_t endRunNumber, Int_t minEntriesGaussFit,  Float_t FPtoMIPratio){
   //
   // Analyze gain - produce the calibration graphs
-  //
+ //
 
   // 1.) try to create MIP spline
   if (fGainMIP) 
@@ -1009,11 +1021,20 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGain(Int_t startRunNumber, Int_t endRun
   fGainArray->AddAt(fGraphMIP,2);
   fGainArray->AddAt(fGraphCosmic,3);
   //
-  // Add HV and PT correction parameterization
+  // 3.) Add HV and PT correction parameterization which was used
   //
   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
   if (param->GetGainSlopesHV())  fGainArray->AddLast(param->GetGainSlopesHV());
   if (param->GetGainSlopesPT())  fGainArray->AddLast(param->GetGainSlopesPT());
+  //
+  // 4.) Add the RecoParam and ClusterParam - for compatibility checks -deffrent sets of paramters can invalidate calibration 
+  //
+  AliTPCClusterParam *clParam =   AliTPCcalibDB::Instance()->GetClusterParam();
+  TObjArray *recoParams = new TObjArray(4) ;
+  for (Int_t i=0;i<4;i++) recoParams->AddAt(AliTPCcalibDB::Instance()->GetRecoParam(i),i);
+  fGainArray->AddLast(clParam);
+  fGainArray->AddLast(recoParams);
+  //
   cout << "fGraphCosmic: " << fGraphCosmic << " fGraphMIP " << fGraphMIP << endl;
   return kTRUE;
 
@@ -1224,9 +1245,9 @@ Bool_t AliTPCPreprocessorOffline::AnalyzeGainChamberByChamber(){
   // get chamber by chamber gain
   //
   if (!fGainMult) return kFALSE;
-  TGraphErrors *grShort  = fGainMult->GetGainPerChamber(0);
-  TGraphErrors *grMedium = fGainMult->GetGainPerChamber(1);
-  TGraphErrors *grLong   = fGainMult->GetGainPerChamber(2);
+  TGraphErrors *grShort  = fGainMult->GetGainPerChamberRobust(0);
+  TGraphErrors *grMedium = fGainMult->GetGainPerChamberRobust(1);
+  TGraphErrors *grLong   = fGainMult->GetGainPerChamberRobust(2);
   if (grShort==0x0 || grMedium==0x0 || grLong==0x0) {
     delete grShort;
     delete grMedium;
index 495502e..3987a4b 100644 (file)
@@ -208,13 +208,20 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa
   // First apply calibration
   //
   //  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
+  TVectorD vec(5, seed->GetParameter());
   for (Int_t irow=0;irow<159;irow++) {
     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
     if (!cluster) continue; 
     AliTPCclusterMI cl0(*cluster);
     Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
     Int_t i[1]={cluster->GetDetector()};
-
+    AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
+    Double_t ty=0,tz=0;
+     
+    if (point){
+      ty = TMath::Abs(point->GetAngleY());
+      tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));      
+    }
     transform->Transform(x,i,0,1);
     //
     // get position correction
@@ -248,7 +255,7 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa
 
 
 
-    if (fStreamLevel>2 && streamCounter<20*fStreamLevel ){
+    if (fStreamLevel>2 && gRandom->Rndm()<0.1 ){
       // dump debug info if required
       TTreeSRedirector *cstream = GetDebugStreamer();
       if (cstream){
@@ -263,6 +270,9 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa
          "cl.="<<cluster<<
          "cy="<<dy<<
          "cz="<<dz<<
+         "ty="<<ty<<
+         "tz="<<tz<<
+         "vec.="<<&vec<<  //track parameters
          "\n";
       }
     }
@@ -441,6 +451,8 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Floa
        "nclIn="<<nclIn<<
        "nclOut="<<nclOut<<
        "ncl="<<ncl<<
+       "seed.="<<seed<<
+       "track.="<<track<<
        "TrIn0.="<<trackInOld<<
        "TrOut0.="<<trackOutOld<<
        "TrIn1.="<<&trackIn<<
index 70db579..6039d59 100644 (file)
@@ -67,6 +67,7 @@ Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
 #include "AliTracker.h"
 #include "AliTPCTransform.h"
 #include "AliTPCROC.h"
+#include "TStatToolkit.h"
 
 ClassImp(AliTPCcalibGainMult)
 
@@ -1951,6 +1952,34 @@ TGraphErrors* AliTPCcalibGainMult::GetGainPerChamber(Int_t padRegion/*=1*/, Bool
   delete histGainSec;
   return gr;
 }
+TGraphErrors* AliTPCcalibGainMult::GetGainPerChamberRobust(Int_t padRegion/*=1*/, Bool_t /*plotQA=kFALSE*/)
+{
+  //
+  // Extract gain variations per chamger for 'padRegion'
+  // Use Robust mean - LTM with 60 % 0 should be similar to the truncated mean 60 %
+  if (padRegion<0||padRegion>2) return 0x0;
+  const Int_t colors[10]={1,2,4,6};
+  const Int_t markers[10]={21,25,22,20};
+  //
+  if (!fHistGainSector) return NULL;
+  if (!fHistGainSector->GetAxis(2)) return NULL;
+  fHistGainSector->GetAxis(2)->SetRangeUser(padRegion,padRegion);
+  if (padRegion==0) fHistGainSector->GetAxis(1)->SetRangeUser(0.0,35.);
+  if (padRegion>0) fHistGainSector->GetAxis(1)->SetRangeUser(36.,71.);
+  //
+  TH2D * histGainSec = fHistGainSector->Projection(0,1);
+  TGraphErrors * gr = TStatToolkit::MakeStat1D(histGainSec, 0, 0.6,2,markers[padRegion],colors[padRegion]);
+  Double_t median = TMath::Median(gr->GetN(),gr->GetY());
+  if (median>0){
+    for (Int_t i=0; i<gr->GetN();i++) {
+      gr->GetY()[i]/=median;
+      gr->SetPointError(i,gr->GetErrorX(i),gr->GetErrorY(i)/median);
+    }    
+  }
+  const char* names[3]={"SHORT","MEDIUM","LONG"};
+  gr->SetNameTitle(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]),Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]));
+  return gr;
+}
 
 // void AliTPCcalibGainMult::Terminate(){
 //   //
index 4e0eefc..cfe4817 100644 (file)
@@ -51,6 +51,7 @@ public:
   TTree *      GetdEdxTree() const {return fdEdxTree;}         // tree for the later minimization
 
   TGraphErrors* GetGainPerChamber(Int_t padRegion=1, Bool_t plotQA=kFALSE);
+  TGraphErrors* GetGainPerChamberRobust(Int_t padRegion=1, Bool_t plotQA=kFALSE);
   //
   void SetMIPvalue(Float_t mip){fMIP = mip;};
   void SetLowerTrunc(Float_t lowerTrunc){fLowerTrunc = lowerTrunc;};
index 7615146..8a897e5 100644 (file)
@@ -802,6 +802,9 @@ void AliTPCcalibSummary::ProcessGain(Int_t irun, Int_t timeStamp){
   //
   static TVectorD vGainQMaxGraphRegion(3);
   static TVectorD vGainQTotGraphRegion(3);
+
+  static TGraphErrors ggrPadEqualMax(36);
+  static TGraphErrors ggrPadEqualTot(36);
   
   vGainGraphIROC.Zero();
   vGainGraphOROCmed.Zero();
@@ -824,6 +827,12 @@ void AliTPCcalibSummary::ProcessGain(Int_t irun, Int_t timeStamp){
     TGraphErrors * graphGainIROC       = (TGraphErrors *) gainSplines->FindObject("TGRAPHERRORS_MEAN_CHAMBERGAIN_SHORT_BEAM_ALL");
     TGraphErrors * graphGainOROCMedium = (TGraphErrors *) gainSplines->FindObject("TGRAPHERRORS_MEAN_CHAMBERGAIN_MEDIUM_BEAM_ALL");
     TGraphErrors * graphGainOROCLong   = (TGraphErrors *) gainSplines->FindObject("TGRAPHERRORS_MEAN_CHAMBERGAIN_LONG_BEAM_ALL");
+    //
+    TGraphErrors *grPadEqualMax = (TGraphErrors * ) gainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
+    TGraphErrors *grPadEqualTot = (TGraphErrors * ) gainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
+    if (grPadEqualMax) ggrPadEqualMax = *grPadEqualMax;
+    if (grPadEqualTot) ggrPadEqualTot = *grPadEqualTot;
+
 
     if (graphGainIROC && graphGainOROCMedium && graphGainOROCLong) {
       Double_t x=0,y=0;
@@ -853,6 +862,8 @@ void AliTPCcalibSummary::ProcessGain(Int_t irun, Int_t timeStamp){
     
   // time dependence of gain 
   (*fPcstream)<<"dcs"<<
+    "grPadEqualMax.=" << &ggrPadEqualMax <<
+    "grPadEqualTot.=" << &ggrPadEqualTot <<
     "rocGainIROC.="            << &vGainGraphIROC        <<          
     "rocGainOROCMedium.="      << &vGainGraphOROCmed     <<
     "rocGainOROCLong.="        << &vGainGraphOROClong    <<