TPC/CalibMacros/CalibEnv.C

   1 /*
   2 //Make a tree dump of TPC calibration:
   3
   4 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
   5 .L $ALICE_ROOT/TPC/CalibMacros/CalibEnv.C+
   6
   7 CalibEnv("run.list");
   8 TFile f("dcsTime.root")
   9 */
  10 #include "TMVA/TSpline1.h"
  11 #include "TMVA/TSpline2.h"
  12 #include <iostream>
  13 #include <fstream>
  14 #include <stdio.h>
  15 #include <AliCDBManager.h>
  16 #include <AliCDBEntry.h>
  17 #include <AliLog.h>
  18 #include <AliMagF.h>
  19 #include "AliTPCcalibDB.h"
  20 #include "AliTPCcalibDButil.h"
  21 #include "AliTPCAltroMapping.h"
  22 #include "AliTPCExB.h"
  23 #include "AliTPCCalROC.h"
  24 #include "AliTPCCalPad.h"
  25 #include "AliTPCSensorTempArray.h"
  26 #include "AliGRPObject.h"
  27 #include "AliTPCTransform.h"
  28 #include "TFile.h"
  29 #include "TKey.h"
  30 #include "TObjArray.h"
  31 #include "TObjString.h"
  32 #include "TString.h"
  33 #include "AliTPCCalPad.h"
  34 #include "AliTPCROC.h"
  35 #include "AliTPCParam.h"
  36 #include "AliTPCCalibPulser.h"
  37 #include "AliTPCCalibPedestal.h"
  38 #include "AliTPCCalibCE.h"
  39 #include "AliTPCExBFirst.h"
  40 #include "TTreeStream.h"
  41 #include "AliTPCTempMap.h"
  42 #include "TVectorD.h"
  43 #include "TMatrixD.h"
  44 #include "AliTPCCalibRaw.h"
  45 #include "AliSplineFit.h"
  46 #include "TGraphErrors.h"
  47
  48 //
  49 AliTPCcalibDB     *calibDB=0;
  50 AliTPCcalibDButil *dbutil =0;
  51 TTree * dcsTree=0;
  52 TString refFile="dummy.root";
  53 TTreeSRedirector *pcstream =0;
  54 void GetNearest(TGraph *graph, Double_t x, Double_t &dx, Double_t &y);
  55 void GetInterpoloationSigma(TGraph *graph, Double_t x, Double_t &dx, Double_t &sy, Double_t deltaX);
  56 //
  57 //
  58 void ProcessRun(Int_t irun, Int_t startTime, Int_t endTime);
  59 void GetProductionInfo(Int_t run, Int_t &nalien, Int_t &nRawAlien, Int_t &nlocal, Int_t &nRawLocal);
  60 void ProcessDrift(Int_t run, Int_t timeStamp);
  61
  62 void CalibEnv(const char * runList, Int_t first=1, Int_t last=-1){
  63   //
  64   // runList - listOfRuns to process
  65   // first   - first run to process
  66   // last    - last  to process
  67   //
  68   refFile=gSystem->Getenv("REF_DATA_FILE");
  69   calibDB = AliTPCcalibDB::Instance();
  70   dbutil= new AliTPCcalibDButil;
  71   //
  72   // make list of runs
  73   //
  74   ifstream in;
  75   in.open(runList);
  76   Int_t irun=0;
  77   TArrayI runArray(100000);
  78   Int_t indexes[100000];
  79   Int_t nruns=0;
  80   while(in.good()) {
  81     in >> irun;
  82     if (in.eof()) break;
  83     if (irun<first) continue;  // process only subset of list
  84     if (last>0 && irun>=last) continue;  // process only subset of list
  85     runArray[nruns]=irun;
  86     nruns++;
  87   }
  88   TMath::Sort(nruns, runArray.fArray, indexes,kFALSE);
  89   pcstream = new TTreeSRedirector("dcsTime.root");
  90   dbutil->SetRefFile(refFile.Data());
  91   Int_t startTime = 0;
  92   Int_t endTime   = 0;
  93   for (Int_t run=0; run<nruns; run++){
  94     Int_t irun=runArray[indexes[run]];
  95     printf("Processing run %d ...\n",irun);
  96     AliTPCcalibDB::Instance()->SetRun(irun);
  97     dbutil->UpdateFromCalibDB();
  98     //
  99     AliDCSSensorArray *arrHV=calibDB->GetVoltageSensors(irun);
 100     if (!arrHV) continue;
 101     for  (Int_t isenHV=0; isenHV<arrHV->NumSensors(); ++isenHV){
 102       AliDCSSensor *senHV=arrHV->GetSensorNum(isenHV);
 103       if (!senHV) {
 104         printf("Not interesting OCDB info\n");
 105         continue;
 106       }
 107       startTime=senHV->GetStartTime();
 108       endTime  =senHV->GetEndTime();
 109       if (startTime>0&&endTime>0) break;
 110     }
 111     ProcessRun(irun, startTime,endTime);
 112   }
 113   delete pcstream;
 114 }
 115
 116
 117 void ProcessRun(Int_t irun, Int_t startTime, Int_t endTime){
 118   //
 119   //
 120   //
 121   AliSplineFit *fitVdrift=0x0;
 122   Int_t startTimeGRP=0, stopTimeGRP=0;
 123   if (calibDB->GetGRP(irun)){
 124     startTimeGRP = AliTPCcalibDB::GetGRP(irun)->GetTimeStart();
 125     stopTimeGRP  = AliTPCcalibDB::GetGRP(irun)->GetTimeEnd();
 126   }
 127
 128   AliTPCSensorTempArray * tempArray = AliTPCcalibDB::Instance()->GetTemperatureSensor(irun);
 129   AliTPCTempMap * tempMap = new AliTPCTempMap(tempArray);
 130   AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(irun);
 131   //
 132   Int_t dtime = TMath::Max((endTime-startTime)/20,10*60);
 133   //
 134   //Goofie statistical data
 135   //
 136   TVectorD vecEntries, vecMean, vecMedian,vecRMS;
 137   dbutil->ProcessGoofie(vecEntries ,vecMedian, vecMean, vecRMS);
 138   //
 139   //CE data processing - see ProcessCEdata function for description of the results
 140   //
 141   TVectorD fitResultsA, fitResultsC;
 142   Int_t nmaskedCE;
 143   Double_t chi2ACE=0,chi2CCE=0;
 144   //     dbutil->ProcessCEdata("(sector<36)++gx++gy++lx++lx**2",fitResultsA,fitResultsC,nmaskedCE);
 145   dbutil->ProcessCEdata("(sector<36)++gx++gy++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",fitResultsA,fitResultsC,nmaskedCE,chi2ACE,chi2CCE);
 146
 147   TVectorD fitCEResultsA(7), fitCEResultsC(7);
 148   Int_t    noutCE;
 149   Double_t chi2CEA=0,chi2CEC=0;
 150   AliTPCCalPad *time0 = dbutil->CreatePadTime0CE(fitCEResultsA, fitCEResultsC, noutCE, chi2CEA, chi2CEC);
 151   delete time0;
 152   //
 153   //
 154   TVectorD vecTEntries, vecTMean, vecTRMS, vecTMedian, vecQEntries, vecQMean, vecQRMS, vecQMedian;
 155   Float_t driftTimeA, driftTimeC;
 156   dbutil->ProcessCEgraphs(vecTEntries, vecTMean, vecTRMS, vecTMedian,
 157                           vecQEntries, vecQMean, vecQRMS, vecQMedian,
 158                           driftTimeA, driftTimeC );
 159   //
 160   //
 161   //
 162   //drift velocity using tracks
 163   //
 164   //     fitVdrift=calibDB->GetVdriftSplineFit("ALISPLINEFIT_MEAN_VDRIFT_COSMICS_ALL",irun);
 165   fitVdrift=calibDB->CreateVdriftSplineFit("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL",irun);
 166   //noise data Processing - see ProcessNoiseData function for description of the results
 167   TVectorD vNoiseMean, vNoiseMeanSenRegions, vNoiseRMS, vNoiseRMSSenRegions;
 168   Int_t nonMaskedZero=0;
 169   dbutil->ProcessNoiseData(vNoiseMean, vNoiseMeanSenRegions, vNoiseRMS, vNoiseRMSSenRegions, nonMaskedZero);
 170   //
 171   // comparisons
 172   //
 173   TVectorF pedestalDeviations;
 174   TVectorF noiseDeviations;
 175   TVectorF pulserQdeviations;
 176   Float_t varQMean;
 177   Int_t npadsOutOneTB;
 178   Int_t npadsOffAdd;
 179   dbutil->ProcessPedestalVariations(pedestalDeviations);
 180   dbutil->ProcessNoiseVariations(noiseDeviations);
 181   dbutil->ProcessPulserVariations(pulserQdeviations,varQMean,npadsOutOneTB,npadsOffAdd);
 182   //
 183   //L3 data
 184   //
 185   Float_t bz=AliTPCcalibDB::GetBz(irun);
 186   Char_t  l3pol=AliTPCcalibDB::GetL3Polarity(irun);
 187   //
 188   //calibration Pulser data processing
 189   //
 190   Int_t nOffChannels=0;
 191   TVectorD vTimePulser;
 192   nOffChannels=dbutil->GetNPulserOutliers();
 193   dbutil->ProcessPulser(vTimePulser);
 194   //
 195   //ALTRO data
 196   //
 197   Int_t nMasked=0;
 198   dbutil->ProcessALTROConfig(nMasked);
 199   //
 200   //Calib RAW data
 201   //
 202   Int_t nFailL1=-1;
 203   if (calibDB->GetCalibRaw()) nFailL1=calibDB->GetCalibRaw()->GetNFailL1Phase();
 204   //
 205   //production information
 206   //
 207   Int_t nalien=0,nRawAlien=0,nlocal=0,nRawLocal=0;
 208   //     GetProductionInfo(irun, nalien, nRawAlien, nlocal,nRawLocal);
 209   //run type
 210   TObjString runType(AliTPCcalibDB::GetRunType(irun).Data());
 211   //
 212   //
 213   //
 214
 215   for (Int_t itime=startTime; itime<endTime; itime+=dtime){
 216     //
 217     TTimeStamp tstamp(itime);
 218     Float_t valuePressure  = calibDB->GetPressure(tstamp,irun,0);
 219     Float_t valuePressure2 = calibDB->GetPressure(tstamp,irun,1);
 220     Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative((UInt_t)itime,irun,0);
 221     Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative((UInt_t)itime,irun,1);
 222     //temperature fits
 223     TLinearFitter * fitter = 0;
 224     TVectorD vecTemp[10];
 225     for (Int_t itype=0; itype<5; itype++)
 226       for (Int_t iside=0; iside<2; iside++){
 227         fitter= tempMap->GetLinearFitter(itype,iside,tstamp);
 228         if (!fitter) continue;
 229         fitter->Eval();
 230         fitter->GetParameters(vecTemp[itype+iside*5]);
 231         delete fitter;
 232       }
 233     //
 234     //measured skirt temperatures
 235     //
 236     TVectorD vecSkirtTempA(18);
 237     TVectorD vecSkirtTempC(18);
 238     Int_t nsenTemp=tempArray->NumSensors();
 239     for (Int_t isenTemp=0;isenTemp<nsenTemp;++isenTemp){
 240       AliTPCSensorTemp *senTemp=(AliTPCSensorTemp*)tempArray->GetSensorNum(isenTemp);
 241       if (senTemp->GetType()!=3) continue;
 242       if (TMath::Sqrt(senTemp->GetX()*senTemp->GetX()+senTemp->GetY()*senTemp->GetY())<100) continue; //only skirt, outer FC vessel
 243       Double_t val=senTemp->GetValue(tstamp);
 244       if (senTemp->GetSide()==0)
 245         vecSkirtTempA[senTemp->GetSector()]=val;
 246       else
 247         vecSkirtTempC[senTemp->GetSector()]=val;
 248     }
 249     //
 250     //goofie data
 251     //
 252     TVectorD vecGoofie;
 253     if (goofieArray){
 254       vecGoofie.ResizeTo(goofieArray->NumSensors());
 255       for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 256         AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 257         if (gsensor){
 258           vecGoofie[isensor] = gsensor->GetValue(tstamp);
 259         }
 260       }
 261     } else {
 262       vecGoofie.ResizeTo(19);
 263     }
 264     //
 265     TVectorD voltagesIROC(36);
 266     TVectorD voltagesOROC(36);
 267     for(Int_t j=1; j<36; j++) voltagesIROC[j-1] = AliTPCcalibDB::Instance()->GetChamberHighVoltage(irun, j,itime);
 268     for(Int_t j=36; j<72; j++) voltagesOROC[j-36] = AliTPCcalibDB::Instance()->GetChamberHighVoltage(irun, j,itime);
 269     Double_t voltIROC = TMath::Median(36, voltagesIROC.GetMatrixArray());
 270     Double_t voltOROC = TMath::Median(36, voltagesOROC.GetMatrixArray());
 271     //
 272     Float_t  coverIA=AliTPCcalibDB::GetCoverVoltage(irun,0,itime);
 273     Float_t  coverIC=AliTPCcalibDB::GetCoverVoltage(irun,18,itime);
 274     Float_t  coverOA=AliTPCcalibDB::GetCoverVoltage(irun,36,itime);
 275     Float_t  coverOC=AliTPCcalibDB::GetCoverVoltage(irun,54,itime);
 276     Float_t  skirtA=AliTPCcalibDB::GetSkirtVoltage(irun,0,itime);
 277     Float_t  skirtC=AliTPCcalibDB::GetSkirtVoltage(irun,18,itime);
 278     Float_t  ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(irun,0,itime);
 279     Float_t  ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(irun,18,itime);
 280     //drift velocity
 281     Float_t dvCorr=-5;
 282     if (fitVdrift) dvCorr=fitVdrift->Eval(itime);
 283     //
 284     // Gain - Alexander Kalveit
 285     //
 286     Float_t  gainCosmic = 0;
 287     Float_t  gainMIP = 0;
 288     TObjArray * gainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(irun);
 289     if (gainSplines) {
 290       TGraphErrors * graphMIP = (TGraphErrors *) gainSplines->FindObject("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL");
 291       TGraphErrors * graphCosmic = (TGraphErrors *) gainSplines->FindObject("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL");
 292       if (graphMIP) gainMIP = graphMIP->Eval(itime);
 293       if (graphCosmic) gainCosmic = graphCosmic->Eval(itime);
 294     }
 295
 296
 297     //tempMap->GetLinearFitter(0,0,itime);
 298     (*pcstream)<<"dcs"<<
 299       "run="<<irun<<
 300       "time="<<itime<<
 301       "startTimeGRP="<<startTimeGRP<<
 302       "stopTimeGRP="<<stopTimeGRP<<
 303       //run type
 304       "runType.="<<&runType<<
 305       // voltage setting
 306       "VIROC.="<<&voltagesIROC<<
 307       "VOROC.="<<&voltagesOROC<<
 308       "medianVIROC="<<voltIROC<<
 309       "medianVOROC="<<voltOROC<<
 310       "coverIA=" << coverIA <<
 311       "coverIC=" << coverIC <<
 312       "coverOA=" << coverOA <<
 313       "coverOC=" << coverOC <<
 314       "skirtA=" << skirtA <<
 315       "skirtC=" << skirtC <<
 316       "ggOffA=" << ggOffA <<
 317       "ggOffC=" << ggOffC <<
 318       //
 319       "ptrel0="<<ptrelative0<<  // deltaTP/TP  - A side
 320       "ptrel1="<<ptrelative1<<  // deltaTP/TPC - C side
 321       "goofie.="<<&vecGoofie<<
 322       "goofieE.="<<&vecEntries<<
 323       "goofieMean.="<<&vecMean<<
 324       "goofieMedian.="<<&vecMedian<<
 325       "goofieRMS.="<<&vecRMS<<
 326       //
 327       "press="<<valuePressure<<
 328       "press2="<<valuePressure2<<
 329       "temp00.="<<&vecTemp[0]<<
 330       "temp10.="<<&vecTemp[1]<<
 331       "temp20.="<<&vecTemp[2]<<
 332       "temp30.="<<&vecTemp[3]<<
 333       "temp40.="<<&vecTemp[4]<<
 334       "temp01.="<<&vecTemp[5]<<
 335       "temp11.="<<&vecTemp[6]<<
 336       "temp21.="<<&vecTemp[7]<<
 337       "temp31.="<<&vecTemp[8]<<
 338       "temp41.="<<&vecTemp[9]<<
 339       "tempSkirtA.="<<&vecSkirtTempA<<
 340       "tempSkirtC.="<<&vecSkirtTempC;
 341     ProcessDrift(irun, itime);
 342     (*pcstream)<<"dcs"<<
 343       //noise data
 344       "meanNoise.="<<&vNoiseMean<<
 345       "meanNoiseSen.="<<&vNoiseMeanSenRegions<<
 346       "rmsNoise.="<<&vNoiseRMS<<
 347       "rmsNoiseSen.="<<&vNoiseRMSSenRegions<<
 348       "zeroNoise="<<nonMaskedZero<<
 349       //pulser data
 350       "timePulser.=" << &vTimePulser <<
 351       "nOffPulser="<<nOffChannels<<
 352       //altro data
 353       "nMasked="<<nMasked<<
 354       //ce data -Jens version
 355       "CEfitA.="<<&fitResultsA<<
 356       "CEfitC.="<<&fitResultsC<<
 357       "nmaskedCE="<<nmaskedCE<<
 358       "chi2ACE="<<chi2ACE<<
 359       "chi2CCE="<<chi2CCE<<
 360       //ce data new - MI version
 361       "CEfitAMI.="<<&fitCEResultsA<<
 362       "CEfitCMI.="<<&fitCEResultsC<<
 363       "chi2CEA="<<chi2CEA<<
 364       "chi2CEC="<<chi2CEC<<
 365       //
 366       //ce graph data
 367       "CEgrTEntries.="<<&vecTEntries<<
 368       "CEgrTMean.="<<&vecTMean<<
 369       "CEgrTRMS.="<<&vecTRMS<<
 370       "CEgrTMedian.="<<&vecTMedian<<
 371       "CEgrQEntries.="<<&vecQEntries<<
 372       "CEgrQMean.="<<&vecQMean<<
 373       "CEgrQRMS.="<<&vecQRMS<<
 374       "CEgrQMedian.="<<&vecQMedian<<
 375       "CEgrDriftA="<<driftTimeA<<
 376       "CEgrDriftC="<<driftTimeC<<
 377       //calib raw data
 378       "nFailL1="<<nFailL1<<
 379       // b field
 380       "Bz="<< bz <<
 381       "L3polarity="<<l3pol<<
 382       // production information
 383       "nalien="<<nalien<<
 384       "nRawAlien="<<nRawAlien<<
 385       "nlocal="<<nlocal<<
 386       "nRawLocal="<<nRawLocal<<
 387       //comparisons with ref data
 388       "pedestalDeviations.="<<&pedestalDeviations<<
 389       "noiseDeviations.="<<&noiseDeviations<<
 390       "pulserQdeviations.="<<&pulserQdeviations<<
 391       //         "pulserVarQMean="<<varQMean<<
 392       "pulserNpadsOutOneTB="<<npadsOutOneTB<<
 393       "pulserNpadsOffAdd="<<npadsOffAdd<<
 394       "driftCorrCosmAll="<<dvCorr<<
 395       // time dependence of gain
 396       "gainMIP="<<gainMIP<<
 397       "gainCosmic="<<gainCosmic<<
 398       "\n";
 399   }//end run loop
 400   //     delete fitVdrift;
 401   //     fitVdrift=0;
 402 }
 403
 404
 405
 406
 407 void GetProductionInfo(Int_t run, Int_t &nalien, Int_t &nRawAlien, Int_t &nlocal, Int_t &nRawLocal){
 408   //
 409   // find number of ESDs in central and local production for run
 410   //
 411
 412   nalien=0;
 413   nRawAlien=0;
 414   nlocal=0;
 415   nRawLocal=0;
 416   TString sNlines;
 417   //find number of ESDs in alien
 418   TString command="alien_find /alice/data/2009 ";
 419   command += Form("%09d",run);
 420   command += " | grep AliESDs.root | wc -l";
 421   sNlines = gSystem->GetFromPipe(command.Data());
 422   nalien=sNlines.Atoi();
 423   //find number of raw files on alien
 424   command="alien_find /alice/data/2009 ";
 425   command += Form("%09d",run);
 426   command += " | grep raw | grep -v tag | wc -l";
 427   sNlines = gSystem->GetFromPipe(command.Data());
 428   nRawAlien=sNlines.Atoi();
 429   //find number of ESDs local
 430   command="find /lustre/alice/alien/alice/data/2009 -name AliESDs.root | grep ";
 431   command += Form("%09d",run);
 432   command += " | wc -l";
 433   sNlines = gSystem->GetFromPipe(command.Data());
 434   nlocal=sNlines.Atoi();
 435   //find number of local raw data files
 436   command="find /lustre/alice/alien/alice/data/2009 -name \"*.root\" | grep ";
 437   command += Form("%09d",run);
 438   command += " | grep raw | grep -v tag | wc -l";
 439   sNlines = gSystem->GetFromPipe(command.Data());
 440   nRawLocal=sNlines.Atoi();
 441 }
 442
 443
 444
 445 void ProcessDrift(Int_t run, Int_t timeStamp){
 446   //
 447   // dump drift calibration data to the tree
 448   //
 449   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
 450   TGraphErrors *laserA[3]={0,0,0};
 451   TGraphErrors *laserC[3]={0,0,0};
 452   TGraphErrors *cosmicAll=0;
 453   static Double_t     vlaserA[3]={0,0,0};
 454   static Double_t     vlaserC[3]={0,0,0};
 455   static Double_t     vcosmicAll=0;
 456
 457
 458   if (array){
 459     laserA[0]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
 460     laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
 461     laserA[2]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_A");
 462     laserC[0]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
 463     laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
 464     laserC[2]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_C");
 465     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
 466   }
 467   if (laserA[0]) vlaserA[0]= laserA[0]->Eval(timeStamp);
 468   if (laserA[1]) vlaserA[1]= laserA[1]->Eval(timeStamp);
 469   if (laserA[2]) vlaserA[2]= laserA[2]->Eval(timeStamp);
 470   if (laserC[0]) vlaserC[0]= laserC[0]->Eval(timeStamp);
 471   if (laserC[1]) vlaserC[1]= laserC[1]->Eval(timeStamp);
 472   if (laserC[2]) vlaserC[2]= laserC[2]->Eval(timeStamp);
 473   if (cosmicAll) vcosmicAll= cosmicAll->Eval(timeStamp);
 474   (*pcstream)<<"dcs"<<
 475     "vlaserA0="<<vlaserA[0]<<
 476     "vlaserA1="<<vlaserA[1]<<
 477     "vlaserA2="<<vlaserA[2]<<
 478     "vlaserC0="<<vlaserC[0]<<
 479     "vlaserC1="<<vlaserC[1]<<
 480     "vlaserC2="<<vlaserC[2]<<
 481     "vcosmicAll="<<vcosmicAll;
 482
 483   //
 484   // define distance to measurement
 485   //
 486   static Double_t dlaserA=0;
 487   static Double_t dlaserC=0;
 488   static Double_t dcosmic=0;
 489   static Double_t slaserA=0;
 490   static Double_t slaserC=0;
 491   static Double_t scosmic=0;
 492   static Double_t  vclaserA[3]={0,0,0};
 493   static Double_t  vclaserC[3]={0,0,0};
 494   static Double_t  vccosmicAll=0;
 495   Double_t dummy;
 496   Double_t deltaX=1800; // +-0.5 hour
 497   for (Int_t i=0;i<3;i++){
 498     if (laserA[i]) GetNearest(laserA[i],timeStamp,dlaserA,vclaserA[i]);
 499     if (laserC[i]) GetNearest(laserC[i],timeStamp,dlaserC,vclaserC[i]);
 500   }
 501   if (laserA[1]) GetInterpoloationSigma(laserA[1],timeStamp,dummy,slaserA,deltaX);
 502   if (laserC[1]) GetInterpoloationSigma(laserC[1],timeStamp,dummy,slaserC,deltaX);
 503
 504   if (cosmicAll) GetNearest(cosmicAll,timeStamp,dcosmic,vccosmicAll);
 505   if (cosmicAll) GetInterpoloationSigma(cosmicAll,timeStamp,dummy,scosmic,deltaX);
 506   (*pcstream)<<"dcs"<<
 507     "vclaserA0="<<vclaserA[0]<<
 508     "vclaserA1="<<vclaserA[1]<<
 509     "vclaserA2="<<vclaserA[2]<<
 510     "vclaserC0="<<vclaserC[0]<<
 511     "vclaserC1="<<vclaserC[1]<<
 512     "vclaserC2="<<vclaserC[2]<<
 513     "vccosmicAll="<<vccosmicAll<<
 514     "dlaserA="<<dlaserA<<
 515     "dlaserC="<<dlaserC<<
 516     "dcosmic="<<dcosmic<<
 517     "slaserA="<<slaserA<<
 518     "slaserC="<<slaserC<<
 519     "scosmic="<<scosmic;
 520 }
 521
 522 void GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
 523   //
 524   // find the closest point to xref  in x  direction
 525   // return dx and value
 526   Int_t index=0;
 527   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
 528   if (index<0) index=0;
 529   if (index>=graph->GetN()-1) index=graph->GetN()-2;
 530   if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
 531   dx = xref-graph->GetX()[index];
 532   y  = graph->GetY()[index];
 533 }
 534
 535
 536 void GetInterpoloationSigma(TGraph *graph, Double_t x, Double_t &/*dx*/, Double_t &sy,Double_t deltaX){
 537   //
 538   //
 539   //
 540   TMVA::TSpline1 * spline1= new TMVA::TSpline1("spline1",graph);
 541   Double_t deltas[1000];
 542   Int_t counter=0;
 543   Int_t index=0;
 544   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),x);
 545   if (index<1) index=1;
 546   if (index>=graph->GetN()-1) index=graph->GetN()-2;
 547
 548   for (Int_t idelta=-10; idelta<=10; idelta++){
 549     Double_t y0=0,dummy=0;
 550     Double_t lx=x+idelta*deltaX/20.;
 551     deltas[counter]=spline1->Eval(lx)-graph->GetY()[index];
 552     counter++;
 553     deltas[counter]=spline1->Eval(lx)-graph->GetY()[index+1];
 554     counter++;
 555   }
 556   sy=TMath::RMS(counter,deltas);
 557 }