3 Macro to perform fits of the Laser Central electrode data
4 Several fit methods implemented
6 1. RebuildData() - transform arbitrary layout of the Input data to the internal format
7 StoreData(); - The data tree expected in file inname (see variable bellow)
8 StoreTree(); - Modify inname and xxside and tcor in order to transform data
10 2. MakeFit(); - Make a fit of the data - already in internal format
14 3. MakeRes(); - Make the final calibration + conbination of different components
16 4. LoadViewer(); - Browse the fit parameters
22 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC -I$ALICE_ROOT/STAT");
23 gSystem->Load("libSTAT.so");
24 .L $ALICE_ROOT/TPC/CalibMacros/AnalyzeLaser.C+
27 Calibration viewer variables:
29 Result - resulting correction
34 out - outlyers not used for fit
35 tcor - offset specified by user before fitting
36 timeF1 - sector time local fit - plane
37 timeF2 - sector time local fit - parabola
38 qF1 - sector q local fit - plane
39 qF2 - sector q local fit - parabola
43 ffit1 - adding common shifts - alpha dependendent
44 ffit2 - adding opposite shifts - alpha dependent
46 fGXY - global fit parameter - XY
47 fInOut - global fit parameter - inner-outer sector matching
48 fLX - global LX dependence
53 // Control variable - check results
56 ffit2~-(timeIn~+tcor~):lx~ - fit value minus input time
59 (timeF2~-ffit2~+fTL~+fInOut~+tcor~):Result~+tcor~
61 timeF2~-Result~:ffit2~-fTL~-fInOut~
68 #include "TStatToolkit.h"
69 #include "AliTPCCalibViewer.h"
70 #include "AliTPCCalibViewerGUI.h"
71 #include "AliTPCPreprocessorOnline.h"
73 //Define interesting variables - file names
75 char * inname = "treeCE08_05-07.root"; // input file with tree
77 // variable name definition in input tree - change it according the input
79 TString qaside("CE_A00_Q_05");
80 TString taside("CE_A00_T_05");
81 TString raside("CE_A00_rms_05");
82 TString qcside("CE_C00_Q_05");
83 TString tcside("CE_C00_T_05");
84 TString rcside("CE_C00_rms_05");
86 // correction variable - usually Pulser time
88 TString tcor("(sector%36>30)*2");
91 char * fname = "treefitCE.root"; // output file with tree
92 char * oname = "fitCE.root"; // output file with CalPads fit
98 AliTPCCalPad *calPadIn = 0; // original time pad
99 AliTPCCalPad *calPadF1 = 0; // original time pad - fit plane
100 AliTPCCalPad *calPadF2 = 0; // original time pad - fit parabola
101 AliTPCCalPad *calPadQIn = 0; // original Q pad
102 AliTPCCalPad *calPadQF1 = 0; // original Q pad
103 AliTPCCalPad *calPadQF2 = 0; // original Q pad
105 AliTPCCalPad *calPadCor = 0; // base correction CalPad
106 AliTPCCalPad *calPadOut = 0; // outlyer CalPad
110 const Float_t tThr=0.5; // max diff in the sector
111 const Float_t qThr0=0.5; // max Q diff in the sector
112 const Float_t qThr1=2; // max Q diff in the sector
117 AliTPCCalPad *calPad0 = 0; // global fit 0 - base
118 AliTPCCalPad *calPad1 = 0; // global fit 1 - common behavior rotation -A-C
119 AliTPCCalPad *calPad2 = 0; // gloabl fit 2 - CE missalign rotation A-C
121 AliTPCCalPad *calPadInOut = 0; // misaalign in-out
122 AliTPCCalPad *calPadLX = 0; // local x missalign
123 AliTPCCalPad *calPadTL = 0; // tan
124 AliTPCCalPad *calPadQ = 0; // time (q) correction
125 AliTPCCalPad *calPadGXY = 0; // global XY missalign (drift velocity grad)
126 AliTPCCalPad *calPadOff = 0; // normalization offset fit
127 AliTPCCalPad *calPadRes = 0; // final calibration
133 AliTPCCalibViewerGUI * viewer=0; //viewerGUI
134 AliTPCCalibViewer *makePad=0; //viewer
135 TTree * tree=0; // working tree
138 void RebuildData(); // transform the input data to the fit format
139 void MakeFit(); // make fits
141 // internal functions
143 void MakeAliases0(); // Make Aliases 0 - for user tree
144 void MakeAliases1(); // Make Aliases 1 - for default tree
145 void LoadData(); // Load data
146 void StoreData(); // store current data
147 void StoreTree(); // store fit data in the output tree
167 TVectorD vec0,vec1,vec2;
173 fstring+="side++"; // offset on 2 different sides //1
174 //fstring+="(1/qp)++"; // Q -threshold effect correction //2
175 //fstring+="(qp)++"; // Q -threshold effect correction //3
176 fstring+="(inn)++"; // inner outer misalign - common //4
177 fstring+="(side*inn)++"; // - opposite //5
179 fstring+="(gyr)++"; // drift velocity gradient - common //6
180 fstring+="(side*gyr)++"; // - opposite //7
181 fstring+="(gxr)++"; // global x tilt - common //8
182 fstring+="(side*gxr)++"; // - opposite //9
184 fstring+="tl^2++"; // local phi correction //10
186 fstring+="(lxr)++"; // zr angle - common //11
187 fstring+="(side*lxr)++"; // - opposite //12
188 fstring+="(inn*lxr)++"; // inner outer angle - common //13
189 fstring+="(side*inn*lxr)++";// - opposite //14
190 fstring+="(lxr^2)++"; // zr second - common //15
191 fstring+="(side*lxr^2)++"; // - opposite //16
193 TString *fit0 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vec0,mat);
194 tree->SetAlias("f0",fit0->Data());
196 // Common "deformation" tendencies
198 fstring+="(sin(atan2(gy.fElements,gx.fElements)))++";
199 fstring+="(cos(atan2(gy.fElements,gx.fElements)))++";
201 fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))++";
202 fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))++";
203 fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))++";
204 fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))++";
206 fstring+="(sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
207 fstring+="(cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
208 fstring+="(sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
209 fstring+="(cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
212 TString *fit1 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&cutCE",chi2,npoints,vec1,mat);
213 tree->SetAlias("f1",fit1->Data());
215 // Central electrode "deformation"
217 fstring+="(side*sin(atan2(gy.fElements,gx.fElements)))++";
218 fstring+="(side*cos(atan2(gy.fElements,gx.fElements)))++";
220 fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))++";
221 fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))++";
222 fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))++";
223 fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))++";
225 fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*2))*lxr++";
226 fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*2))*lxr++";
227 fstring+="(side*sin(atan2(gy.fElements,gx.fElements)*3))*lxr++";
228 fstring+="(side*cos(atan2(gy.fElements,gx.fElements)*3))*lxr++";
230 TString *fit2 =stat.FitPlane(tree,"dt",fstring.Data(),"cutF&&abs(dt-f0)<0.7&&cutCE",chi2,npoints,vec2,mat);
231 tree->SetAlias("f2",fit2->Data());
235 TString tmpstr = fstring;
236 TObjArray *arr = tmpstr.Tokenize("++");
237 TString fitQ("0"); // q correction
238 TString fitLX("0"); // lx correction
239 TString fitInOut("0"); // inner-outer - match
240 TString fitGXY("0"); // global xy fit
241 TString fitOff("0"); // side offsets
242 TString fitTL("0"); // side offsets
249 for(Int_t i=0;i<arr->GetEntriesFast();i++){
250 if (!arr->At(i)) continue;
251 TString *fitstr = new TString(arr->At(i)->GetName());
253 //Bool_t isQ = fitstr->Contains("qp)");
254 Bool_t isRot = fitstr->Contains("sin(")+fitstr->Contains("cos(");
255 Bool_t isLX = fitstr->Contains("lxr");
256 Bool_t isIn = fitstr->Contains("inn");
257 Bool_t isGXY = fitstr->Contains("gxr")+fitstr->Contains("gyr");
258 if (fitstr->Contains("tl^2")){
260 fitTL+=(*fitstr)+"*";
265 fitGXY+=(*fitstr)+"*";
271 // fitQ+=(*fitstr)+"*";
275 if (isLX&&!isRot&&!isIn){
277 fitLX+=(*fitstr)+"*";
283 fitInOut+=(*fitstr)+"*";
289 tree->SetAlias("fInOut",fitInOut.Data());
290 tree->SetAlias("fLX",fitLX.Data());
291 tree->SetAlias("fGXY",fitGXY.Data());
292 tree->SetAlias("fOff",fitOff.Data());
293 //tree->SetAlias("fQ",fitQ.Data());
294 tree->SetAlias("fTL",fitTL.Data());
299 calPad0 = makePad->GetCalPad("f0","1", "ffit0");
300 calPad1 = makePad->GetCalPad("f1","1", "ffit1");
301 calPad2 = makePad->GetCalPad("f2","1", "ffit2");
302 calPadInOut = makePad->GetCalPad("fInOut","1", "fInOut");
303 calPadLX = makePad->GetCalPad("fLX","1", "fLX");
304 calPadTL = makePad->GetCalPad("fTL","1", "fTL");
305 //calPadQ = makePad->GetCalPad("fQ","1", "fQ");
306 calPadGXY = makePad->GetCalPad("fGXY","1", "fGXY");
307 calPadOff = makePad->GetCalPad("fOff","1", "fOff");
314 TObjArray * array = AliTPCCalibViewerGUI::ShowGUI(fname);
315 viewer = (AliTPCCalibViewerGUI*)array->At(0);
316 makePad = viewer->GetViewer();
317 tree = viewer->GetViewer()->GetTree();
328 // transform the input data to the fit format
330 makePad = new AliTPCCalibViewer(inname);
331 tree = makePad->GetTree();
334 calPadCor = makePad->GetCalPad("tcor","1", "tcor");
335 calPadOut = makePad->GetCalPad("1","!((cutA||cutC)&&abs(ly.fElements/lx.fElements)<0.155)", "out");
336 calPadIn = makePad->GetCalPad("dt-tcor","(cutA||cutC)&&abs(ly.fElements/lx.fElements)<0.155","timeIn");
337 calPadF1 = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
338 calPadQIn = makePad->GetCalPad("qa*(sector%36<18)+qc*(sector%36>17)","1","qIn");
340 // Update outlyer maps
342 for (Int_t isector=0;isector<72; isector++){
343 for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
344 Float_t val0= calPadIn->GetCalROC(isector)->GetValue(ich);
345 Float_t val1= calPadF1->GetCalROC(isector)->GetValue(ich);
346 if (TMath::Abs(val0-val1)>tThr) calPadOut->GetCalROC(isector)->SetValue(ich,1);
349 calPadF1 = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
350 calPadF2 = calPadIn->GlobalFit("timeF2", calPadOut,kTRUE,1,0.9);
351 calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kTRUE,0,0.9);
352 calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1);
354 // Update outlyer maps
356 for (Int_t isector=0;isector<72; isector++){
357 for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
358 Float_t val0= calPadQIn->GetCalROC(isector)->GetValue(ich);
359 Float_t val1= calPadQF2->GetCalROC(isector)->GetValue(ich);
361 calPadOut->GetCalROC(isector)->SetValue(ich,1);
364 if (TMath::Abs(val0/val1)<qThr0) calPadOut->GetCalROC(isector)->SetValue(ich,1);
365 if (TMath::Abs(val0/val1)>qThr1) calPadOut->GetCalROC(isector)->SetValue(ich,1);
368 calPadF1 = calPadIn->GlobalFit("timeF1", calPadOut,kTRUE,0,0.9);
369 calPadF2 = calPadIn->GlobalFit("timeF2", calPadOut,kTRUE,1,0.9);
370 calPadQF1 = calPadQIn->GlobalFit("qF1", calPadOut,kTRUE,0,0.9);
371 calPadQF2 = calPadQIn->GlobalFit("qF2", calPadOut,kFALSE,1);
379 calPadIn = (AliTPCCalPad*)f.Get("timeIn"); // original time pad
380 calPadF1 = (AliTPCCalPad*)f.Get("timeF1"); // original time pad - fit plane
381 calPadF2 = (AliTPCCalPad*)f.Get("timeF2"); // original time pad - fit parabola
383 calPadQIn = (AliTPCCalPad*)f.Get("qIn"); // original time pad
384 calPadQF1 = (AliTPCCalPad*)f.Get("qF1"); // original time pad - fit plane
385 calPadQF2 = (AliTPCCalPad*)f.Get("qF2"); // original time pad - fit parabola
387 calPadCor = (AliTPCCalPad*)f.Get("tcor"); // base correction CalPad
388 calPadOut = (AliTPCCalPad*)f.Get("out"); // outlyer CalPad
390 calPad0 = (AliTPCCalPad*)f.Get("ffit0"); // global fit 0 - base
391 calPad1 = (AliTPCCalPad*)f.Get("ffit1"); // global fit 1 - common behavior rotation -A-C
392 calPad2 = (AliTPCCalPad*)f.Get("ffit2"); // gloabl fit 2 - CE missalign rotation A-C
393 calPadInOut = (AliTPCCalPad*)f.Get("fInOut");// misaalign in-out
394 calPadLX = (AliTPCCalPad*)f.Get("fLX"); // local x missalign
395 calPadTL = (AliTPCCalPad*)f.Get("fTL"); // local y/x missalign
396 calPadQ = (AliTPCCalPad*)f.Get("fQ"); // time (q) correction
397 calPadGXY = (AliTPCCalPad*)f.Get("fGXY"); // global XY missalign (drift velocity grad)
398 calPadOff = (AliTPCCalPad*)f.Get("fOff"); // normalization offset fit
399 calPadRes = (AliTPCCalPad*)f.Get("Result"); //resulting calibration
406 TFile * fstore = new TFile(oname,"recreate");
407 if (calPadIn) calPadIn->Write("timeIn"); // original time pad
408 if (calPadF1) calPadF1->Write("timeF1"); // original time pad - fit plane
409 if (calPadF2) calPadF2->Write("timeF2"); // original time pad - fit parabola
411 if (calPadQIn) calPadQIn->Write("qIn"); // original time pad
412 if (calPadQF1) calPadQF1->Write("qF1"); // original time pad - fit plane
413 if (calPadQF2) calPadQF2->Write("qF2"); // original time pad - fit parabola
415 if (calPadCor) calPadCor->Write("tcor"); // base correction CalPad
416 if (calPadOut) calPadOut->Write("out"); // outlyer CalPad
418 if (calPad0) calPad0->Write("ffit0"); // global fit 0 - base
419 if (calPad1) calPad1->Write("ffit1"); // global fit 1 - common behavior rotation -A-C
420 if (calPad2) calPad2->Write("ffit2"); // gloabl fit 2 - CE missalign rotation A-C
421 if (calPadInOut)calPadInOut->Write("fInOut"); // misaalign in-out
422 if (calPadLX) calPadLX->Write("fLX"); // local x missalign
423 if (calPadTL) calPadTL->Write("fTL"); // local y/x missalign
424 if (calPadQ) calPadQ->Write("fQ"); // time (q) correction
425 if (calPadGXY) calPadGXY->Write("fGXY"); // global XY missalign (drift velocity grad)
426 if (calPadOff) calPadOff->Write("fOff"); // normalization offset fit
427 if (calPadRes) calPadRes->Write("Result"); //resulting calibration
436 AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
438 if (calPadIn) preprocesor->AddComponent(calPadIn->Clone());
439 if (calPadF1) preprocesor->AddComponent(calPadF1->Clone());
440 if (calPadF2) preprocesor->AddComponent(calPadF2->Clone());
442 if (calPadQIn) preprocesor->AddComponent(calPadQIn->Clone());
443 if (calPadQF1) preprocesor->AddComponent(calPadQF1->Clone());
444 if (calPadQF2) preprocesor->AddComponent(calPadQF2->Clone());
446 if (calPadCor) preprocesor->AddComponent(calPadCor->Clone());
447 if (calPadOut) preprocesor->AddComponent(calPadOut->Clone());
449 if (calPad0) preprocesor->AddComponent(calPad0->Clone());
450 if (calPad1) preprocesor->AddComponent(calPad1->Clone());
451 if (calPad2) preprocesor->AddComponent(calPad2->Clone());
452 if (calPadInOut)preprocesor->AddComponent(calPadInOut->Clone());
453 if (calPadLX) preprocesor->AddComponent(calPadLX->Clone());
454 if (calPadTL) preprocesor->AddComponent(calPadTL->Clone());
455 if (calPadQ) preprocesor->AddComponent(calPadQ->Clone());
456 if (calPadGXY) preprocesor->AddComponent(calPadGXY->Clone());
457 if (calPadOff) preprocesor->AddComponent(calPadOff->Clone());
458 if (calPadRes) preprocesor->AddComponent(calPadRes->Clone());
459 preprocesor->DumpToFile(fname);
466 // Define variables and selection of outliers - for user defined tree
468 tree->SetAlias("tcor",tcor.Data()); // correction variable
469 tree->SetAlias("ta",taside+".fElements");
470 tree->SetAlias("tc",tcside+".fElements");
471 tree->SetAlias("qa",qaside+".fElements");
472 tree->SetAlias("qc",qcside+".fElements");
473 tree->SetAlias("ra",raside+".fElements");
474 tree->SetAlias("rc",rcside+".fElements");
475 tree->SetAlias("side","1-(sector%36>17)*2");
476 tree->SetAlias("dt","(ta)*(sector%36<18)+(tc)*(sector%36>17)+tcor");
477 tree->SetAlias("cutA","qa>30&&qa<400&&abs(ta)<2&&ra>0.5&&ra<2");
478 tree->SetAlias("cutC","qc>30&&qc<400&&abs(tc)<2&&rc>0.5&&rc<2");
479 tree->SetAlias("cutF","(pad.fElements%4==0)&&(row.fElements%3==0)");
480 tree->SetAlias("cutCE","V.out.fElements");
484 tree->SetAlias("inn","sector<36");
485 tree->SetAlias("gxr","(gx.fElements/250.)"); //
486 tree->SetAlias("gyr","(gy.fElements/250.)"); //
487 tree->SetAlias("lxr","(lx.fElements-133.41)/250.");
488 tree->SetAlias("qp","((sector%36<18)*sqrt(qa)/10.+(sector%36>17)*sqrt(qc)/10.)"); //
489 tree->SetAlias("tl","(ly.fElements/lx.fElements)/0.17");
495 // Define variables and selection of outliers -for default usage
497 tree->SetAlias("tcor","tcor.fElements"); // correction variable
498 tree->SetAlias("side","1-(sector%36>17)*2");
499 tree->SetAlias("dt","timeIn.fElements+tcor");
501 tree->SetAlias("cutA","out.fElements==1");
502 tree->SetAlias("cutC","out.fElements==1");
503 tree->SetAlias("cutF","(pad.fElements%5==0)&&(row.fElements%4==0)");
504 tree->SetAlias("cutCE","out.fElements<0.5");
508 tree->SetAlias("inn","sector<36");
509 tree->SetAlias("gxr","(gx.fElements/250.)"); //
510 tree->SetAlias("gyr","(gy.fElements/250.)"); //
511 tree->SetAlias("lxr","(lx.fElements-133.41)/250.");
512 tree->SetAlias("qp","(sqrt(qIn.fElements)/10.+(out.fElements>0.5))"); //
513 tree->SetAlias("tl","(ly.fElements/lx.fElements)/0.17");
520 // make final calibration
522 AliTPCCalPad * calPadRes0 =new AliTPCCalPad(*calPadIn);
523 calPadRes0->Add(calPadCor); // add correction
524 calPadRes0->Add(calPad2,-1); // remove global fit
525 calPadRes = calPadRes0->GlobalFit("Result", calPadOut,kTRUE,1,0.9);
529 Float_t tlmedian = calPadTL->GetMedian();
530 for (Int_t isector=0;isector<72; isector++){
531 for (UInt_t ich=0;ich<calPadIn->GetCalROC(isector)->GetNchannels();ich++){
534 Float_t val0 = calPadRes->GetCalROC(isector)->GetValue(ich);
535 if (TMath::Abs(val0)>0.5) calPadRes->GetCalROC(isector)->SetValue(ich,0);
536 Float_t tl = calPadTL->GetCalROC(isector)->GetValue(ich);
537 Float_t inOut = calPadInOut->GetCalROC(isector)->GetValue(ich);
538 calPadRes->GetCalROC(isector)->SetValue(ich,calPadRes->GetCalROC(isector)->GetValue(ich)+tl+inOut);
543 calPadRes->Add(calPadCor,-1); // remove back correction (e.g Pulser time 0)