41ff7bad3d47852b99c9ca38fb71879d163cccbf
[u/mrichter/AliRoot.git] / 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 }