AliTPCcalibDB.cxx - Calculate the distance to the closest measurement
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 12:04:26 +0000 (12:04 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 12:04:26 +0000 (12:04 +0000)
AliTPCcalibDButil.h AliTPCcalibDButil.cxx  - graph filtering functions
                                           - calculation of drift velocity CE
CalibMacros/CalibEnv.C                     - dump all drift velocity option

Marian

TPC/AliTPCcalibDB.cxx
TPC/AliTPCcalibDButil.cxx
TPC/AliTPCcalibDButil.h
TPC/CalibMacros/CalibEnv.C

index bf34840..07a8489 100644 (file)
@@ -1571,9 +1571,9 @@ Double_t AliTPCcalibDB::GetVDriftCorrectionTime(Int_t timeStamp, Int_t run, Int_
   //
   Double_t result;
   // mode TPC crossing and laser 
-  if (mode==1) {
-    result=AliTPCcalibDButil::GetVDriftTPC(run,timeStamp);
-    
+  Double_t deltaT=0;
+  if (mode==1) {   
+    result=AliTPCcalibDButil::GetVDriftTPC(deltaT,run,timeStamp);    
   }
 
   return result;
index 26fc9be..dbfed98 100644 (file)
@@ -30,8 +30,9 @@
 #include <TGraph.h>
 #include <TFile.h>
 #include <TDirectory.h>
-
+#include <AliCDBStorage.h>
 #include <AliDCSSensorArray.h>
+#include <AliTPCSensorTempArray.h>
 #include <AliDCSSensor.h>
 #include "AliTPCcalibDB.h"
 #include "AliTPCCalPad.h"
 #include "AliTPCParam.h"
 #include "AliTPCCalibRaw.h"
 #include "TGraphErrors.h"
-
+#include "AliLog.h"
 #include "AliTPCcalibDButil.h"
 #include "AliTPCPreprocessorOnline.h"
+#include "AliTPCCalibVdrift.h"
 
 ClassImp(AliTPCcalibDButil)
 AliTPCcalibDButil::AliTPCcalibDButil() :
@@ -78,7 +80,10 @@ AliTPCcalibDButil::AliTPCcalibDButil() :
   fCETmaxLimitAbs(1.5),
   fPulTmaxLimitAbs(1.5),
   fPulQmaxLimitAbs(5),
-  fPulQminLimit(11)
+  fPulQminLimit(11),
+  fRuns(0),                         // run list with OCDB info
+  fRunsStart(0),                    // start time for given run
+  fRunsStop(0)                     // stop time for given run
 {
   //
   // Default ctor
@@ -1354,12 +1359,13 @@ Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Dou
   return (valType==0) ? median:mean;
 }
 
-Double_t  AliTPCcalibDButil::GetVDriftTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
+Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
   //
   // Get the correction of the drift velocity
   // combining information from the laser track calibration 
   // and from cosmic calibration
   //
+  // dist      - return value - distance to closest point in graph
   // run       - run number
   // timeStamp - tim stamp in seconds
   // deltaT    - integration period to calculate time0 offset 
@@ -1381,8 +1387,11 @@ Double_t  AliTPCcalibDButil::GetVDriftTPC(Int_t run, Int_t timeStamp, Double_t d
   TGraphErrors *cosmicAll=0;
   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
   if (!cosmicAll) return 0;
+  Double_t grY=0;
+  AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
+
   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
-  Double_t vcosmic=  cosmicAll->Eval(timeStamp);
+  Double_t vcosmic=  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=timeStamp>cosmicAll->GetY()[cosmicAll->GetN()-1];
   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
   return  vcosmic+t0;
@@ -1407,3 +1416,536 @@ Double_t  AliTPCcalibDButil::GetVDriftTPC(Int_t run, Int_t timeStamp, Double_t d
   
 }
 
+Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
+  //
+  // Get the correction of the drift velocity using the laser tracks calbration
+  //
+  // run       - run number
+  // timeStamp - tim stamp in seconds
+  // deltaT    - integration period to calculate time0 offset 
+  // side      - 0 - A side,  1 - C side, 2 - mean from both sides
+  // Note in case no data form both A and C side - the value from active side used
+  TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
+  TGraphErrors *grlaserA=0;
+  TGraphErrors *grlaserC=0;
+  Double_t vlaserA=0, vlaserC=0;
+  if (!array) return 0;
+  grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
+  grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
+  Double_t deltaY;
+  if (grlaserA) {
+    AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
+    if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
+    else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
+  }
+  if (grlaserC) {
+    AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
+    if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
+    else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
+  }
+  if (side==0) return vlaserA;
+  if (side==1) return vlaserC;
+  Double_t mdrift=(vlaserA+vlaserC)*0.5;
+  if (!grlaserA) return vlaserC;
+  if (!grlaserC) return vlaserA;
+  return mdrift;
+}
+
+
+Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
+  //
+  // Get the correction of the drift velocity using the CE laser data
+  // combining information from the CE,  laser track calibration
+  // and P/T calibration 
+  //
+  // run       - run number
+  // timeStamp - tim stamp in seconds
+  // deltaT    - integration period to calculate time0 offset 
+  // side      - 0 - A side,  1 - C side, 2 - mean from both sides
+  TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
+  AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
+  TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData(); 
+  AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
+  //
+  //
+  Double_t corrPTA = 0, corrPTC=0;
+  Double_t ltime0A = 0, ltime0C=0;
+  Double_t gry=0;
+  Double_t corrA=0, corrC=0;
+  Double_t timeA=0, timeC=0;
+  TGraph *graphA = (TGraph*)arrT->At(72);
+  TGraph *graphC = (TGraph*)arrT->At(73);
+  if (!graphA && !graphC) return 0.;
+  if (graphA &&graphA->GetN()>0) {
+    AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
+    timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
+    Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
+    ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
+    if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
+    corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-3.*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
+    corrA-=corrPTA;
+  }
+  if (graphC&&graphC->GetN()>0){
+    AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
+    timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
+    Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
+    ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
+    if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
+    corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-3.*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
+    corrC-=corrPTC;
+  }
+  
+  if (side ==0 ) return corrA;
+  if (side ==1 ) return corrC;
+  Double_t corrM= (corrA+corrC)*0.5;
+  if (!graphA) corrM=corrC; 
+  if (!graphC) corrM=corrA; 
+  return corrM;
+}
+
+
+
+
+Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
+  //
+  // VERY obscure method - we need something in framework
+  // Find the TPC runs with temperature OCDB entry
+  // cache the start and end of the run
+  //
+  AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
+  if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
+  if (!storage) return 0;
+  TString path=storage->GetURI(); 
+  TString runsT;
+  {    
+    TString command;
+    if (path.Contains("local")){  // find the list if local system
+      path.ReplaceAll("local://","");
+      path+="TPC/Calib/Temperature";
+      command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
+    }
+    runsT=gSystem->GetFromPipe(command);
+  }
+  TObjArray *arr= runsT.Tokenize("r");
+  if (!arr) return 0;
+  //
+  TArrayI indexes(arr->GetEntries());
+  TArrayI runs(arr->GetEntries());
+  Int_t naccept=0;
+  {for (Int_t irun=0;irun<arr->GetEntries();irun++){
+      Int_t irunN = atoi(arr->At(irun)->GetName());
+      if (irunN<startRun) continue;
+      if (irunN>stopRun) continue;
+      runs[naccept]=irunN;
+      naccept++;
+    }}
+  fRuns.Set(naccept);
+  fRunsStart.Set(fRuns.fN);
+  fRunsStop.Set(fRuns.fN);
+  TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
+  for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
+  
+  //
+  AliCDBEntry * entry = 0;
+  {for (Int_t irun=0;irun<fRuns.fN; irun++){
+      entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
+      if (!entry) continue;
+      AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
+      if (!tmpRun) continue;
+      fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
+      fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
+      //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
+    }}
+  return fRuns.fN;
+}
+
+
+Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
+  //
+  // binary search - find the run for given time stamp
+  //
+  Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
+  Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
+  Int_t cindex  = -1;
+  for (Int_t index=index0; index<=index1; index++){
+    if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
+    if (debug) {
+      printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
+    }
+  }
+  if (cindex<0) cindex =(index0+index1)/2;
+  if (cindex<0) {
+    return 0; 
+  }
+  return fRuns[cindex];
+}
+
+
+
+
+
+TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
+  //
+  // filter outlyer measurement
+  // Only points around median +- sigmaCut filtered 
+  if (!graph) return  0;
+  Int_t kMinPoints=2;
+  Int_t npoints0 = graph->GetN();
+  Int_t npoints=0;
+  Float_t  rmsY=0;
+  Double_t *outx=new Double_t[npoints0];
+  Double_t *outy=new Double_t[npoints0];
+  //
+  //
+  if (npoints0<kMinPoints) return 0;
+  for (Int_t iter=0; iter<3; iter++){
+    npoints=0;
+    for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
+      if (graph->GetY()[ipoint]==0) continue;
+      if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;  
+      outx[npoints]  = graph->GetX()[ipoint];
+      outy[npoints]  = graph->GetY()[ipoint];
+      npoints++;
+    }
+    if (npoints<=1) break;
+    medianY  =TMath::Median(npoints,outy);
+    rmsY   =TMath::RMS(npoints,outy);
+  }
+  TGraph *graphOut=0;
+  if (npoints>1) graphOut= new TGraph(npoints,outx,outy); 
+  return graphOut;
+}
+
+
+TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
+  //
+  // filter outlyer measurement
+  // Only points around median +- cut filtered 
+  if (!graph) return  0;
+  Int_t kMinPoints=2;
+  Int_t npoints0 = graph->GetN();
+  Int_t npoints=0;
+  Float_t  rmsY=0;
+  Double_t *outx=new Double_t[npoints0];
+  Double_t *outy=new Double_t[npoints0];
+  //
+  //
+  if (npoints0<kMinPoints) return 0;
+  for (Int_t iter=0; iter<3; iter++){
+    npoints=0;
+    for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
+      if (graph->GetY()[ipoint]==0) continue;
+      if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;  
+      outx[npoints]  = graph->GetX()[ipoint];
+      outy[npoints]  = graph->GetY()[ipoint];
+      npoints++;
+    }
+    if (npoints<=1) break;
+    medianY  =TMath::Median(npoints,outy);
+    rmsY   =TMath::RMS(npoints,outy);
+  }
+  TGraph *graphOut=0;
+  if (npoints>1) graphOut= new TGraph(npoints,outx,outy); 
+  return graphOut;
+}
+
+
+
+TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
+  //
+  // filter outlyer measurement
+  // Only points with normalized errors median +- sigmaCut filtered
+  //
+  Int_t kMinPoints=10;
+  Int_t npoints0 = graph->GetN();
+  Int_t npoints=0;
+  Float_t  medianErr=0, rmsErr=0;
+  Double_t *outx=new Double_t[npoints0];
+  Double_t *outy=new Double_t[npoints0];
+  Double_t *erry=new Double_t[npoints0];
+  Double_t *nerry=new Double_t[npoints0];
+  Double_t *errx=new Double_t[npoints0];
+  //
+  //
+  if (npoints0<kMinPoints) return 0;
+  for (Int_t iter=0; iter<3; iter++){
+    npoints=0;
+    for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
+      nerry[npoints]  = graph->GetErrorY(ipoint);
+      if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;  
+      erry[npoints]  = graph->GetErrorY(ipoint);
+      outx[npoints]  = graph->GetX()[ipoint];
+      outy[npoints]  = graph->GetY()[ipoint];
+      errx[npoints]  = graph->GetErrorY(ipoint);
+      npoints++;
+    }
+    if (npoints==0) break;
+    medianErr=TMath::Median(npoints,erry);
+    medianY  =TMath::Median(npoints,outy);
+    rmsErr   =TMath::RMS(npoints,erry);
+  }
+  TGraphErrors *graphOut=0;
+  if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry); 
+  delete []outx;
+  delete []outy;
+  delete []errx;
+  delete []erry;
+  return graphOut;
+}
+
+
+void AliTPCcalibDButil::Sort(TGraph *graph){
+  //
+  // sort array - neccessay for approx
+  //
+  Int_t npoints = graph->GetN();
+  Int_t *indexes=new Int_t[npoints];
+  Double_t *outx=new Double_t[npoints];
+  Double_t *outy=new Double_t[npoints];
+  TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
+  for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
+  for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
+  for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
+  for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
+}
+void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
+  //
+  // smmoth graph - mean on the interval
+  //
+  Sort(graph);
+  Int_t npoints = graph->GetN();
+  Double_t *outy=new Double_t[npoints];
+  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+    Double_t lx=graph->GetX()[ipoint];
+    Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
+    Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
+    if (index0<0) index0=0;
+    if (index1>=npoints-1) index1=npoints-1;
+    if ((index1-index0)>1){
+      outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
+    }else{
+      outy[ipoint]=graph->GetY()[ipoint];
+    }
+  }
+  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+    graph->GetY()[ipoint] = outy[ipoint];
+  }
+  delete[] outy;
+}
+
+Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
+  //
+  // Use constant interpolation outside of range 
+  //
+  if (!graph) {
+    printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
+    return 0;
+  }
+  if (graph->GetN()<1){
+    printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
+    return 0;
+  }
+  if (xref<graph->GetX()[0]) return graph->GetY()[0];
+  if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1]; 
+  return graph->Eval( xref);
+}
+
+
+void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
+  //
+  // Filter CE data
+  // 0.  remove outlyers
+  //     0.1 absolute cut
+  //     0.2 nsigma cut
+  // 1.  smooth the graphs
+  //
+  const Int_t kMinPoints=5;  // minimal number of points to define the CE
+  TObjArray *arrT=fCalibDB->GetCErocTtime();
+  Double_t medianY=0;
+  TObjArray*  cearray =   fCalibDB->GetCEData(); 
+  if (!cearray) return;
+  AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
+  AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernPressure");
+  if ( tempMapCE && cavernPressureCE){
+    // recalculate P/T correction map for time of the CE
+    AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
+    driftCalib->SetName("driftPTCE");
+    driftCalib->SetTitle("driftPTCE");
+    cearray->AddLast(driftCalib);
+  }
+
+  for (Int_t i=0; i<arrT->GetEntries();i++){
+    TGraph *graph= (TGraph*)arrT->At(i);
+    if (!graph) continue;
+    if (graph->GetN()<kMinPoints){
+      arrT->AddAt(0,i);
+      delete graph;  // delete empty graph
+      continue;
+    }
+    TGraph* graphTS0= AliTPCcalibDButil::FilterGraphMedianAbs(graph,cutAbs,medianY);
+    if (!graphTS0) continue;
+    if (graphTS0->GetN()<kMinPoints) {
+      delete graphTS0;  
+      delete graph;
+      arrT->AddAt(0,i);
+      continue;
+    }
+    TGraph* graphTS= AliTPCcalibDButil::FilterGraphMedian(graphTS0,cutSigma,medianY);
+    
+    graphTS->Sort();
+    AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);      
+    if (pcstream){
+      Int_t run = AliTPCcalibDB::Instance()->GetRun();
+      (*pcstream)<<"filterCE"<<
+       "run="<<run<<
+       "isec="<<i<<
+       "mY="<<medianY<<
+       "graph.="<<graph<<
+       "graphTS0.="<<graphTS0<<
+       "graphTS.="<<graphTS<<
+       "\n";
+    }
+    delete graphTS0;
+    if (!graphTS) continue;
+    arrT->AddAt(graphTS,i);
+    delete graph;
+  }
+}
+
+
+void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
+  //
+  // Filter Drift velocity measurement using the tracks
+  // 0.  remove outlyers - error based
+  //     cutSigma      
+  //
+  //
+  const Int_t kMinPoints=1;  // minimal number of points to define value
+  TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
+  Double_t medianY=0;
+  if (!arrT) return;
+  for (Int_t i=0; i<arrT->GetEntries();i++){
+    TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
+    if (!graph) continue;
+    if (graph->GetN()<kMinPoints){
+      delete graph;
+      arrT->AddAt(0,i);
+      continue;
+    }
+    TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
+    if (!graph2) {
+      delete graph; arrT->AddAt(0,i); continue;
+    }
+    if (graph2->GetN()<1) {
+      delete graph; arrT->AddAt(0,i); continue;
+    }
+    graph2->SetName(graph->GetName());
+    graph2->SetTitle(graph->GetTitle());
+    arrT->AddAt(graph2,i);
+    if (pcstream){
+      (*pcstream)<<"filterTracks"<<
+       "run="<<run<<
+       "isec="<<i<<
+       "mY="<<medianY<<
+       "graph.="<<graph<<
+       "graph2.="<<graph2<<
+       "\n";
+    }
+    delete graph;
+  }
+}
+
+
+void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, TTreeSRedirector *pcstream){
+  //
+  // Filter Goofie data
+  // 0.  remove outlyers      - cutSigma around median
+  // 1.  smooth the graphs    - deltaT - smoothing window
+  //  
+  for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
+    //
+    //
+    AliDCSSensor *sensor = goofieArray->GetSensor(isensor);       
+    Double_t medianY;    
+    if (sensor &&  sensor->GetGraph()){
+      TGraph * graph = sensor->GetGraph();
+      if (isensor==3 && graph->GetN()>1){ // drift velocity
+       TGraph * graphv = AliTPCcalibDButil::FilterGraphMedianAbs(graph,0.2,medianY);
+       delete graph;
+       graph=graphv;
+      }
+      TGraph * graph2 = AliTPCcalibDButil::FilterGraphMedian(graph,cutSigma,medianY);
+      if (!graph2) continue;
+      AliTPCcalibDButil::SmoothGraph(graph2,deltaT);
+      if (pcstream){
+       Int_t run = AliTPCcalibDB::Instance()->GetRun();
+       (*pcstream)<<"filterG"<<
+         "run="<<run<<
+         "isensor="<<isensor<<
+         "mY="<<medianY<<
+         "graph.="<<graph<<
+         "graph2.="<<graph2<<
+         "\n";
+      }      
+      if (graph2->GetN()<2) {delete graph2; continue;}
+      delete graph;
+      sensor->SetGraph(graph2);      
+    }
+  }
+}
+
+
+
+
+
+
+
+Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
+  //
+  //
+  // get laser time offset 
+  // median around timeStamp+-deltaT   
+  // QA - chi2 needed for later usage - to be added
+  //    - currently cut on error
+  //
+  Int_t kMinPoints=1;
+  Double_t kMinDelay=0.01;
+  Double_t kMinDelayErr=0.0001;
+
+  TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
+  if (!array) return 0;
+  TGraphErrors *tlaser=0;
+  if (array){
+    if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
+    if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
+  }
+  if (!tlaser) return 0;
+  Int_t npoints0= tlaser->GetN();
+  if (npoints0==0) return 0;
+  Double_t *xlaser = new Double_t[npoints0];
+  Double_t *ylaser = new Double_t[npoints0];
+  Int_t npoints=0;
+  for (Int_t i=0;i<npoints0;i++){
+    //printf("%d\n",i);
+    if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros  
+    if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
+    xlaser[npoints]=tlaser->GetX()[npoints];
+    ylaser[npoints]=tlaser->GetY()[npoints];
+    npoints++;
+  }
+  //
+  //
+  Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
+  Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
+  //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
+  if (index0<0) index0=0;
+  if (index1>=npoints-1) index1=npoints-1;
+  if (index1-index0<kMinPoints) return 0;
+  //
+  //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
+    Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
+  delete [] ylaser;
+  delete [] xlaser;
+  return mean;
+}
index 315fd42..8b89585 100644 (file)
@@ -12,6 +12,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 
 #include <TObject.h>
+#include <TArrayI.h>
 
 class AliDCSSensorArray;
 class AliTPCcalibDB;
@@ -19,6 +20,8 @@ class AliTPCCalPad;
 class AliTPCmapper;
 class AliTPCCalibRaw;
 class TGraph;
+class TGraphErrors;
+class TTreeSRedirector;
 
 class AliTPCcalibDButil : public TObject
 {
@@ -86,9 +89,26 @@ public:
   //
   // graph tools
   //
-  static Double_t  GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT=86400, Double_t deltaTLaser=3600, Int_t valType=0);
-  static Double_t  GetVDriftTPC(Int_t run, Int_t timeStamp, Double_t deltaT=86400, Double_t deltaTLaser=3600, Int_t valType=0);
+  static Double_t GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side);
+  static TGraph* FilterGraphMedian(TGraph * graph, Float_t sigmaCut, Double_t &medianY);
+  static TGraph* FilterGraphMedianAbs(TGraph * graph, Float_t cut, Double_t &medianY);
+  static TGraphErrors* FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY);
+  //
+  static void Sort(TGraph *graph);
+  static void SmoothGraph(TGraph *graph, Double_t delta);
   static Int_t     GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y);
+  static Double_t EvalGraphConst(TGraph *graph, Double_t xref);
+  
+  void FilterCE(Double_t deltaT=100, Double_t cutAbs=10, Double_t cutSigma=4., TTreeSRedirector *pcstream=0);
+  void FilterTracks(Int_t run, Double_t cutSigma=20., TTreeSRedirector *pcstream=0);
+
+  void FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT=2, Double_t cutSigma=4., TTreeSRedirector *pcstream=0);
+  static Double_t  GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT=86400, Double_t deltaTLaser=3600, Int_t valType=0);
+  static Double_t  GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT=86400, Double_t deltaTLaser=3600, Int_t valType=0);
+  static Double_t  GetVDriftTPCLaserTracks(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT=43200, Int_t side=2);
+  static Double_t  GetVDriftTPCCE(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT=43200, Int_t side=2);
+  Int_t MakeRunList(Int_t startRun, Int_t stopRun); // find the list of usable runs
+  Int_t FindRunTPC(Int_t    itime, Bool_t debug=kFALSE);
 private:
   AliTPCcalibDB *fCalibDB;            //pointer to calibDB object
   AliTPCCalPad  *fPadNoise;           //noise information
@@ -128,7 +148,17 @@ private:
   Float_t fPulTmaxLimitAbs;              //maximum variation of Pulser Signals (time) before pads will be treated as outliers
   Float_t fPulQmaxLimitAbs;              //maximum variation of Pulser Signals (charge) before pads will be treated as outliers
   Float_t fPulQminLimit;                 //minimum charge value for Pulser Signals before pads will be treated as outliers
+
+  //
+  // helpers to get the run number for given time stamps
+  //
+  // filters  
   
+public:
+  TArrayI fRuns;                         // run list with OCDB info
+  TArrayI fRunsStart;                    // start time for given run
+  TArrayI fRunsStop;                     // stop time for given run
+private:
   AliTPCcalibDButil (const AliTPCcalibDButil& );
   AliTPCcalibDButil& operator= (const AliTPCcalibDButil& );
 
index 41ff7ba..4daf732 100644 (file)
@@ -51,14 +51,14 @@ AliTPCcalibDButil *dbutil =0;
 TTree * dcsTree=0;
 TString refFile="dummy.root"; 
 TTreeSRedirector *pcstream =0;
-void GetNearest(TGraph *graph, Double_t x, Double_t &dx, Double_t &y);
-void GetInterpoloationSigma(TGraph *graph, Double_t x, Double_t &dx, Double_t &sy, Double_t deltaX);
 //
 //
 void ProcessRun(Int_t irun, Int_t startTime, Int_t endTime);
 void GetProductionInfo(Int_t run, Int_t &nalien, Int_t &nRawAlien, Int_t &nlocal, Int_t &nRawLocal);
 void ProcessDrift(Int_t run, Int_t timeStamp);
-
+void ProcessDriftCE(Int_t run, Int_t timeStamp);
+void ProcessDriftAll(Int_t run, Int_t timeStamp);
+void ProcessKryptonTime(Int_t run, Int_t timeStamp);
 void CalibEnv(const char * runList, Int_t first=1, Int_t last=-1){
   //
   // runList - listOfRuns to process
@@ -108,6 +108,11 @@ void CalibEnv(const char * runList, Int_t first=1, Int_t last=-1){
       endTime  =senHV->GetEndTime();
       if (startTime>0&&endTime>0) break;
     }
+    dbutil->FilterCE(120., 3., 4.,pcstream);
+    dbutil->FilterTracks(irun, 10.,pcstream);
+    AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(irun);                 
+    //if (goofieArray) dbutil->FilterGoofie(goofieArray,2,4,pcstream);
+    // don't filter goofie for the moment
     ProcessRun(irun, startTime,endTime);
   }
   delete pcstream;  
@@ -339,6 +344,9 @@ void ProcessRun(Int_t irun, Int_t startTime, Int_t endTime){
       "tempSkirtA.="<<&vecSkirtTempA<<
       "tempSkirtC.="<<&vecSkirtTempC;
     ProcessDrift(irun, itime);
+    ProcessDriftCE(irun,itime);
+    ProcessDriftAll(irun,itime);
+    ProcessKryptonTime(irun,itime);
     (*pcstream)<<"dcs"<<       
       //noise data
       "meanNoise.="<<&vNoiseMean<<
@@ -450,10 +458,14 @@ void ProcessDrift(Int_t run, Int_t timeStamp){
   TGraphErrors *laserA[3]={0,0,0};
   TGraphErrors *laserC[3]={0,0,0};
   TGraphErrors *cosmicAll=0;
+  TGraphErrors *laserAE[3]={0,0,0};
+  TGraphErrors *laserCE[3]={0,0,0};
+  TGraphErrors *cosmicAllE=0;
   static Double_t     vlaserA[3]={0,0,0};
   static Double_t     vlaserC[3]={0,0,0};
   static Double_t     vcosmicAll=0;
-  
+  static Double_t     vdrift1=0;
+  vdrift1=AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(timeStamp,run,0,1);
 
   if (array){
     laserA[0]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
@@ -464,13 +476,13 @@ void ProcessDrift(Int_t run, Int_t timeStamp){
     laserC[2]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_C");
     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
   }
-  if (laserA[0]) vlaserA[0]= laserA[0]->Eval(timeStamp);
-  if (laserA[1]) vlaserA[1]= laserA[1]->Eval(timeStamp);
-  if (laserA[2]) vlaserA[2]= laserA[2]->Eval(timeStamp);
-  if (laserC[0]) vlaserC[0]= laserC[0]->Eval(timeStamp);
-  if (laserC[1]) vlaserC[1]= laserC[1]->Eval(timeStamp);
-  if (laserC[2]) vlaserC[2]= laserC[2]->Eval(timeStamp);
-  if (cosmicAll) vcosmicAll= cosmicAll->Eval(timeStamp); 
+  if (laserA[0]) vlaserA[0]= AliTPCcalibDButil::EvalGraphConst(laserA[0],timeStamp);
+  if (laserA[1]) vlaserA[1]= AliTPCcalibDButil::EvalGraphConst(laserA[1],timeStamp);
+  if (laserA[2]) vlaserA[2]= AliTPCcalibDButil::EvalGraphConst(laserA[2],timeStamp);
+  if (laserC[0]) vlaserC[0]= AliTPCcalibDButil::EvalGraphConst(laserC[0],timeStamp);
+  if (laserC[1]) vlaserC[1]= AliTPCcalibDButil::EvalGraphConst(laserC[1],timeStamp);
+  if (laserC[2]) vlaserC[2]= AliTPCcalibDButil::EvalGraphConst(laserC[2],timeStamp);
+  if (cosmicAll) vcosmicAll= AliTPCcalibDButil::EvalGraphConst(cosmicAll,timeStamp); 
   (*pcstream)<<"dcs"<<
     "vlaserA0="<<vlaserA[0]<<
     "vlaserA1="<<vlaserA[1]<<
@@ -478,7 +490,8 @@ void ProcessDrift(Int_t run, Int_t timeStamp){
     "vlaserC0="<<vlaserC[0]<<
     "vlaserC1="<<vlaserC[1]<<
     "vlaserC2="<<vlaserC[2]<<
-    "vcosmicAll="<<vcosmicAll;
+    "vcosmicAll="<<vcosmicAll<<
+    "vdrift1="<<vdrift1;
 
   //
   // define distance to measurement
@@ -492,17 +505,11 @@ void ProcessDrift(Int_t run, Int_t timeStamp){
   static Double_t  vclaserA[3]={0,0,0};
   static Double_t  vclaserC[3]={0,0,0};
   static Double_t  vccosmicAll=0;
-  Double_t dummy;
-  Double_t deltaX=1800; // +-0.5 hour
   for (Int_t i=0;i<3;i++){
-    if (laserA[i]) GetNearest(laserA[i],timeStamp,dlaserA,vclaserA[i]);
-    if (laserC[i]) GetNearest(laserC[i],timeStamp,dlaserC,vclaserC[i]);
+    if (laserA[i]) AliTPCcalibDButil::GetNearest(laserA[i],timeStamp,dlaserA,vclaserA[i]);
+    if (laserC[i]) AliTPCcalibDButil::GetNearest(laserC[i],timeStamp,dlaserC,vclaserC[i]);
   }  
-  if (laserA[1]) GetInterpoloationSigma(laserA[1],timeStamp,dummy,slaserA,deltaX);
-  if (laserC[1]) GetInterpoloationSigma(laserC[1],timeStamp,dummy,slaserC,deltaX);
-
-  if (cosmicAll) GetNearest(cosmicAll,timeStamp,dcosmic,vccosmicAll);
-  if (cosmicAll) GetInterpoloationSigma(cosmicAll,timeStamp,dummy,scosmic,deltaX);
+  if (cosmicAll) AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dcosmic,vccosmicAll);
   (*pcstream)<<"dcs"<<
     "vclaserA0="<<vclaserA[0]<<
     "vclaserA1="<<vclaserA[1]<<
@@ -519,39 +526,139 @@ void ProcessDrift(Int_t run, Int_t timeStamp){
     "scosmic="<<scosmic;
 }
 
-void GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
-  //
-  // find the closest point to xref  in x  direction
-  // return dx and value 
-  Int_t index=0;
-  index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
-  if (index<0) index=0;
-  if (index>=graph->GetN()-1) index=graph->GetN()-2;
-  if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
-  dx = xref-graph->GetX()[index];
-  y  = graph->GetY()[index];
-}
+void ProcessDriftCE(Int_t run,Int_t timeStamp){
+  //
+  // dump drift calibration data CE
+  //
+  TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
+  AliTPCParam *param=AliTPCcalibDB::Instance()->GetParameters();
+  static TVectorD tdriftCE(74);
+  static TVectorD tndriftCE(74);
+  static TVectorD vdriftCE(74);
+  static TVectorD tcdriftCE(74);
+  static TVectorD tddriftCE(74);
+  static Double_t ltime0A;
+  static Double_t ltime0C;
+  //
+  //
+  //
+  ltime0A  = dbutil->GetLaserTime0(run,timeStamp,36000,0);
+  ltime0C  = dbutil->GetLaserTime0(run,timeStamp,36000,1);
+  //
+  for (Int_t i=0; i<arrT->GetEntries();i++){
+    tdriftCE[i]=0;
+    vdriftCE[i]=0;
+    TGraph *graph = (TGraph*)arrT->At(i);
+    if (!graph) continue;
+    tdriftCE[i]=AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
+    Double_t deltaT,gry;
+    AliTPCcalibDButil::GetNearest(graph,timeStamp,deltaT,gry);
+    tndriftCE[i]=graph->GetN();
+    tcdriftCE[i]=gry;         
+    tddriftCE[i]=deltaT;              
+    if (i%36<18){
+      vdriftCE[i] =(param->GetZLength(i)/(tdriftCE[i]*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()))/param->GetDriftV();
+    }else{
+      vdriftCE[i] =(param->GetZLength(i)/(tdriftCE[i]*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()))/param->GetDriftV();
+    }
+  }
+
+  (*pcstream)<<"dcs"<<  
+    "tdriftCE.="<<&tdriftCE<<
+    "vdriftCE.="<<&vdriftCE<<
+    "tndriftCE.="<<&tndriftCE<<
+    "tcdriftCE.="<<&tcdriftCE<<
+    "tddriftCE.="<<&tddriftCE<<
+    "ltime0A="<<ltime0A<<
+    "ltime0C="<<ltime0C;
+   }
 
 
-void GetInterpoloationSigma(TGraph *graph, Double_t x, Double_t &/*dx*/, Double_t &sy,Double_t deltaX){
+void ProcessDriftAll(Int_t run,Int_t timeStamp){
+  //
+  // dump drift calibration data  all calibrations form DB util
+  // test of utils
+  static Double_t vdriftCEA=0, vdriftCEC=0, vdriftCEM=0;
+  static Double_t vdriftLTA=0, vdriftLTC=0, vdriftLTM=0;
+  static Double_t dce=0, dla=0,dlc=0,dlm=0;
+  static Double_t ltime0A;
+  static Double_t ltime0C;
+  static Double_t     vdrift1=0;
+  vdrift1=AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(timeStamp,run,0,1);
+  
   //
+  vdriftCEA= dbutil->GetVDriftTPCCE(dla,run,timeStamp,36000,0);
+  vdriftCEC= dbutil->GetVDriftTPCCE(dlc,run,timeStamp,36000,1);
+  vdriftCEM= dbutil->GetVDriftTPCCE(dlm,run,timeStamp,36000,2);
   //
+  vdriftLTA= dbutil->GetVDriftTPCLaserTracks(dla,run,timeStamp,36000,0);
+  vdriftLTC= dbutil->GetVDriftTPCLaserTracks(dlc,run,timeStamp,36000,1);
+  vdriftLTM= dbutil->GetVDriftTPCLaserTracks(dlm,run,timeStamp,36000,2);
   //
-  TMVA::TSpline1 * spline1= new TMVA::TSpline1("spline1",graph);
-  Double_t deltas[1000];
-  Int_t counter=0;
-  Int_t index=0;
-  index = TMath::BinarySearch(graph->GetN(), graph->GetX(),x);
-  if (index<1) index=1;
-  if (index>=graph->GetN()-1) index=graph->GetN()-2;
-  
-  for (Int_t idelta=-10; idelta<=10; idelta++){
-    Double_t y0=0,dummy=0;
-    Double_t lx=x+idelta*deltaX/20.;
-    deltas[counter]=spline1->Eval(lx)-graph->GetY()[index];
-    counter++;
-    deltas[counter]=spline1->Eval(lx)-graph->GetY()[index+1];
-    counter++;
+  ltime0A  = dbutil->GetLaserTime0(run,timeStamp,36000,0);
+  ltime0C  = dbutil->GetLaserTime0(run,timeStamp,36000,1);
+
+  (*pcstream)<<"dcs"<<  
+    //
+    "vdriftCEA="<<vdriftCEA<<
+    "vdriftCEC="<<vdriftCEC<<
+    "vdriftCEM="<<vdriftCEM<<
+    "dce="<<dce<<
+    "vdriftLTA="<<vdriftLTA<<
+    "vdriftLTC="<<vdriftLTC<<
+    "vdriftLTM="<<vdriftLTM<<
+    "dla="<<dla<<
+    "dlc="<<dlc<<
+    "dlm="<<dlm<<
+    "vdrift1="<<vdrift1;
+
+}
+
+
+
+void ProcessKryptonTime(Int_t run, Int_t timeStamp){
+  //
+  // Dumping  krypton calibration results
+  //
+  static TObjArray * krArray=0;
+  if (!krArray) {
+    AliCDBEntry* entry = AliCDBManager::Instance()->Get("TPC/Calib/TimeGainKrypton", run);
+    if (entry){
+      krArray = (TObjArray*)entry->GetObject();
+    }
   }
-  sy=TMath::RMS(counter,deltas);
+  static TVectorD krMean(74);
+  static TVectorD krErr(74);
+  static TVectorD krDist(74);
+  TGraphErrors *gr=0;
+  Double_t deltaT=0,gry=0;
+  if (krArray){
+    for (Int_t isec=0; isec<72; isec++){
+      krMean[isec]=0;
+      krDist[isec]=0;
+      krErr[isec]=0;
+      gr=(TGraphErrors*)krArray->At(isec);
+      if (gr) {
+       krMean[isec]=gr->Eval(timeStamp);
+       AliTPCcalibDButil::GetNearest(gr, timeStamp,deltaT,gry);
+       krDist[isec]=deltaT;
+      }     
+      if (72+isec<krArray->GetEntries()) {
+       gr=(TGraphErrors*)krArray->At(72+isec);
+       if (gr) krErr[isec]=gr->Eval(timeStamp);
+      }
+    }
+    krMean[72]= TMath::Median(36,krMean.GetMatrixArray());
+    krMean[73]= TMath::Median(36,&(krMean.GetMatrixArray()[36]));
+    krErr[72]= TMath::Median(36,krErr.GetMatrixArray());
+    krErr[73]= TMath::Median(36,&(krErr.GetMatrixArray()[36]));
+  }
+  (*pcstream)<<"dcs"<<
+    "krMean.="<<&krMean<<
+    "krErr.="<<&krErr<<
+    "krDist.="<<&krDist;
 }
+
+
+
+