]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibTimeVdrift.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding OCDB configuration for different calibration use cases
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibTimeVdrift.C
CommitLineData
191306e9 1/* small macro to generate and update OCDB entries for a given run:
15e6432c 2
191306e9 3// How to use it...
4//
5// 1. go to the directory with the "CalibObjectsTrain1.root" file and start your aliroot
6//
7// 2. copy and paste the following commands
15e6432c 8//
15e6432c 9gSystem->Load("libSTEER");
10gSystem->Load("libANALYSIS");
11gSystem->Load("libSTAT");
12gSystem->Load("libTPCcalib");
13gSystem->AddIncludePath("-I$ALICE_ROOT/STEER");
14gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
15e6432c 15.L $ALICE_ROOT/TPC/CalibMacros/CalibTimeVdrift.C+
16
191306e9 17//
18// 3. setup your OCDB output path e.g:
19//
20ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
21//
22// 4. simple execution
23//
24CalibTimeVdriftGlobal()
25//
26// 5. try to visualize new entry
27//
28ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
29AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
30AliCDBManager::Instance()->SetSpecificStorage("/TPC/Calib/TimeDrift",ocdbStorage.Data());
31Int_t run=90000; //
32AliCDBEntry* entry = AliCDBManager::Instance()->Get("/TPC/Calib/TimeDrift",run)
33TObjArray * arr = (TObjArray*)entry->GetObject();
34TObjArray *picArray = new TObjArray;
35MakeDefaultPlots(arr,picArray);
36
37
15e6432c 38
15e6432c 39*/
40
191306e9 41#include "TMap.h"
15e6432c 42#include "TGraphErrors.h"
43#include "AliExternalTrackParam.h"
44#include "TFile.h"
45#include "TGraph.h"
191306e9 46#include "TMultiGraph.h"
47#include "TCanvas.h"
48#include "THnSparse.h"
49#include "TLegend.h"
50#include "TPad.h"
51
15e6432c 52
53#include "AliTPCcalibTime.h"
54#include "AliSplineFit.h"
55#include "AliCDBMetaData.h"
56#include "AliCDBId.h"
57#include "AliCDBManager.h"
58#include "AliCDBStorage.h"
191306e9 59#include "AliTPCcalibBase.h"
60#include "AliTPCcalibDB.h"
61#include "AliTPCcalibDButil.h"
62#include "AliRelAlignerKalman.h"
15e6432c 63
191306e9 64TString ocdbStorage="dummy";
65//
66//
67//
68void AddHistoGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift, Int_t minEntries);
69void AddAlignmentGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift);
70void AddLaserGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift);
71//
15e6432c 72
191306e9 73TGraphErrors* FilterGraphMedianAbs(TGraphErrors * graph, Float_t cut,Double_t &medianY);
74TGraphErrors* FilterGraphDrift(TGraphErrors * graph, Float_t errSigmaCut, Float_t medianCutAbs);
15e6432c 75
15e6432c 76
15e6432c 77
15e6432c 78
191306e9 79void CalibTimeVdriftGlobal(Char_t* file="CalibObjectsTrain1.root"){
15e6432c 80
191306e9 81 const Int_t kMinEntries=500; // minimal number of entries
82 //
83 // 1. Initialization and run range setting
84 TFile fcalib(file);
85 AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
86 Int_t startRun=0, endRun=0;
87 Int_t startTime=0, endTime=0;
88 //Int_t startPT=-1, endPT=-1;
89 //
90 // find the fist and last run
91 //
92 TObjArray *hisArray =timeDrift->GetHistoDrift();
93 {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
94 THnSparse* addHist=(THnSparse*)hisArray->UncheckedAt(i);
95 if (addHist->GetEntries()<kMinEntries) continue;
96 if (!addHist) continue;
97 TH1D* histo =addHist->Projection(3);
98 TH1D* histoTime=addHist->Projection(0);
99 if (startRun==0){
100 startRun=histo->FindFirstBinAbove(0);
101 endRun =histo->FindLastBinAbove(0);
102 }else{
103 startRun=TMath::Min(histo->FindFirstBinAbove(0),startRun);
104 endRun =TMath::Max(histo->FindLastBinAbove(0),endRun);
105 }
106 if (startTime==0){
107 startTime=histoTime->FindFirstBinAbove(0);
108 endTime =histoTime->FindLastBinAbove(0);
109 }else{
110 startTime=TMath::Min(histoTime->FindFirstBinAbove(0),startTime);
111 endTime =TMath::Max(histoTime->FindLastBinAbove(0),endTime);
112 }
113 delete histo;
114 delete histoTime;
115 }}
116 {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
117 THnSparse* addHist=(THnSparse*)hisArray->At(i);
118 if (!addHist) continue;
119 addHist->GetAxis(0)->SetRange(startTime-1,endTime+1);
120 addHist->GetAxis(3)->SetRange(startRun-1,endRun+1);
121 }
122 }
123 printf("Run range :\t%d-%d\n", startRun, endRun);
124 printf("Time range :\t%d-%d\n", startTime, endTime);
125 //
126 // 2. extraction of the information
127 //
128 TObjArray * vdriftArray = new TObjArray();
129 AddAlignmentGraphs(vdriftArray,timeDrift);
130 AddHistoGraphs(vdriftArray,timeDrift,kMinEntries);
131 AddLaserGraphs(vdriftArray,timeDrift);
132 //
133 //
134 // 3. update of OCDB
15e6432c 135 //
15e6432c 136 //
191306e9 137 AliCDBMetaData *metaData= new AliCDBMetaData();
138 metaData->SetObjectClassName("TObjArray");
139 metaData->SetResponsible("Dag Toppe Larsen");
140 metaData->SetBeamPeriod(1);
141 metaData->SetAliRootVersion("05-25-01"); //root version
142 metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
143 AliCDBId* id1=NULL;
144 id1=new AliCDBId("TPC/Calib/TimeDrift", startRun, AliCDBRunRange::Infinity());
145 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
146 gStorage->Put(vdriftArray, (*id1), metaData);
147
148}
149
150
151void CalibTimeVdrift(Int_t runNumber, AliTPCcalibTime* vdrift, Int_t end=false){
152 TObjArray* vdriftArray = new TObjArray();
153 TObjArray * array=vdrift->GetHistoDrift();
154 if(!array) return;
155 TIterator* iterator = array->MakeIterator();
156 iterator->Reset();
157 THnSparse* hist=NULL;
158 while((hist=(THnSparseF*)iterator->Next())){
159 if(!hist) continue;
160 hist->Print();
161 if(end) hist->GetAxis(3)->SetRangeUser(runNumber-0.5, end+0.5);
162 else hist->GetAxis(3)->SetRangeUser(runNumber-0.5, runNumber+0.5);
163 TString name=hist->GetName();
164 Int_t dim[4]={0,1,2,3};
165 THnSparse* newHist=hist->Projection(4,dim);
166 newHist->SetName(name);
167 vdriftArray->Add(newHist);
168 TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
169 printf("name=%s graph=%i\n", name.Data(), graph==0);
170 if(!graph || !graph->GetN()) continue;
171 printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
172 Int_t pos=name.Index("_");
173 name=name(pos,name.Capacity()-pos);
174 TString graphName=graph->ClassName();
175 graphName+=name;
176 graphName.ToUpper();
177 graph->SetName(graphName);
178 printf("name=%s\n", graphName.Data());
179 vdriftArray->Add(graph);
180// AliSplineFit* fit=new AliSplineFit();
181// fit->SetGraph(graph);
182// fit->SetMinPoints(graph->GetN()+1);
183// fit->InitKnots(graph,2,0,0.001);
184// fit->SplineFit(0);
185// TString fiName=fit->ClassName();
186// fiName+=type;
187// fiName+=trigger;
188// fiName.ToUpper();
189// fit->SetName(fiName.Data());
190// printf("name=%s\n", fiName.Data());
191// vdriftArray->Add(fit);
192 }
193 THnSparse* laserHist=NULL;
194 TGraphErrors* laserGraph=NULL;
195 TString laserName="";
196
197 laserHist=vdrift->GetHistVdriftLaserA(1);
198 laserName=laserHist->ClassName();
199 laserName+="_MEAN_DRIFT_LASER_ALL_A";
200 laserName.ToUpper();
201 laserHist->SetName(laserName);
202 vdriftArray->Add(laserHist);
203 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
204 if(laserGraph && laserGraph->GetN()){
205 laserName=laserGraph->GetName();
206 laserName+="_MEAN_DRIFT_LASER_ALL_A";
207 laserName.ToUpper();
208 laserGraph->SetName(laserName);
209 vdriftArray->Add(laserGraph);
210 }
211
212 laserHist=vdrift->GetHistVdriftLaserC(1);
213 laserName=laserHist->ClassName();
214 laserName+="_MEAN_DRIFT_LASER_ALL_C";
215 laserName.ToUpper();
216 laserHist->SetName(laserName);
217 vdriftArray->Add(laserHist);
218 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
219 if(laserGraph && laserGraph->GetN()){
220 laserName=laserGraph->GetName();
221 laserName+="_MEAN_DRIFT_LASER_ALL_C";
222 laserName.ToUpper();
223 laserGraph->SetName(laserName);
224 vdriftArray->Add(laserGraph);
225 }
15e6432c 226
227 AliCDBMetaData *metaData= new AliCDBMetaData();
228 metaData->SetObjectClassName("TObjArray");
229 metaData->SetResponsible("Dag Toppe Larsen");
230 metaData->SetBeamPeriod(1);
191306e9 231 metaData->SetAliRootVersion("05-25-01"); //root version
15e6432c 232 metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
191306e9 233 AliCDBId* id1=NULL;
234 if(end) id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, end);
235 else id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, runNumber);
236 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
237 gStorage->Put(vdriftArray, (*id1), metaData);
238 printf("done runNumber=%i, end=%i\n", runNumber, end);
239}
240
241void CalibrateAll(Char_t* file="CalibObjectsTrain1.root"){
242 TFile fcalib(file);
243 AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
244 THnSparse* addHist=timeDrift->GetHistoDrift("all");
245 TH1D* histo=addHist->Projection(3);
246 Int_t start=histo->FindFirstBinAbove(0);
247 Int_t end =histo->FindLastBinAbove(0);
248 CalibTimeVdrift(start-1, timeDrift, end-1);
249 for(Int_t i=start;i<end;i++){
250 if(histo->GetBinContent(i)){
251 printf("start=%i i=%i end=%i i=%i center=%f", start, i, end, i, histo->GetBinCenter(i));
252 CalibTimeVdrift(i, timeDrift);
253 }
254 }
255}
256
257void MakePlot(Char_t* file="CalibObjects.root", Char_t* trigger="", /*Int_t useError=true,*/ Int_t separateFits=false, Int_t pointLimit=0, Int_t showHisto=true, Int_t otherColour=false, Int_t spline=true, Double_t yMin=0.0, Double_t yMax=0.0){
258 TFile fcalib(file);
259 AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
15e6432c 260
191306e9 261//Get graphs
262 TObjArray* grArray=new TObjArray();
263 TObjArray* hiArray=new TObjArray();
264 Int_t counter=0;
265 Int_t pos=counter+1;
266 Int_t colour=2;
267 TObjArray * addArray=timeDrift->GetHistoDrift();
268 if(!addArray) return;
269 TIterator* iterator = addArray->MakeIterator();
270 iterator->Reset();
271 THnSparse* addHist=NULL;
15e6432c 272
191306e9 273 THnSparse* laserHist=NULL;
274 TString laserName="";
275
276 laserHist=timeDrift->GetHistVdriftLaserA(1);
277 laserName=laserHist->ClassName();
278 laserName+="_MEAN_VDRIFT_LASER_ALL_A";
279 laserName.ToUpper();
280 laserHist->SetName(laserName);
281 addArray->Add(laserHist);
282
283 laserHist=timeDrift->GetHistVdriftLaserC(1);
284 laserName=laserHist->ClassName();
285 laserName+="_MEAN_VDRIFT_LASER_ALL_C";
286 laserName.ToUpper();
287 laserHist->SetName(laserName);
288 addArray->Add(laserHist);
289
290 while((addHist=(THnSparseF*)iterator->Next())){
291 if (!addHist) continue;
292 TString histName=addHist->GetName();
293 if(!histName.Contains(trigger) && !histName.Contains("COSMICS_ALL")) continue;
294 addHist->Print();
295 TGraph* graph=AliTPCcalibBase::FitSlices(addHist,2,0,400,100,0.05,0.95, kTRUE);
296 TH2D* histo=addHist->Projection(2,0);
297 if(!graph || !histo || !graph->GetN()) continue;
298 if(spline){
299 AliSplineFit fit;
300 fit.SetGraph(graph);
301 fit.SetMinPoints(graph->GetN()+1);
302 fit.InitKnots(graph,2,0,0.001);
303 fit.SplineFit(0);
304 graph=fit.MakeGraph(graph->GetX()[0],graph->GetX()[graph->GetN()-1],50000,0);
305 }
306 graph->SetNameTitle(addHist->GetName(),"Relative drift velocity change");
307 graph->GetXaxis()->SetTitle("time");
308 graph->GetXaxis()->SetTimeFormat();
309 graph->GetXaxis()->SetTimeDisplay(1);
310 graph->GetYaxis()->SetTitleOffset(1.2);
311 graph->GetXaxis()->SetLabelSize(0.035);
312 graph->GetYaxis()->SetLabelSize(0.04);
313 graph->GetYaxis()->SetTitle("relative drift velocity change");
314 if(yMin) graph->SetMinimum(yMin);
315 if(yMax) graph->SetMaximum(yMax);
316 printf("bins=%i name=%s couter=%i pos=%i colour=%i Ymax=%f Ymin=%f Xmax=%f Xmin=%f\n", graph->GetXaxis()->GetNbins(), addHist->GetName(), counter, pos, colour, graph->GetYaxis()->GetXmax(), graph->GetYaxis()->GetXmin(), graph->GetXaxis()->GetXmax(), graph->GetXaxis()->GetXmin());
317 if(graph->GetXaxis()->GetNbins()<=pointLimit) continue;
318// if(!useError&&!spline) for(Int_t j=0; j<graph->GetXaxis()->GetNbins(); j++) graph->SetPointError(j,0,0);
319// if( *((Int_t*)(graph->GetName()))==*((Int_t*)("all")) ){
320 TString name=graph->GetName();
321 if(name.Contains("COSMICS_ALL")){
322 graph->SetLineColor(0+1);
323 graph->SetMarkerColor(0+1);
324 grArray->AddAtAndExpand(graph,0);
325 histo->SetLineColor(1);
326 histo->SetMarkerColor(1);
327 hiArray->AddAtAndExpand(histo,0);
328 pos=counter;
329 }
330 else{
331 graph->SetLineColor(colour);
332 graph->SetMarkerColor(colour);
333 grArray->AddAtAndExpand(graph,pos);
334 histo->SetLineColor(colour);
335 histo->SetMarkerColor(colour);
336 hiArray->AddAtAndExpand(histo,pos);
337 colour++;
338 if(colour==10) colour++;
339 }
340 pos++;
341 counter++;
342 graph=0;
343 addHist=0;
344 }
345
346 if(!separateFits){
347 TCanvas* plot = new TCanvas("plot", "Relative drift velocity change", 1500, 900);
348 //Draw graphs
349 TGraph* grTemp=0;
350 TH2D* hiTemp=0;
351 for(Int_t i=0; i<counter; i++){
352 if(!(*grArray)[i] || !(*hiArray)[i]) continue;
353 grTemp=(TGraph*)(*grArray)[i];
354 hiTemp=(TH2D*)(*hiArray)[i];
355 if(i==0) grTemp->Draw("alp");
356 else grTemp->Draw("lp");
357 if(showHisto) hiTemp->Draw("same");
358 }
359 //Add legends
360 TLegend* leg=NULL;
361 if(pointLimit){
362 leg = new TLegend(0.4,0.75,1,0.95,NULL,"brNDC");
363 leg->SetTextSize(0.03);
364 }
365 else{
366 leg = new TLegend(0.4,0.65,1,1,NULL,"brNDC");
367 leg->SetTextSize(0.02);
368 }
369 leg->SetLineColor(0);
370 leg->SetLineWidth(0);
371 leg->SetFillColor(0);
372 TLegendEntry* entry=0;
373 for(Int_t i=0; i<counter; i++){
374 if(!(*grArray)[i]) continue;
375 TString name=(*grArray)[i]->GetName();
376 Int_t pos=name.Index("_")+1+12;
377 name=name(pos,name.Capacity()-pos);
378 entry = leg->AddEntry((*grArray)[i],name,"l");
379 leg->Draw();
380 }
381 plot->SaveAs("relative_drift_velocity_change.pdf");
382 plot->SaveAs("relative_drift_velocity_change.png");
383 }
384 else{
385 //Draw graphs
386 TGraph* grTemp=0;
387 TH2D* hiTemp=0;
388 for(Int_t i=0; i<counter; i++){
389 if(!(*grArray)[i] || !(*hiArray)[i]) continue;
390 grTemp=(TGraph*)(*grArray)[i];
391 hiTemp=(TH2D*)(*hiArray)[i];
392 TString name=grTemp->GetName();
393 Int_t pos=name.Index("_")+1+12;
394 name=name(pos,name.Capacity()-pos);
395 grTemp->SetTitle(name);
396 if(otherColour){
397 grTemp->SetLineColor(2);
398 grTemp->SetMarkerColor(2);
399 hiTemp->SetLineColor(1);
400 hiTemp->SetMarkerColor(1);
401 }
402 TCanvas* plot = new TCanvas(name, name, 1500, 900);
403 grTemp->Draw("alp");
404 if(showHisto) hiTemp->Draw("same");
405 TString name2=name;
406 name+=".pdf";
407 name2+=".png";
408 plot->SaveAs(name);
409 plot->SaveAs(name2);
410 }
411 }
15e6432c 412}
413
191306e9 414void MakeQaPlot(Char_t* file="CalibObjects.root"){
415 TFile fcalib(file);
416 AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
417 TCanvas* plotQa = new TCanvas("plotQa", "Relative drift velocity change QA", 1600, 600);
418 plotQa->Divide(5,2);
419 for(Int_t i=0; i<10; i++){
420 plotQa->cd(i+1);
421 timeDrift->GetCosmiMatchingHisto(i)->Draw();
422 }
423 plotQa->cd(0);
424 plotQa->SaveAs("relative_drift_velocity_change_qa.pdf");
425}
15e6432c 426
191306e9 427void MakeHistPlot(Char_t* file="CalibObjects.root", Char_t* trigger="COSMICS_ALL"){
428 TFile fcalib(file);
429 AliTPCcalibTime* timeDrift=(AliTPCcalibTime*)fcalib.Get("calibTime");
430 THnSparse* addHist=timeDrift->GetHistoDrift(trigger);
431 TH2D* histo=addHist->Projection(2,0);
432
433 histo->SetTitle("Relative drift velocity change");
434 histo->GetXaxis()->SetTitle("time");
435 histo->GetXaxis()->SetTimeFormat();
436 histo->GetXaxis()->SetTimeDisplay(1);
437 histo->GetYaxis()->SetTitleOffset(1.3);
438 histo->GetYaxis()->SetTitle("relative drift velocity change");
439 histo->SetLineColor(1);
440 histo->SetMarkerColor(1);
441 TGraphErrors* graph=AliTPCcalibBase::FitSlices(addHist,2,0,400,100,0.05,0.95, kTRUE);
442 if(!graph) return;
443 graph->SetNameTitle(addHist->GetName(),"Relative drift velocity change");
444 graph->GetXaxis()->SetTitle("time");
445 graph->GetXaxis()->SetTimeFormat();
446 graph->GetXaxis()->SetTimeDisplay(1);
447 graph->GetYaxis()->SetTitleOffset(1.3);
448 graph->GetYaxis()->SetTitle("relative drift velocity change");
449 graph->SetMaximum(0.003);
450 graph->SetMinimum(-0.003);
451 graph->SetLineColor(2);
452 graph->SetMarkerColor(2);
453 TCanvas* histo_fit = new TCanvas("histo_fit", "Relative drift velocity change histogram/fit", 800, 600);
454 histo->Draw();
455 graph->Draw("lp");
456 histo_fit->SaveAs("relative_drift_velocity_change_histo_fit.pdf");
457 histo_fit->SaveAs("relative_drift_velocity_change_histo_fit.png");
458}
459
460
461
462void PrintArray(TObjArray *array){
463 //
464 //
465 //
466 Int_t entries = array->GetEntries();
467 for (Int_t i=0; i<entries; i++){
468 printf("%s\n", array->At(i)->GetName());
469 }
470}
15e6432c 471
191306e9 472
473TMultiGraph * MakeMultiGraph(TObjArray *array, const char *mask, Int_t color, Int_t style){
474 //
475 // Make a multi graph with graphs selected
476 //
477 if (!array) return 0;
478 TMultiGraph *mgraph=new TMultiGraph("MultiGraph", "MultiGraph");
479 TObjArray * maskArray=0;
480 if (mask){
481 TString grmask(mask);
482 maskArray = grmask.Tokenize("*");
483 }
484 Double_t ymin=0,ymax=-1;
485 for (Int_t i=0; i<array->GetEntries();i++){
486 TGraphErrors *graph= (TGraphErrors*)array->At(i);
487 if (!graph) continue;
488 if (maskArray){
489 Bool_t isOK=kTRUE;
490 TString str(graph->GetName());
491 for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
492 if (!str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
493 if (!isOK) continue;
494 }
495 if (ymax<ymin){
496 ymin=graph->GetMinimum();
497 ymax=graph->GetMaximum();
498 }
499 ymin = TMath::Min(ymin,graph->GetMinimum());
500 ymax = TMath::Max(ymax,graph->GetMaximum());
501 }
502 ymin=ymin-(ymax-ymin)*0.2;
503 ymax=ymax+(ymax-ymin)*0.5;
504 Bool_t first=kTRUE;
505 for (Int_t i=0; i<array->GetEntries();i++){
506 TGraphErrors *graph= (TGraphErrors*)array->At(i);
507 if (!graph) continue;
508 if (maskArray){
509 Bool_t isOK=kTRUE;
510 TString str(graph->GetName());
511 for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
512 if (str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
513 if (!isOK) continue;
514 }
515 graph->SetMarkerColor(color);
516 graph->SetMarkerStyle(style);
517 graph->SetLineColor(color);
518 graph->SetMaximum(ymax);
519 graph->SetMinimum(ymin);
520 if (first) {first=kFALSE;}
521 mgraph->Add(graph);
522 }
523 return mgraph;
524}
525
526
527TGraphErrors* FilterGraphDrift(TGraphErrors * graph, Float_t errSigmaCut, Float_t medianCutAbs){
528 // 2 filters:
529 // 1. filter graph - error cut errSigmaCut
530 // 2. filter graph - medianCutAbs around median
531 //
532 // errSigmaCut - cut on error
533 // medianCutAbs - cut on value around median
534 Double_t dummy=0; //
535 //
536 // 1. filter graph - error cut errSigmaCut
537 //
538 TGraphErrors *graphF;
539 graphF = AliTPCcalibDButil::FilterGraphMedianErr(graph,errSigmaCut,dummy);
540 delete graph;
541 if (!graphF) return 0;
542 graph = AliTPCcalibDButil::FilterGraphMedianErr(graphF,errSigmaCut,dummy);
543 delete graphF;
544 if (!graph) return 0;
545 //
546 // filter graph - kMedianCutAbs around median
547 //
548 graphF=FilterGraphMedianAbs(graph, medianCutAbs,dummy);
549 delete graph;
550 if (!graphF) return 0;
551 graph=FilterGraphMedianAbs(graphF, medianCutAbs,dummy);
552 delete graphF;
553 if (!graph) return 0;
554 return graph;
555}
556
557
558
559TGraphErrors* FilterGraphMedianAbs(TGraphErrors * graph, Float_t cut,Double_t &medianY){
560 //
561 // filter outlyer measurement
562 // Only points around median +- cut filtered
563 //
564 if (!graph) return 0;
565 Int_t kMinPoints=2;
566 Int_t npoints0 = graph->GetN();
567 Int_t npoints=0;
568 Float_t rmsY=0;
569 Double_t *outx=new Double_t[npoints0];
570 Double_t *outy=new Double_t[npoints0];
571 Double_t *errx=new Double_t[npoints0];
572 Double_t *erry=new Double_t[npoints0];
573 //
574 //
575 if (npoints0<kMinPoints) return 0;
576 for (Int_t iter=0; iter<3; iter++){
577 npoints=0;
578 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
579 if (graph->GetY()[ipoint]==0) continue;
580 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
581 outx[npoints] = graph->GetX()[ipoint];
582 outy[npoints] = graph->GetY()[ipoint];
583 errx[npoints] = graph->GetErrorX(ipoint);
584 erry[npoints] = graph->GetErrorY(ipoint);
585 npoints++;
586 }
587 if (npoints<=1) break;
588 medianY =TMath::Median(npoints,outy);
589 rmsY =TMath::RMS(npoints,outy);
590 }
591 TGraphErrors *graphOut=0;
592 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
593 return graphOut;
594}
595
596
597void AddHistoGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift, Int_t minEntries){
598 //
599 // Add graphs corresponding to the alignment
600 //
601 const Double_t kErrSigmaCut=5; // error sigma cut - for filtering
602 const Double_t kMedianCutAbs=0.03; // error sigma cut - for filtering
603 //
604 TObjArray * array=timeDrift->GetHistoDrift();
605 if (array){
606 THnSparse* hist=NULL;
607 // 2.a) cosmics with different triggers
608 for (Int_t i=0; i<array->GetEntriesFast();i++){
609 hist=(THnSparseF*)array->UncheckedAt(i);
610 if(!hist) continue;
611 if (hist->GetEntries()<minEntries) continue;
612 //hist->Print();
613 TString name=hist->GetName();
614 Int_t dim[4]={0,1,2,3};
615 THnSparse* newHist=hist->Projection(4,dim);
616 newHist->SetName(name);
617 TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
618 printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
619 Int_t pos=name.Index("_");
620 name=name(pos,name.Capacity()-pos);
621 TString graphName=graph->ClassName();
622 graphName+=name;
623 graphName.ToUpper();
624 //
625 graph = FilterGraphDrift(graph, kErrSigmaCut, kMedianCutAbs);
626 if (!graph) {
627 printf("Graph =%s filtered out\n", name.Data());
628 continue;
629 }
630 //
631 graph->SetMarkerStyle(i%8+20);
632 graph->SetMarkerColor(i%7);
633 graph->GetXaxis()->SetTitle("Time");
634 graph->GetYaxis()->SetTitle("v_{dcor}");
635 graph->SetName(graphName);
636 graph->SetTitle(graphName);
637 printf("Graph %d\t=\t%s\n", i, graphName.Data());
638 vdriftArray->Add(graph);
639 }
640 }
641}
642
643
644
645
646void AddAlignmentGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift){
647 //
648 // Add graphs corresponding to alignment to the object array
649 //
650 TObjArray *arrayITS=0;
651 TObjArray *arrayTOF=0;
652 TObjArray *arrayTRD=0;
653 TMatrixD *mstatITS=0;
654 TMatrixD *mstatTOF=0;
655 TMatrixD *mstatTRD=0;
656 //
657 arrayITS=timeDrift->GetAlignITSTPC();
658 arrayTRD=timeDrift->GetAlignTRDTPC();
659 arrayTOF=timeDrift->GetAlignTOFTPC();
660 mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.9,50,0.025);
661 mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.9,1000,0.025);
662 mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.9,50,0.025);
663 //
664 TObjArray * arrayITSP= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 0, 5.);
665 TObjArray * arrayITSM= AliTPCcalibDButil::SmoothRelKalman(arrayITS,*mstatITS, 1, 5.);
666 TObjArray * arrayTOFP= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,*mstatTOF, 0, 5.);
667 TObjArray * arrayTOFM= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,*mstatTOF, 1, 5.);
668 TObjArray * arrayTRDP= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,*mstatTRD, 0, 5.);
669 TObjArray * arrayTRDM= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,*mstatTRD, 1, 5.);
670 //
671 //
672 Int_t entries=TMath::Max(arrayITS->GetEntriesFast(),arrayTOF->GetEntriesFast());
673 TObjArray *arrays[9]={arrayITS, arrayITSP, arrayITSM,
674 arrayTRD, arrayTRDP, arrayTRDM,
675 arrayTOF, arrayTOFP, arrayTOFM};
676 TString grnames[9]={"ALIGN_ITS", "ALIGN_ITSP", "ALIGN_ITSM",
677 "ALIGN_TRD", "ALIGN_TRDP", "ALIGN_TRDM",
678 "ALIGN_TOF", "ALIGN_TOFP", "ALIGN_TOFM"};
679 TString grpar[9]={"DELTAPSI", "DELTATHETA", "DELTAPHI",
680 "DELTAX", "DELTAY", "DELTAZ",
681 "DRIFTVD", "T0", "VDGY"};
682
683
684 TVectorD vX(entries);
685 TVectorD vY(entries);
686 TVectorD vEx(entries);
687 TVectorD vEy(entries);
688 TObjArray *arr=0;
689 for (Int_t iarray=0; iarray<9; iarray++){
690 arr = arrays[iarray];
691 if (arr==0) continue;
692 for (Int_t ipar=0; ipar<9; ipar++){
693 Int_t counter=0;
694 for (Int_t itime=0; itime<arr->GetEntriesFast(); itime++){
695 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) arr->UncheckedAt(itime);
696 if (!kalman) continue;
697 vX[counter]=kalman->GetTimeStamp();
698 vY[counter]=(*(kalman->GetState()))[ipar];
699 if (ipar==6) vY[counter]=1./(*(kalman->GetState()))[ipar]-1;
700 vEx[counter]=0;
701 vEy[counter]=TMath::Sqrt((*(kalman->GetStateCov()))(ipar,ipar));
702 counter++;
703 }
704
705 TGraphErrors * graph=new TGraphErrors(counter, vX.GetMatrixArray(),
706 vY.GetMatrixArray(),
707 vEx.GetMatrixArray(),
708 vEy.GetMatrixArray());
709 TString grName=grnames[iarray];
710 grName+="_TPC_";
711 grName+=grpar[ipar];
712 graph->SetName(grName.Data());
713 vdriftArray->AddLast(graph);
714 }
715 }
716}
717
718
719
720void AddLaserGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift){
721 //
722 // Add laser graphs
723 //
724 THnSparse* laserHist=NULL;
725 TGraphErrors* laserGraph=NULL;
726 TString laserName="";
727 laserHist=timeDrift->GetHistVdriftLaserA(1);
728 laserName=laserHist->ClassName();
729 laserName+="_MEAN_DRIFT_LASER_ALL_A";
730 laserName.ToUpper();
731 laserHist->SetName(laserName);
732 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,300,10,0.05,0.95, kTRUE);
733 laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
734 //
735 if(laserGraph && laserGraph->GetN()){
736 laserName=laserGraph->GetName();
737 laserName+="_MEAN_DRIFT_LASER_ALL_A";
738 laserName.ToUpper();
739 laserGraph->SetName(laserName);
740 vdriftArray->Add(laserGraph);
741 }
742
743 laserHist=NULL;
744 laserGraph=NULL;
745 laserName="";
746
747 laserHist=timeDrift->GetHistVdriftLaserC(1);
748 laserName=laserHist->ClassName();
749 laserName+="_MEAN_DRIFT_LASER_ALL_C";
750 laserName.ToUpper();
751 laserHist->SetName(laserName);
752 //vdriftArray->Add(laserHist);// really add the whole histogram ??
753 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,300,10,0.05,0.95, kTRUE);
754 laserGraph = FilterGraphDrift(laserGraph, 5,0.04);
755 if(laserGraph && laserGraph->GetN()){
756 laserName=laserGraph->GetName();
757 laserName+="_MEAN_DRIFT_LASER_ALL_C";
758 laserName.ToUpper();
759 laserGraph->SetName(laserName);
760 vdriftArray->Add(laserGraph);
761 }
762
763 // 2.c) additional informaion from laser for A and C side : DELAY
764 laserHist=NULL;
765 laserGraph=NULL;
766 laserName="";
767
768 laserHist=timeDrift->GetHistVdriftLaserA(0);
769 laserName=laserHist->ClassName();
770 laserName+="_MEAN_DELAY_LASER_ALL_A";
771 laserName.ToUpper();
772 laserHist->SetName(laserName);
773 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
15e6432c 774
191306e9 775 if(laserGraph && laserGraph->GetN()){
776 laserName=laserGraph->GetName();
777 laserName+="_MEAN_DELAY_LASER_ALL_A";
778 laserName.ToUpper();
779 laserGraph->SetName(laserName);
780 vdriftArray->Add(laserGraph);
781 }
782
783 laserHist=NULL;
784 laserGraph=NULL;
785 laserName="";
786
787 laserHist=timeDrift->GetHistVdriftLaserC(0);
788 laserName=laserHist->ClassName();
789 laserName+="_MEAN_DELAY_LASER_ALL_C";
790 laserName.ToUpper();
791 laserHist->SetName(laserName);
792 //vdriftArray->Add(laserHist);// really add the whole histogram ??
793 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
794 if(laserGraph && laserGraph->GetN()){
795 laserName=laserGraph->GetName();
796 laserName+="_MEAN_DELAY_LASER_ALL_C";
797 laserName.ToUpper();
798 laserGraph->SetName(laserName);
799 vdriftArray->Add(laserGraph);
800 }
801
802 // 2.d) additional informaion from laser for A and C side : GlobalY Gradient
803 laserHist=NULL;
804 laserGraph=NULL;
805 laserName="";
806
807 laserHist=timeDrift->GetHistVdriftLaserA(2);
808 laserName=laserHist->ClassName();
809 laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_A";
810 laserName.ToUpper();
811 laserHist->SetName(laserName);
812 //vdriftArray->Add(laserHist); // really add the whole histogram ??
813 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
814 laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
815 if(laserGraph && laserGraph->GetN()){
816 laserName=laserGraph->GetName();
817 laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_A";
818 laserName.ToUpper();
819 laserGraph->SetName(laserName);
820 vdriftArray->Add(laserGraph);
821 }
822
823 laserHist=NULL;
824 laserGraph=NULL;
825 laserName="";
826
827 laserHist=timeDrift->GetHistVdriftLaserC(2);
828 laserName=laserHist->ClassName();
829 laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_C";
830 laserName.ToUpper();
831 laserHist->SetName(laserName);
832 laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,100,10,0.05,0.95, kTRUE);
833 laserGraph=FilterGraphDrift(laserGraph, 5,0.04);
834 if(laserGraph && laserGraph->GetN()){
835 laserName=laserGraph->GetName();
836 laserName+="_MEAN_GLOBALYGRADIENT_LASER_ALL_C";
837 laserName.ToUpper();
838 laserGraph->SetName(laserName);
839 vdriftArray->Add(laserGraph);
840 }
841}
842
843void SetDefaultGraphDrift(TGraph *graph, Int_t color, Int_t style){
844 //
845 //
846 //
847 graph->GetXaxis()->SetTimeDisplay(kTRUE);
848 graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
849 graph->SetMaximum( 0.025);
850 graph->SetMinimum(-0.025);
851 graph->GetXaxis()->SetTitle("Time");
852 graph->GetYaxis()->SetTitle("v_{dcorr}");
853 //
854 graph->GetYaxis()->SetLabelSize(0.03);
855 graph->GetXaxis()->SetLabelSize(0.03);
856 //
857 graph->GetXaxis()->SetNdivisions(10,5,0);
858 graph->GetYaxis()->SetNdivisions(10,5,0);
859 //
860 graph->GetXaxis()->SetLabelOffset(0.02);
861 graph->GetYaxis()->SetLabelOffset(0.005);
862 //
863 graph->GetXaxis()->SetTitleOffset(1.3);
864 graph->GetYaxis()->SetTitleOffset(1.2);
865 //
866 graph->SetMarkerColor(color);
867 graph->SetLineColor(color);
868 graph->SetMarkerStyle(style);
869}
870
871void SetPadStyle(TPad *pad, Float_t mx0, Float_t mx1, Float_t my0, Float_t my1){
872 //
873 //
874 pad->SetTicks(1,1);
875 pad->SetMargin(mx0,mx1,my0,my1);
876}
877
878
879void MakeDefaultPlots(TObjArray * arr, TObjArray *picArray){
880 //
881 //
882 //
883 // margins
884 Float_t mx0=0.12, mx1=0.1, my0=0.15, my1=0.1;
15e6432c 885 //
191306e9 886 TGraphErrors* laserA =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
887 TGraphErrors* laserC =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
888 TGraphErrors* cosmic =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
889 TGraphErrors* cross =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_VDRIFT_CROSS_ALL");
890 TGraphErrors* itstpcP =(TGraphErrors*)arr->FindObject("ALIGN_ITSP_TPC_DRIFTVD");
891 TGraphErrors* itstpcM =(TGraphErrors*)arr->FindObject("ALIGN_ITSM_TPC_DRIFTVD");
15e6432c 892 //
191306e9 893 if (laserA) SetDefaultGraphDrift(laserA,2,25);
894 if (laserC) SetDefaultGraphDrift(laserC,4,26);
895 if (cosmic) SetDefaultGraphDrift(cosmic,3,27);
896 if (cross) SetDefaultGraphDrift(cross,4,28);
897 if (itstpcP) SetDefaultGraphDrift(itstpcP,2,29);
898 if (itstpcM) SetDefaultGraphDrift(itstpcM,4,30);
899 //
900 //
901 TPad *pad=0;
902 //
903 // Laser-Laser
904 //
905 if (laserA&&laserC){
906 pad = new TCanvas("TPCLaserVDrift","TPCLaserVDrift");
907 laserA->Draw("alp");
908 SetPadStyle(pad,mx0,mx1,my0,my1);
909 laserA->Draw("apl");
910 laserC->Draw("p");
911 TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
912 legend->AddEntry(laserA,"Laser A side");
913 legend->AddEntry(laserC,"Laser C side");
914 legend->Draw();
915 picArray->AddLast(pad);
916 }
917
918 if (itstpcP&&itstpcM){
919 pad = new TCanvas("ITSTPC","ITSTPC");
920 itstpcP->Draw("alp");
921 SetPadStyle(pad,mx0,mx1,my0,my1);
922 itstpcP->Draw("alp");
923 gPad->Clear();
924 itstpcM->Draw("apl");
925 itstpcP->Draw("p");
926 TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
927 legend->AddEntry(itstpcP,"ITS-TPC smooth forward");
928 legend->AddEntry(itstpcM,"ITS-TPC smooth back");
929 legend->Draw();
930 picArray->AddLast(pad);
931 }
932
933 if (itstpcP&&laserA){
934 pad = new TCanvas("ITSTPC_LASER","ITSTPC_LASER");
935 SetPadStyle(pad,mx0,mx1,my0,my1);
936 laserA->Draw("alp");
937 itstpcM->Draw("p");
938 TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
939 legend->AddEntry(laserA,"TPC laser");
940 legend->AddEntry(itstpcM,"ITS-TPC smooth back");
941 legend->Draw();
942 picArray->AddLast(pad);
943 }
944
945 if (itstpcP&&cross){
946 pad = new TCanvas("ITSTPC_CROSS","ITSTPC_CROSS");
947 SetPadStyle(pad,mx0,mx1,my0,my1);
948 itstpcP->Draw("alp");
949 pad->Clear();
950 cross->Draw("ap");
951 itstpcP->Draw("p");
952 //
953 TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
954
955 legend->AddEntry(cross,"TPC cross tracks");
956 legend->AddEntry(itstpcP,"ITS-TPC smooth back");
957 legend->Draw();
958 picArray->AddLast(pad);
959 }
960 if (itstpcP&&cosmic){
961 pad = new TCanvas("ITSTPC_COSMIC","ITSTPC_COSMIC");
962 SetPadStyle(pad,mx0,mx1,my0,my1);
963 itstpcP->Draw("alp");
964 pad->Clear();
965 cosmic->Draw("ap");
966 itstpcP->Draw("p");
967 //
968 TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
969
970 legend->AddEntry(cosmic,"TPC cross tracks0 up-down");
971 legend->AddEntry(itstpcP,"ITS-TPC smooth back");
972 legend->Draw();
973 picArray->AddLast(pad);
974 }
975}
976
977
978TObjArray * SelectEntries(TObjArray* array, char * mask){
979 //
980 //
981 //
982 TObjArray * arraySelect = new TObjArray;
983 TObjArray * maskArray=0;
984 if (mask){
985 TString grmask(mask);
986 maskArray = grmask.Tokenize("*");
987 }
988 for (Int_t i=0; i<array->GetEntries();i++){
989 TObject *graph= array->At(i);
990 if (!graph) continue;
991 if (maskArray){
992 Bool_t isOK=kTRUE;
993 TString str(graph->GetName());
994 for (Int_t imask=0; imask<maskArray->GetEntries();imask++)
995 if (str.Contains(maskArray->At(imask)->GetName())==0) isOK=kFALSE;
996 if (!isOK) continue;
997 if (isOK) {
998 printf("%s\n",graph->GetName());
999 arraySelect->AddLast(graph->Clone());
1000 }
1001 }
1002 }
1003 return arraySelect;
1004}
1005
1006
1007void Substract(TGraphErrors *refgr, TObjArray *arraySelect){
15e6432c 1008 //
15e6432c 1009 //
191306e9 1010 /*
1011 TObjArray *arraySelect = SelectEntries(arr,"COSMIC");
1012 TGraphErrors *refgr= (TGraphErrors*)arr->FindObject("ALIGN_ITSP_TPC_DRIFTVD");
1013 Substract(refgr,arraySelect);
1014 */
1015 Float_t mx0=0.12, mx1=0.1, my0=0.15, my1=0.5;
1016 {
1017 for (Int_t igr=0; igr<arraySelect->GetEntries(); igr++){
1018 TGraphErrors *gr = (TGraphErrors *)arraySelect->At(igr);
1019 Int_t ngr = gr->GetN();
1020 for (Int_t i=0; i<ngr;i++){
1021 gr->GetY()[i]-=refgr->Eval(gr->GetX()[i]);
1022 gr->GetY()[i]*=250/0.264;
1023 }
1024 gr->GetXaxis()->SetRangeUser(refgr->GetXaxis()->GetXmin(),refgr->GetXaxis()->GetXmax());
1025 SetDefaultGraphDrift(gr,((igr)),20+igr);
1026 gr->SetMaximum(8.);
1027 gr->SetMinimum(4.);
1028 gr->GetYaxis()->SetTitle("#Delta_{t}(Time Bin)");
1029 }
15e6432c 1030 //
191306e9 1031 TPad * pad = new TCanvas("Delays","Delays",1000,800);
1032 SetPadStyle(pad,mx0,mx1,my0,my1);
1033 TLegend *legend = new TLegend(mx0+0.01,1-my1, 1-mx1, 1-0.01, "Time Offset");
1034 {
1035 arraySelect->At(0)->Draw("ap");
1036 for (Int_t igr=0; igr<arraySelect->GetEntries(); igr++){
1037 TGraphErrors *gr = (TGraphErrors *)arraySelect->At(igr);
1038 if (gr->GetN()<5) continue;
1039 if (gr) gr->Draw("p");
1040 legend->AddEntry(gr,gr->GetName());
1041 }
1042 }
1043 legend->Draw();
1044 }
15e6432c 1045}