Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 :
17 :
18 : /*
19 : Responsible: marian.ivanov@cern.ch
20 : Code to analyze the TPC calibration and to produce OCDB entries
21 :
22 :
23 : .x ~/rootlogon.C
24 : gSystem->Load("libANALYSIS");
25 : gSystem->Load("libTPCcalib");
26 :
27 : AliTPCPreprocessorOffline proces;
28 : TString ocdbPath="local:////"
29 : ocdbPath+=gSystem->GetFromPipe("pwd");
30 :
31 : proces.CalibTimeGain("CalibObjects.root",run0,run1,ocdbPath);
32 : proces.CalibTimeVdrift("CalibObjects.root",run0,run1,ocdbPath);
33 : // take the raw calibration data from the file CalibObjects.root
34 : // and make a OCDB entry with run validity run0-run1
35 : // results are stored at the ocdbPath - local or alien ...
36 : // default storage ""- data stored at current working directory
37 :
38 : e.g.
39 : gSystem->Load("libANALYSIS");
40 : gSystem->Load("libTPCcalib");
41 : AliTPCPreprocessorOffline proces;
42 : proces.CalibTimeGain("TPCMultObjects.root",114000,140040,0);
43 : TFile oo("OCDB/TPC/Calib/TimeGain/Run114000_121040_v0_s0.root")
44 : TObjArray * arr = AliCDBEntry->GetObject()
45 : arr->At(4)->Draw("alp")
46 :
47 : */
48 : #include "Riostream.h"
49 : #include <fstream>
50 : #include "TMap.h"
51 : #include "TGraphErrors.h"
52 : #include "AliExternalTrackParam.h"
53 : #include "TROOT.h"
54 : #include "TFile.h"
55 : #include "TDirectory.h"
56 : #include "TGraph.h"
57 : #include "TMultiGraph.h"
58 : #include "TCanvas.h"
59 : #include "THnSparse.h"
60 : #include "TLegend.h"
61 : #include "TPad.h"
62 : #include "TH2D.h"
63 : #include "TH3D.h"
64 : #include "AliTPCROC.h"
65 : #include "AliTPCCalROC.h"
66 : //#include "AliESDfriend.h"
67 : #include "AliTPCcalibTime.h"
68 : #include "AliSplineFit.h"
69 : #include "AliCDBMetaData.h"
70 : #include "AliCDBId.h"
71 : #include "AliCDBManager.h"
72 : #include "AliCDBStorage.h"
73 : #include "AliTPCcalibBase.h"
74 : #include "AliTPCcalibDB.h"
75 : #include "AliTPCcalibDButil.h"
76 : #include "AliRelAlignerKalman.h"
77 : #include "AliTPCParamSR.h"
78 : #include "AliTPCcalibTimeGain.h"
79 : #include "AliTPCcalibGainMult.h"
80 : #include "AliTPCcalibAlign.h"
81 : #include "AliSplineFit.h"
82 : #include "AliTPCComposedCorrection.h"
83 : #include "AliTPCExBTwist.h"
84 : #include "AliTPCCalibGlobalMisalignment.h"
85 : #include "TStatToolkit.h"
86 : #include "TChain.h"
87 : #include "TCut.h"
88 : #include "AliTrackerBase.h"
89 : #include "AliTracker.h"
90 : #include "AliTPCPreprocessorOffline.h"
91 : #include "AliTPCCorrectionFit.h"
92 : #include "AliCDBEntry.h"
93 :
94 : #include "AliTPCClusterParam.h"
95 : #include "AliTPCRecoParam.h"
96 :
97 : using std::endl;
98 : using std::cout;
99 :
100 6 : ClassImp(AliTPCPreprocessorOffline)
101 : //_____________________________________________________________________________
102 : AliTPCPreprocessorOffline::AliTPCPreprocessorOffline():
103 0 : TNamed("TPCPreprocessorOffline","TPCPreprocessorOffline"),
104 0 : fNormaliseQA(kTRUE),
105 0 : fGainCalibrationType(kFullGainCalib), // gain calibration type
106 0 : fMinEntries(500), // minimal number of entries for fit
107 0 : fStartRun(0), // start Run - used to make fast selection in THnSparse
108 0 : fEndRun(0), // end Run - used to make fast selection in THnSparse
109 0 : fStartTime(0), // fStartTime - used to make fast selection in THnSparse
110 0 : fEndTime(0), // fEndTime - used to make fast selection in THnSparse
111 0 : fOCDBstorage(0), // OCDB storage
112 0 : fVdriftArray(new TObjArray),
113 0 : fTimeDrift(0),
114 0 : fGraphMIP(0), // graph time dependence of MIP
115 0 : fGraphCosmic(0), // graph time dependence at Plateu
116 0 : fGraphAttachmentMIP(0),
117 0 : fFitMIP(0), // fit of dependence - MIP
118 0 : fFitCosmic(0), // fit of dependence - Plateu
119 0 : fGainArray(new TObjArray), // array to be stored in the OCDB
120 0 : fGainArrayCombined(0x0),
121 0 : fArrQAhist(0x0),
122 0 : fGainMIP(0), // calibration component for MIP
123 0 : fGainCosmic(0), // calibration component for cosmic
124 0 : fGainMult(0),
125 0 : fAlignTree(0), // alignment tree
126 0 : fSwitchOnValidation(kFALSE), // flag to switch on validation of OCDB parameters
127 0 : fMinGain(1.5),
128 0 : fMaxGain(4.5),
129 0 : fMaxVdriftCorr(0.03),
130 0 : fNtracksVdrift(0),
131 0 : fMinTracksVdrift(0),
132 0 : fNeventsVdrift(0),
133 0 : fMinEventsVdrift(0),
134 0 : fCalibrationStatus(0),
135 0 : fDriftCDBentry(NULL)
136 0 : {
137 : //
138 : // default constructor
139 : //
140 0 : }
141 :
142 : //_____________________________________________________________________________
143 0 : AliTPCPreprocessorOffline::~AliTPCPreprocessorOffline() {
144 : //
145 : // Destructor
146 : //
147 0 : delete fDriftCDBentry;
148 0 : }
149 :
150 : //_____________________________________________________________________________
151 : void AliTPCPreprocessorOffline::GetRunRange(AliTPCcalibTime * const timeDrift){
152 : //
153 : // find the fist and last run
154 : //
155 0 : TObjArray *hisArray =timeDrift->GetHistoDrift();
156 0 : {for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
157 0 : THnSparse* addHist=(THnSparse*)hisArray->UncheckedAt(i);
158 0 : if (!addHist) continue;
159 0 : if (addHist->GetEntries()<fMinEntries) continue;
160 0 : TH1D* histo =addHist->Projection(3);
161 0 : TH1D* histoTime=addHist->Projection(0);
162 0 : AliInfo(Form("%s\t%f\t%d\t%d",histo->GetName(), histo->GetEntries(),histo->FindFirstBinAbove(0),histo->FindLastBinAbove(0)));
163 :
164 0 : if (fStartRun<=0){
165 0 : fStartRun=histo->FindFirstBinAbove(0);
166 0 : fEndRun =histo->FindLastBinAbove(0);
167 0 : }else{
168 0 : fStartRun=TMath::Min(histo->FindFirstBinAbove(0),fStartRun);
169 0 : fEndRun =TMath::Max(histo->FindLastBinAbove(0),fEndRun);
170 : }
171 0 : if (fStartTime==0){
172 0 : fStartTime=histoTime->FindFirstBinAbove(0);
173 0 : fEndTime =histoTime->FindLastBinAbove(0);
174 0 : }else{
175 0 : fStartTime=TMath::Min(histoTime->FindFirstBinAbove(0),fStartTime);
176 0 : fEndTime =TMath::Max(histoTime->FindLastBinAbove(0),fEndTime);
177 : }
178 0 : delete histo;
179 0 : delete histoTime;
180 0 : }}
181 0 : if (fStartRun<0) fStartRun=0;
182 0 : if (fEndRun<0) fEndRun=100000000;
183 0 : AliInfo(Form("Run range :\t%d-%d", fStartRun, fEndRun));
184 0 : AliInfo(Form("Time range :\t%d-%d", fStartTime, fEndTime));
185 :
186 0 : }
187 :
188 : //_____________________________________________________________________________
189 : Int_t AliTPCPreprocessorOffline::CalibTimeVdrift(AliTPCcalibTime* timeDrift, Int_t ustartRun, Int_t uendRun)
190 : {
191 : // make calibration of the drift velocity
192 : // Input parameters:
193 : // timeDrift - the calibration object
194 : // ustartRun, uendRun - run validity period
195 : // return 0 on success, 1 on failure
196 :
197 0 : fTimeDrift=timeDrift;
198 0 : fStartRun=ustartRun;
199 0 : fEndRun=ustartRun;
200 0 : GetRunRange(fTimeDrift);
201 :
202 : //extract statistics
203 0 : fNtracksVdrift = TMath::Nint(fTimeDrift->GetResHistoTPCITS(0)->GetEntries());
204 : //if we have 0 ITS TPC matches it means we have no ITS tracks and we try to use TPC-TOF matching for calibration
205 0 : if (fNtracksVdrift==0) fNtracksVdrift=TMath::Nint(fTimeDrift->GetResHistoTPCTOF(0)->GetEntries());
206 0 : fNeventsVdrift = TMath::Nint(fTimeDrift->GetTPCVertexHisto(0)->GetEntries());
207 :
208 0 : TObjArray *hisArray =fTimeDrift->GetHistoDrift();
209 0 : for (Int_t i=0; i<hisArray->GetEntriesFast(); i++){
210 0 : THnSparse* addHist=(THnSparse*)hisArray->At(i);
211 0 : if (!addHist) continue;
212 0 : if (fStartTime<fEndTime) addHist->GetAxis(0)->SetRange(fStartTime-1,fEndTime+1);
213 0 : if (fStartRun<fEndRun) addHist->GetAxis(3)->SetRange(fStartRun-1,fEndRun+1);
214 0 : }
215 : //
216 : //
217 : // 2. extraction of the information
218 : //
219 0 : if (fVdriftArray)
220 : {
221 0 : fVdriftArray->Delete();
222 0 : delete fVdriftArray;
223 : }
224 0 : fVdriftArray = new TObjArray();
225 0 : AddAlignmentGraphs(fVdriftArray,fTimeDrift);
226 0 : AddHistoGraphs(fVdriftArray,fTimeDrift,fMinEntries);
227 0 : AddLaserGraphs(fVdriftArray,fTimeDrift);
228 :
229 : //
230 : // 3. Append QA plots
231 : //
232 0 : MakeDefaultPlots(fVdriftArray,fVdriftArray);
233 :
234 : //
235 : // 4. validate OCDB entries
236 : //
237 0 : if(fSwitchOnValidation==kTRUE && ValidateTimeDrift()==kFALSE) {
238 0 : AliWarning("TPC time drift OCDB parameters out of range!");
239 0 : return 1;
240 : }
241 : //
242 : //4.b make alignment
243 : //
244 0 : MakeFitTime();
245 0 : TFile * ftime= TFile::Open("fitITSVertex.root");
246 0 : if (ftime){
247 0 : TObject * alignmentTime=ftime->Get("FitCorrectionTime");
248 0 : if (alignmentTime) fVdriftArray->AddLast(alignmentTime);
249 0 : }
250 : //
251 : // 5.) Add the RecoParam and ClusterParam - for compatibility checks -different sets of parameters can invalidate calibration
252 : //
253 0 : AliTPCClusterParam *clParam = AliTPCcalibDB::Instance()->GetClusterParam();
254 0 : TObjArray *recoParams = new TObjArray(4) ;
255 0 : for (Int_t i=0;i<4;i++) recoParams->AddAt(AliTPCcalibDB::Instance()->GetRecoParam(i),i);
256 0 : fVdriftArray->AddLast(clParam);
257 0 : fVdriftArray->AddLast(recoParams);
258 : return 0;
259 0 : }
260 :
261 : //_____________________________________________________________________________
262 : void AliTPCPreprocessorOffline::CalibTimeVdrift(const Char_t* file, Int_t ustartRun, Int_t uendRun, AliCDBStorage* pocdbStorage){
263 : //
264 : // make calibration of the drift velocity
265 : // Input parameters:
266 : // file - the location of input file
267 : // ustartRun, uendRun - run validity period
268 : // pocdbStorage - path to hte OCDB storage
269 : // - if empty - local storage 'pwd' uesed
270 0 : if (pocdbStorage) fOCDBstorage=pocdbStorage;
271 : else {
272 0 : TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
273 0 : fOCDBstorage=AliCDBManager::Instance()->GetStorage(localStorage.Data());
274 0 : }
275 :
276 : //
277 : // 1. Extract the calibration object form file, may have a number of layouts
278 0 : TFile fcalib(file);
279 0 : TObject* obj = dynamic_cast<TObject*>(fcalib.Get("TPCCalib"));
280 0 : TObjArray* array = dynamic_cast<TObjArray*>(obj);
281 0 : TDirectory* dir = dynamic_cast<TDirectory*>(obj);
282 : AliTPCcalibTime* timeDrift = NULL;
283 0 : if (dir) {
284 0 : timeDrift = dynamic_cast<AliTPCcalibTime*>(dir->Get("calibTime"));
285 0 : }
286 0 : else if (array){
287 0 : timeDrift = (AliTPCcalibTime *)array->FindObject("calibTime");
288 0 : } else {
289 0 : timeDrift = (AliTPCcalibTime*)fcalib.Get("calibTime");
290 : }
291 0 : if(!timeDrift) return;
292 :
293 : //calculate the calibration
294 0 : if (CalibTimeVdrift(timeDrift,ustartRun,uendRun)!=0) return;
295 :
296 : //
297 : //
298 : // 6. update of OCDB
299 : //
300 : //
301 0 : UpdateOCDBDrift(ustartRun,uendRun,fOCDBstorage);
302 0 : }
303 :
304 : //_____________________________________________________________________________
305 : AliCDBEntry* AliTPCPreprocessorOffline::CreateDriftCDBentryObject(Int_t ustartRun, Int_t uendRun)
306 : {
307 : //create the CDB entry (and cache internally)
308 :
309 : //reset if we dont have anything
310 0 : if (!fVdriftArray)
311 : {
312 0 : delete fDriftCDBentry; fDriftCDBentry=NULL;
313 0 : return NULL;
314 : }
315 :
316 : //will be owned by the cdb entry once attached
317 0 : AliCDBMetaData* metaData = new AliCDBMetaData;
318 0 : metaData->SetObjectClassName("TObjArray");
319 0 : metaData->SetResponsible("Marian Ivanov");
320 0 : metaData->SetBeamPeriod(1);
321 0 : metaData->SetAliRootVersion("05-25-01"); //root version
322 0 : metaData->SetComment("Calibration of the time dependence of the drift velocity");
323 :
324 0 : AliCDBId id1("TPC/Calib/TimeDrift", ustartRun, uendRun);
325 :
326 : //now the entry owns the metadata, but NOT the data
327 0 : delete fDriftCDBentry;
328 0 : fDriftCDBentry=new AliCDBEntry(fVdriftArray,id1,metaData,kFALSE);
329 :
330 : return fDriftCDBentry;
331 0 : }
332 :
333 : //_____________________________________________________________________________
334 : void AliTPCPreprocessorOffline::UpdateOCDBDrift( Int_t ustartRun, Int_t uendRun, AliCDBStorage* storage ){
335 : //
336 : // Update OCDB
337 : //
338 : Bool_t status=kFALSE;
339 0 : if (CreateDriftCDBentryObject(ustartRun,uendRun))
340 : {
341 0 : status=storage->Put(fDriftCDBentry);
342 0 : }
343 0 : if (status==kFALSE) fCalibrationStatus|=kCalibFailedExport ;
344 0 : }
345 :
346 : void AliTPCPreprocessorOffline::TakeOwnershipDriftCDBEntry(){
347 0 : fVdriftArray = NULL;
348 0 : fDriftCDBentry = NULL;
349 0 : }
350 :
351 : //_____________________________________________________________________________
352 : Bool_t AliTPCPreprocessorOffline::ValidateTimeGain()
353 : {
354 : //
355 : // Validate time gain corrections
356 : //
357 0 : AliInfo("ValidateTimeGain..." );
358 0 : Float_t minGain = fMinGain;
359 0 : Float_t maxGain = fMaxGain;
360 :
361 0 : TGraphErrors *gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL");
362 :
363 : // ===| treat the case if the combined calibration should be used |==========
364 0 : if (fGainCalibrationType==kCombinedGainCalib) {
365 0 : gr = (TGraphErrors*)fGainArrayCombined->FindObject("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL");
366 0 : }
367 :
368 0 : if (!gr) {
369 0 : gr = (TGraphErrors*)fGainArray->FindObject("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL");
370 0 : if (fGainCalibrationType==kCombinedGainCalib) {
371 0 : gr = (TGraphErrors*)fGainArrayCombined->FindObject("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL");
372 0 : }
373 0 : if (!gr)
374 : {
375 0 : fCalibrationStatus |= kCalibFailedTimeGain;
376 0 : return kFALSE;
377 : }
378 0 : AliInfo("Assuming given run is a cosmic run. Using gain calibration from Fermi-plateau muons.");
379 0 : }
380 0 : if(gr->GetN()<1)
381 : {
382 0 : fCalibrationStatus |= kCalibFailedTimeGain;
383 0 : return kFALSE;
384 : }
385 :
386 : // check whether gain in the range
387 0 : for(Int_t iPoint=0; iPoint<gr->GetN(); iPoint++)
388 : {
389 0 : if(gr->GetY()[iPoint] < minGain || gr->GetY()[iPoint] > maxGain)
390 : {
391 0 : fCalibrationStatus |= kCalibFailedTimeGain;
392 0 : AliError(TString::Format("Gain point outside of range %f != (%f,%f)",gr->GetY()[iPoint], minGain, maxGain ).Data());
393 0 : return kFALSE;
394 : }
395 : }
396 :
397 0 : AliInfo("Validation successful");
398 0 : return kTRUE;
399 0 : }
400 :
401 : //_____________________________________________________________________________
402 : AliTPCPreprocessorOffline::EGainCalibType AliTPCPreprocessorOffline::GetGainCalibrationTypeFromString(const TString& type)
403 : {
404 : //
405 : // return the gain calibration type analysing the string 'type'
406 : // if an error occurs, return kNGainCalibTypes
407 : //
408 0 : if (type.IsNull()) {
409 0 : ::Warning("AliTPCPreprocessorOffline::GetGainCalibrationTypeFromString","Empty gain calibration type string");
410 0 : return kNGainCalibTypes;
411 : }
412 :
413 0 : if (!type.IsDigit()) {
414 0 : ::Warning("AliTPCPreprocessorOffline::GetGainCalibrationTypeFromString","Gain calibration string '%s' is not a digit", type.Data());
415 0 : return kNGainCalibTypes;
416 : }
417 :
418 0 : const Int_t gainType = type.Atoi();
419 0 : if (gainType<0 || gainType>=Int_t(kNGainCalibTypes)) {
420 0 : ::Warning("AliTPCPreprocessorOffline::GetGainCalibrationTypeFromString","Gain calibration string '%s' is not a valid gain calibration type (0-%d)", type.Data(), kNGainCalibTypes-1);
421 0 : return kNGainCalibTypes;
422 : }
423 :
424 0 : return EGainCalibType(gainType);
425 0 : }
426 :
427 : //_____________________________________________________________________________
428 : Bool_t AliTPCPreprocessorOffline::SetGainCalibrationType(const TString& type)
429 : {
430 : //
431 : // set the gain calibration type analysing the string 'type'
432 : //
433 0 : EGainCalibType gainType=GetGainCalibrationTypeFromString(type);
434 0 : if (type==kNGainCalibTypes) return kFALSE;
435 :
436 0 : fGainCalibrationType=gainType;
437 :
438 0 : return kTRUE;
439 0 : }
440 :
441 : //_____________________________________________________________________________
442 : Bool_t AliTPCPreprocessorOffline::ProduceCombinedGainCalibration()
443 : {
444 : //
445 : // This function will produce the combined calibratin of the objects presently stored in the OCDB
446 : // and the calibration extracted in this calibration
447 : //
448 :
449 : // ===| write some info |=====================================================
450 0 : AliCDBManager *cdbMan = AliCDBManager::Instance();
451 0 : TString calibTimeGainStorage=cdbMan->GetDefaultStorage()->GetURI();
452 0 : const AliCDBEntry *e=cdbMan->Get("TPC/Calib/TimeGain");
453 0 : const TString timeGainID=e->GetId().ToString();
454 0 : if (cdbMan->GetSpecificStorage("TPC/Calib/TimeGain")) {
455 0 : calibTimeGainStorage=cdbMan->GetSpecificStorage("TPC/Calib/TimeGain")->GetURI();
456 : }
457 0 : AliInfoF("Using TimeGain '%s','%s' for combined gain calibration", calibTimeGainStorage.Data(), timeGainID.Data());
458 :
459 : // ===| get latest gain calibration from OCDB |===============================
460 0 : TObjArray *gainOCDB = AliTPCcalibDB::Instance()->GetTimeGainSplines();
461 0 : if (!gainOCDB) {
462 0 : AliError("Could not retrieve gain calibration from OCDB. Cannot perform combined calibration.");
463 0 : return kFALSE;
464 : }
465 :
466 0 : const Int_t nOCDB=gainOCDB->GetEntriesFast();
467 0 : const Int_t nThis=fGainArray->GetEntriesFast();
468 :
469 : // // ===| check consistency of entries |=======================================
470 : // TString type("this");
471 : // if ( nOCDB != nThis ) {
472 : // AliError(TString::Format("Entries in present OCDB calibration and this pass differ: %d != %d", nOCDB, nThis));
473 : // TObjArray *a1=gainOCDB->GetEntriesFast();
474 : // TObjArray *a2=fGainArray->GetEntriesFast();
475 : // if (nThis>nOCDB) {
476 : // a2=gainOCDB->GetEntriesFast();
477 : // a1=fGainArray->GetEntriesFast();
478 : // type="OCDB";
479 : // }
480 : //
481 : // for (Int_t icalib=0; icalib<a1->GetEntriesFast(); ++icalib) {
482 : // const TObject *o=a1->At(icalib);
483 : // if (!o) continue;
484 : // if (!a2->FindObject(o->GetName())) {
485 : // AliError(TString::Format("Could not find '%s' in %s calibration", o->GetName(), type.Data()));
486 : // }
487 : // }
488 : // return kFALSE;
489 : // }
490 :
491 : // ===| create combined gain array if needed and reset |=====================
492 0 : if (!fGainArrayCombined) fGainArrayCombined=new TObjArray(nThis);
493 : //fGainArrayCombined->SetOwner(); // either not owner, or all steering objects must be cloned
494 0 : fGainArrayCombined->Clear();
495 :
496 : // ===| explicitly treat entries |===========================================
497 0 : TGraphErrors *grOCDB = 0x0;
498 0 : TGraphErrors *grThis = 0x0;
499 : TGraphErrors *grCombined = 0x0;
500 : AliSplineFit *splineFit = 0x0;
501 :
502 : Bool_t error=kFALSE;
503 :
504 : // ---| gain vs time |--------------------------------------------------------
505 0 : GetGraphs("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL", grOCDB, grThis);
506 0 : AliInfo(Form("Graphs: %p, %p: %s, %s", grOCDB, grThis, grOCDB?grOCDB->GetName():"", grThis?grThis->GetName():""));
507 0 : if (!grOCDB || !grThis) {
508 0 : AliError("Gain vs. time for beam cannot be processed");
509 0 : if (fGainMIP) {
510 : error=kTRUE;
511 0 : }
512 : } else {
513 0 : grCombined = CombineGraphs(grOCDB, grThis, 1 );
514 0 : AliInfo(Form("Combined: %p: %s", grCombined, grCombined?grCombined->GetName():""));
515 0 : if (grCombined) {
516 0 : splineFit = AliTPCcalibTimeGain::MakeSplineFit(grCombined);
517 0 : fGainArrayCombined->AddAt(splineFit ,0);
518 0 : fGainArrayCombined->AddAt(grCombined,2);
519 : } else {
520 0 : AliError("Gain vs. time for beam cannot be processed");
521 : error=kTRUE;
522 : }
523 : }
524 :
525 0 : GetGraphs("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL", grOCDB, grThis);
526 0 : AliInfo(Form("Graphs: %p, %p: %s, %s", grOCDB, grThis, grOCDB?grOCDB->GetName():"", grThis?grThis->GetName():""));
527 0 : if (!grOCDB || !grThis) {
528 0 : AliError("Gain vs. time from cosmics cannot be processed");
529 0 : if (fGainCosmic) {
530 : error=kTRUE;
531 0 : }
532 : } else {
533 0 : grCombined = CombineGraphs(grOCDB, grThis, 1);
534 0 : AliInfo(Form("Combined: %p: %s", grCombined, grCombined?grCombined->GetName():""));
535 0 : if (grCombined) {
536 0 : splineFit = AliTPCcalibTimeGain::MakeSplineFit(grCombined);
537 0 : fGainArrayCombined->AddAt(splineFit ,1);
538 0 : fGainArrayCombined->AddAt(grCombined,3);
539 : } else {
540 0 : AliError("Gain vs. time from cosmics cannot be processed");
541 : error=kTRUE;
542 : }
543 :
544 : }
545 :
546 : // ===| steering objects |====================================================
547 0 : TObjArray steeringObjectNames;
548 0 : steeringObjectNames.Add(new TNamed("GainSlopesHV","1"));
549 0 : steeringObjectNames.Add(new TNamed("GainSlopesPT","1"));
550 0 : steeringObjectNames.Add(new TNamed("AliTPCClusterParam","1"));
551 0 : steeringObjectNames.Add(new TNamed("TObjArray","1"));
552 :
553 0 : for (Int_t isteer=0; isteer<steeringObjectNames.GetEntriesFast(); ++isteer) {
554 0 : TObject *objName = steeringObjectNames.At(isteer);
555 0 : TObject *steerObj = fGainArray->FindObject(objName->GetName());
556 0 : if (!steerObj) {
557 0 : AliErrorF("%s cannot be processed", objName->GetName());
558 : // ---| check if missing object should produce an error |---
559 0 : if (TString(objName->GetTitle()).Atoi()) {
560 : error=kTRUE;
561 0 : }
562 0 : continue;
563 : }
564 0 : fGainArrayCombined->AddLast(steerObj);
565 0 : }
566 :
567 : // ===| attachement is not applied, simply copy |=============================
568 0 : fGainArrayCombined->AddLast(fGainArray->FindObject("TGRAPHERRORS_MEAN_ATTACHMENT_BEAM_ALL"));
569 :
570 : // ===| multiplicative corrections |==========================================
571 0 : TObjArray graphsToCombine;
572 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL","1"));
573 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL","1"));
574 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL","0"));
575 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL","0"));
576 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEAN_CHAMBERGAIN_SHORT_BEAM_ALL","1"));
577 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEAN_CHAMBERGAIN_MEDIUM_BEAM_ALL","1"));
578 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_MEAN_CHAMBERGAIN_LONG_BEAM_ALL","1"));
579 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QMAX_DIPANGLE_SHORT_BEAM_ALL","1"));
580 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QTOT_DIPANGLE_SHORT_BEAM_ALL","1"));
581 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QMAX_DIPANGLE_MEDIUM_BEAM_ALL","1"));
582 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QTOT_DIPANGLE_MEDIUM_BEAM_ALL","1"));
583 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QMAX_DIPANGLE_LONG_BEAM_ALL","1"));
584 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QTOT_DIPANGLE_LONG_BEAM_ALL","1"));
585 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QMAX_DIPANGLE_ABSOLUTE_BEAM_ALL","1"));
586 0 : graphsToCombine.Add(new TNamed("TGRAPHERRORS_QTOT_DIPANGLE_ABSOLUTE_BEAM_ALL","1"));
587 :
588 0 : for (Int_t igraph=0; igraph<graphsToCombine.GetEntriesFast(); ++igraph) {
589 0 : TObject *graphObj=graphsToCombine.At(igraph);
590 0 : GetGraphs(graphObj->GetName(), grOCDB, grThis);
591 :
592 : // --- check graphs
593 0 : if (!grOCDB || !grThis) {
594 0 : AliError(TString::Format("%s cannot be processed missing object(s) OCDB: %d, This: %d", graphObj->GetName(), grOCDB!=0x0, grThis!=0x0));
595 0 : if (TString(graphObj->GetTitle()).Atoi()) {
596 : error=kTRUE;
597 0 : }
598 0 : continue;
599 : }
600 :
601 : // --- combine graphs
602 0 : grCombined = CombineGraphs(grOCDB, grThis);
603 0 : fGainArrayCombined->AddLast(grCombined);
604 :
605 : // --- normalize pad region gain
606 0 : if ( TString(graphObj->GetName()).Contains("PADREGIONGAIN") ) {
607 : // first normalize to the weighted mean using the number of
608 : // rows in the different pad regions
609 0 : NormaliseYToWeightedMeandEdx(grCombined);
610 : // now correct for the effect the truncated mean in the dE/dx
611 : // calculation would have
612 0 : NormaliseYToTruncateddEdx(grCombined);
613 : }
614 :
615 : // --- for dip angle graphs normalize and do fit
616 0 : if ( TString(graphObj->GetName()).Contains("DIPANGLE") ) {
617 0 : NormaliseYToMean(grCombined);
618 0 : TF1 * fun= new TF1("","1++abs(x)++abs(x*x)");
619 0 : grCombined->Fit(fun,"w","rob=0.9",-0.8,0.8);
620 0 : TString funName(graphObj->GetName());
621 0 : funName.ReplaceAll("TGRAPHERRORS","TF1");
622 0 : fun->SetNameTitle(funName,funName);
623 0 : fGainArrayCombined->AddLast(fun);
624 0 : }
625 0 : }
626 :
627 0 : return !error;
628 0 : }
629 :
630 : //_____________________________________________________________________________
631 : void AliTPCPreprocessorOffline::GetGraphs(const char* name, TGraphErrors* &grOCDB, TGraphErrors* &grThis)
632 : {
633 : //
634 : // get graphs from OCDB gain array and local gain array
635 : //
636 :
637 0 : grOCDB=0x0;
638 0 : grThis=0x0;
639 :
640 0 : TObjArray *gainOCDB = AliTPCcalibDB::Instance()->GetTimeGainSplines();
641 0 : if (!gainOCDB) return;
642 :
643 0 : grOCDB = dynamic_cast<TGraphErrors*>(gainOCDB ->FindObject(name));
644 0 : grThis = dynamic_cast<TGraphErrors*>(fGainArray->FindObject(name));
645 :
646 0 : if (!grOCDB) {
647 0 : AliError(TString::Format("Could not find graph '%s' in OCDB",name));
648 0 : }
649 :
650 0 : if (!grThis) {
651 0 : AliError(TString::Format("Could not find graph '%s' in present calibration",name));
652 0 : }
653 0 : }
654 :
655 : //_____________________________________________________________________________
656 : TGraphErrors* AliTPCPreprocessorOffline::CombineGraphs(TGraphErrors *grOCDB, TGraphErrors *grThis, const Int_t type/*=0*/, const Bool_t multiply/*=kTRUE*/)
657 : {
658 : //
659 : // Combine two graphs.
660 : // type
661 : // 0: Combine point by point, only do interpolation using Eval
662 : // 1: Combine using EvalConst
663 : // 2: Combine using Eval
664 : //
665 : // If mult is true, multiply the data points, otherwise add
666 : //
667 :
668 0 : Double_t x1,y1,x2,y2,x1Err,y1Err,y2Err,yNew,yNewErr;
669 0 : x1=y1=x2=y2=x1Err=y1Err=y2Err=yNew=yNewErr=0;
670 :
671 0 : const Int_t nOCDB=grOCDB->GetN();
672 0 : const Int_t nThis=grThis->GetN();
673 :
674 : // ===| output graph |========================================================
675 : // copy from gThis to retain the attributes
676 0 : TGraphErrors *grCombined=new TGraphErrors(*grThis);
677 0 : grCombined->Set(0);
678 :
679 : // ===| sort graphs for better comparison |===================================
680 0 : grOCDB->Sort();
681 0 : grThis->Sort();
682 :
683 0 : const Double_t xMinOCDB = grOCDB->GetX()[0];
684 0 : const Double_t xMaxOCDB = grOCDB->GetX()[nOCDB-1];
685 0 : const Double_t xMinThis = grThis->GetX()[0];
686 0 : const Double_t xMaxThis = grThis->GetX()[nThis-1];
687 :
688 : // ===========================================================================
689 : // ===| treat combination types |=============================================
690 : //
691 : const Double_t kVerySmall=1e-10;
692 : Int_t ninterpol=0;
693 :
694 0 : const Bool_t interpol = (type==1) || (type==2);
695 0 : const Bool_t evalConst = type!=2;
696 :
697 : Bool_t fallback=kFALSE;
698 :
699 0 : if (type==0) {
700 0 : if ( (xMinOCDB+kVerySmall < xMinThis) || (xMaxOCDB-kVerySmall > xMaxThis) ) {
701 0 : AliErrorF("Point by point combination of %s requested, but ranges are not compatible: [%.2f, %.2f] != [%.2f,%.2f]", grThis->GetName(), xMinOCDB, xMaxOCDB, xMinThis, xMaxThis);
702 0 : return 0x0;
703 : }
704 : }
705 :
706 0 : for (Int_t ipoint=0; ipoint<nThis; ++ipoint) {
707 0 : x2=y2=0;
708 0 : grThis->GetPoint(ipoint, x1, y1);
709 0 : x1Err=grThis->GetErrorX(ipoint);
710 0 : y1Err=grThis->GetErrorY(ipoint);
711 :
712 0 : if (type==0 && !fallback) {
713 0 : if (ipoint<nOCDB) {
714 0 : grOCDB->GetPoint(ipoint, x2, y2);
715 0 : y2Err=grOCDB->GetErrorY(ipoint);
716 0 : fallback=!TMath::AreEqualAbs(x1, x2, kVerySmall);
717 0 : }
718 : else {
719 : fallback=kTRUE;
720 : }
721 : }
722 :
723 0 : if (interpol || fallback) {
724 0 : GetPointWithError(grOCDB, x1, y2, y2Err, evalConst);
725 0 : ++ninterpol;
726 0 : }
727 :
728 0 : if (multiply) {
729 0 : yNew = y1*y2;
730 : yNewErr=0.;
731 0 : if (!TMath::AreEqualAbs(yNew, 0., kVerySmall)) {
732 0 : yNewErr = TMath::Sqrt( (y1Err*y1Err)/(y1*y1) + (y2Err*y2Err)/(y2*y2) ) * yNew;
733 0 : } else {
734 0 : AliErrorF("Cannot combined points from graph %s This y = %.4g +- %.4g, OCDB y = %.4g +- %.4g", grCombined->GetName(), y1, y1Err, y2, y2Err);
735 : }
736 : }
737 : else {
738 0 : yNew = y1+y2;
739 0 : yNewErr = TMath::Sqrt( (y1Err*y1Err) + (y2Err*y2Err) );
740 : }
741 :
742 0 : const Int_t point=grCombined->GetN();
743 0 : grCombined->SetPoint (point, x1 , yNew );
744 0 : grCombined->SetPointError(point, x1Err, yNewErr);
745 : }
746 :
747 0 : if ( (type==0) && ninterpol ) AliWarningF("%d number of points were interpolated for %s, although point-by-point combination was requested",ninterpol, grCombined->GetName());
748 :
749 0 : return grCombined;
750 0 : }
751 :
752 : //_____________________________________________________________________________
753 : Bool_t AliTPCPreprocessorOffline::GetPointWithError(const TGraphErrors *gr, const Double_t xPos, Double_t &y, Double_t &ey, Bool_t evalConst/*=kTRUE*/)
754 : {
755 : //
756 : // The function assumes the points to be sorted in x
757 : //
758 : // Evaluate the graph at xPos, do linear interpolation between neighbouring points
759 : // if 'evalConst' is true, the first/last point will be returned in case xPos is outside the graph range
760 : //
761 : // Errors in y are also calculated by linear interpolation between neighbouring points
762 : // Errors in x are not treated
763 : //
764 : // Return value indicates if xPos was inside the range
765 :
766 : // ===| reset input values |==================================================
767 0 : y=ey=0.;
768 :
769 : // ===| treat case of 0 point |===============================================
770 0 : if (!gr || gr->GetN()==0) return kFALSE;
771 :
772 : // ===| init variables |======================================================
773 0 : Double_t x1,x2,y1,y2, ey1,ey2;
774 0 : x1=x2=y1=y2=ey1=ey2=0.;
775 :
776 0 : const Int_t npoints=gr->GetN();
777 0 : const Double_t xmin=gr->GetX()[0];
778 0 : const Double_t xmax=gr->GetX()[npoints-1];
779 :
780 0 : const Bool_t belowBound = xPos<xmin;
781 0 : const Bool_t aboveBound = xPos>xmax;
782 0 : const Bool_t returnValue = !(belowBound || aboveBound);
783 :
784 : // ===| treat case of 1 point |===============================================
785 0 : if (npoints==1) {
786 0 : y = gr->GetY()[0];
787 0 : ey = gr->GetErrorY(0);
788 0 : return returnValue;
789 : }
790 :
791 : // ===| treat eval const case |===============================================
792 0 : if (evalConst) {
793 0 : if (belowBound) {
794 0 : y = gr->GetY()[0];
795 0 : ey = gr->GetErrorY(0);
796 0 : return returnValue;
797 : }
798 :
799 0 : if (aboveBound) {
800 0 : y = gr->GetY()[npoints-1];
801 0 : ey = gr->GetErrorY(npoints-1);
802 0 : return returnValue;
803 : }
804 : }
805 :
806 : // ===| 2 and more points |===================================================
807 0 : Int_t point = TMath::BinarySearch(npoints, gr->GetX(), xPos);
808 0 : Printf("n, i: %d, %d", npoints, point);
809 0 : if (point==-1) point=0;
810 0 : if (point==npoints-1) --point;
811 :
812 0 : gr->GetPoint(point, x1, y1);
813 0 : ey1 = gr->GetErrorY(point);
814 :
815 0 : gr->GetPoint(point+1, x2, y2);
816 0 : ey2 = gr->GetErrorY(point+1);
817 :
818 0 : Printf("%d, (%.2f, %.2f), (%.2f, %.2f)", point, x1, y1, x2, y2);
819 :
820 0 : if ( !(x2>x1) ) {
821 0 : AliTPCPreprocessorOffline p;
822 0 : p.Error("GetPointWithError","Graph not sorted, or error in extracting points");
823 : return kFALSE;
824 0 : }
825 :
826 0 : y = y1 + (y2-y1) /(x2-x1) * (xPos-x1);
827 0 : ey=ey1 + (ey2-ey1)/(x2-x1) * (xPos-x1);
828 :
829 : // ===| linear error increase outside bounds |================================
830 0 : if (belowBound) {
831 0 : ey=ey1 + (ey1)/(x2-x1) * TMath::Abs(xPos-x1);
832 0 : }
833 0 : else if (aboveBound) {
834 0 : ey=ey2 + (ey1)/(x2-x1) * TMath::Abs(xPos-x2);
835 0 : }
836 :
837 0 : return returnValue;
838 0 : }
839 :
840 : //_____________________________________________________________________________
841 : Bool_t AliTPCPreprocessorOffline::ValidateTimeDrift()
842 : {
843 : //
844 : // Validate time drift velocity corrections
845 : //
846 0 : AliInfo("ValidateTimeDrift..." );
847 :
848 0 : Float_t maxVDriftCorr = fMaxVdriftCorr;
849 :
850 0 : TGraphErrors* gr = (TGraphErrors*)fVdriftArray->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
851 0 : AliInfo(Form("ALIGN_ITSB_TPC_DRIFTVD graph = %p",gr));
852 0 : if (!gr)
853 : {
854 0 : gr = (TGraphErrors*)fVdriftArray->FindObject("ALIGN_TOFB_TPC_DRIFTVD");
855 0 : AliInfo(Form("ALIGN_TOFB_TPC_DRIFTVD graph = %p",gr));
856 0 : }
857 :
858 0 : if(!gr)
859 : {
860 0 : fCalibrationStatus|=kCalibFailedTimeDrift;
861 0 : return kFALSE;
862 : }
863 :
864 : // for now we validate even with low statistics
865 : ////check if we have enough statistics
866 : //if (fNtracksVdrift<fMinTracksVdrift)
867 : //{
868 : // fCalibrationStatus|=kCalibFailedTimeDrift;
869 : // return kFALSE;
870 : //}
871 :
872 0 : if(gr->GetN()<1) {
873 0 : AliInfo(Form("ALIGN_ITSB_TPC_DRIFTVD number of points = %d",gr->GetN()));
874 : {
875 0 : fCalibrationStatus|=kCalibFailedTimeDrift;
876 0 : return kFALSE;
877 : }
878 : }
879 :
880 : // check whether drift velocity corrections in the range
881 0 : for(Int_t iPoint = 0; iPoint<gr->GetN(); iPoint++)
882 : {
883 : //AliInfo(Form("Y value from the graph: %f",TMath::Abs(gr->GetY()[iPoint])));
884 0 : if(TMath::Abs(gr->GetY()[iPoint]) > maxVDriftCorr)
885 : {
886 0 : fCalibrationStatus|=kCalibFailedTimeDrift;
887 0 : return kFALSE;
888 : }
889 : }
890 :
891 0 : return kTRUE;
892 0 : }
893 :
894 : //_____________________________________________________________________________
895 : void AliTPCPreprocessorOffline::UpdateDriftParam(AliTPCParam *param, TObjArray *const arr, Int_t lstartRun){
896 : //
897 : // update the OCDB entry for the nominal time0
898 : //
899 : //
900 : // AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
901 0 : AliTPCParam *paramNew = (AliTPCParam *)param->Clone();
902 0 : TGraphErrors *grT = (TGraphErrors *)arr->FindObject("ALIGN_ITSM_TPC_T0");
903 0 : Double_t deltaTcm = TMath::Median(grT->GetN(),grT->GetY());
904 0 : Double_t deltaT = deltaTcm/param->GetDriftV();
905 0 : paramNew->SetL1Delay(param->GetL1Delay()-deltaT);
906 0 : paramNew->Update();
907 :
908 0 : AliCDBMetaData *metaData= new AliCDBMetaData();
909 0 : metaData->SetObjectClassName("TObjArray");
910 0 : metaData->SetResponsible("Marian Ivanov");
911 0 : metaData->SetBeamPeriod(1);
912 0 : metaData->SetAliRootVersion("05-25-02"); //root version
913 0 : metaData->SetComment("Updated calibration of nominal time 0");
914 : AliCDBId* id1=NULL;
915 0 : id1=new AliCDBId("TPC/Calib/Parameters", lstartRun, AliCDBRunRange::Infinity());
916 0 : Bool_t status = fOCDBstorage->Put(param, (*id1), metaData);
917 0 : if (status==kFALSE) fCalibrationStatus|=kCalibFailedExport ;
918 0 : }
919 :
920 : //_____________________________________________________________________________
921 : void AliTPCPreprocessorOffline::PrintArray(TObjArray *array){
922 : //
923 : // Print the names of the entries in array
924 : //
925 0 : Int_t entries = array->GetEntries();
926 0 : for (Int_t i=0; i<entries; i++){
927 0 : if (!array->At(i)) continue;
928 0 : Printf("%d\t %s", i, array->At(i)->GetName());
929 0 : }
930 0 : }
931 :
932 : //_____________________________________________________________________________
933 : TGraphErrors* AliTPCPreprocessorOffline::FilterGraphDrift(TGraphErrors * graph, Float_t errSigmaCut, Float_t medianCutAbs){
934 : // 2 filters:
935 : // 1. filter graph - error cut errSigmaCut
936 : // 2. filter graph - medianCutAbs around median
937 : //
938 : // errSigmaCut - cut on error
939 : // medianCutAbs - cut on value around median
940 0 : Double_t dummy=0; //
941 : //
942 : // 1. filter graph - error cut errSigmaCut
943 : //
944 : TGraphErrors *graphF;
945 0 : graphF = AliTPCcalibDButil::FilterGraphMedianErr(graph,errSigmaCut,dummy);
946 0 : delete graph;
947 0 : if (!graphF) return 0;
948 0 : graph = AliTPCcalibDButil::FilterGraphMedianErr(graphF,errSigmaCut,dummy);
949 0 : delete graphF;
950 0 : if (!graph) return 0;
951 : //
952 : // filter graph - kMedianCutAbs around median
953 : //
954 0 : graphF=FilterGraphMedianAbs(graph, medianCutAbs,dummy);
955 0 : delete graph;
956 0 : if (!graphF) return 0;
957 0 : graph=FilterGraphMedianAbs(graphF, medianCutAbs,dummy);
958 0 : delete graphF;
959 0 : if (!graph) return 0;
960 0 : return graph;
961 0 : }
962 :
963 : //_____________________________________________________________________________
964 : TGraphErrors* AliTPCPreprocessorOffline::FilterGraphMedianAbs(TGraphErrors * graph, Float_t cut,Double_t &medianY){
965 : //
966 : // filter outlyer measurement
967 : // Only points around median +- cut filtered
968 : //
969 0 : if (!graph) return 0;
970 : Int_t kMinPoints=2;
971 0 : Int_t npoints0 = graph->GetN();
972 : Int_t npoints=0;
973 : Float_t rmsY=0;
974 0 : Double_t *outx=new Double_t[npoints0];
975 0 : Double_t *outy=new Double_t[npoints0];
976 0 : Double_t *errx=new Double_t[npoints0];
977 0 : Double_t *erry=new Double_t[npoints0];
978 : //
979 : //
980 0 : if (npoints0<kMinPoints) {
981 0 : delete []outx;
982 0 : delete []outy;
983 0 : delete []errx;
984 0 : delete []erry;
985 0 : return 0;
986 : }
987 0 : for (Int_t iter=0; iter<3; iter++){
988 : npoints=0;
989 0 : for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
990 0 : if (graph->GetY()[ipoint]==0) continue;
991 0 : if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
992 0 : outx[npoints] = graph->GetX()[ipoint];
993 0 : outy[npoints] = graph->GetY()[ipoint];
994 0 : errx[npoints] = graph->GetErrorX(ipoint);
995 0 : erry[npoints] = graph->GetErrorY(ipoint);
996 0 : npoints++;
997 0 : }
998 0 : if (npoints<=1) break;
999 0 : medianY =TMath::Median(npoints,outy);
1000 0 : rmsY =TMath::RMS(npoints,outy);
1001 : }
1002 : TGraphErrors *graphOut=0;
1003 0 : if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
1004 0 : delete []outx;
1005 0 : delete []outy;
1006 0 : delete []errx;
1007 0 : delete []erry;
1008 : return graphOut;
1009 0 : }
1010 :
1011 : //_____________________________________________________________________________
1012 : void AliTPCPreprocessorOffline::AddHistoGraphs( TObjArray * vdriftArray, AliTPCcalibTime * const timeDrift, Int_t minEntries){
1013 : //
1014 : // Add graphs corresponding to the alignment
1015 : //
1016 : const Double_t kErrSigmaCut=5; // error sigma cut - for filtering
1017 : const Double_t kMedianCutAbs=0.03; // error sigma cut - for filtering
1018 : //
1019 0 : TObjArray * array=timeDrift->GetHistoDrift();
1020 0 : if (array){
1021 : THnSparse* hist=NULL;
1022 : // 2.a) cosmics with different triggers
1023 0 : for (Int_t i=0; i<array->GetEntriesFast();i++){
1024 0 : hist=(THnSparseF*)array->UncheckedAt(i);
1025 0 : if(!hist) continue;
1026 0 : if (hist->GetEntries()<minEntries) continue;
1027 : //hist->Print();
1028 0 : TString name=hist->GetName();
1029 0 : Int_t dim[4]={0,1,2,3};
1030 0 : THnSparse* newHist=hist->Projection(4,dim);
1031 0 : newHist->SetName(name);
1032 0 : TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
1033 0 : delete newHist;
1034 0 : if (!graph) {
1035 0 : AliInfo(Form("Graph =%s filtered out", name.Data()));
1036 0 : continue;
1037 : }
1038 0 : AliInfo(Form("name=%s graph=%i, N=%i", name.Data(), graph==0, graph->GetN()));
1039 0 : Int_t pos=name.Index("_");
1040 0 : name=name(pos,name.Capacity()-pos);
1041 0 : TString graphName=graph->ClassName();
1042 0 : graphName+=name;
1043 0 : graphName.ToUpper();
1044 : //
1045 0 : graph = FilterGraphDrift(graph, kErrSigmaCut, kMedianCutAbs);
1046 : //
1047 0 : if (graph){
1048 0 : graph->SetMarkerStyle(i%8+20);
1049 0 : graph->SetMarkerColor(i%7);
1050 0 : graph->GetXaxis()->SetTitle("Time");
1051 0 : graph->GetYaxis()->SetTitle("v_{dcor}");
1052 0 : graph->SetName(graphName);
1053 0 : graph->SetTitle(graphName);
1054 0 : AliInfo(Form("Graph %d\t=\t%s", i, graphName.Data()));
1055 0 : vdriftArray->Add(graph);
1056 : }
1057 0 : }
1058 0 : }
1059 0 : }
1060 :
1061 : //_____________________________________________________________________________
1062 : void AliTPCPreprocessorOffline::AddAlignmentGraphs( TObjArray * vdriftArray, AliTPCcalibTime *const timeDrift){
1063 : //
1064 : // Add graphs corresponding to alignment to the object array
1065 : //
1066 : TObjArray *arrayITS=0;
1067 : TObjArray *arrayTOF=0;
1068 : TObjArray *arrayTRD=0;
1069 : TMatrixD *mstatITS=0;
1070 : TMatrixD *mstatTOF=0;
1071 : TMatrixD *mstatTRD=0;
1072 : //
1073 0 : arrayITS=timeDrift->GetAlignITSTPC();
1074 0 : arrayTRD=timeDrift->GetAlignTRDTPC();
1075 0 : arrayTOF=timeDrift->GetAlignTOFTPC();
1076 :
1077 0 : if (arrayITS->GetEntries()>0) mstatITS= AliTPCcalibDButil::MakeStatRelKalman(arrayITS,0.7,50,fMaxVdriftCorr);
1078 0 : if (arrayTOF->GetEntries()>0) mstatTOF= AliTPCcalibDButil::MakeStatRelKalman(arrayTOF,0.7,1000,fMaxVdriftCorr);
1079 0 : if (arrayTRD->GetEntries()>0) mstatTRD= AliTPCcalibDButil::MakeStatRelKalman(arrayTRD,0.7,50,fMaxVdriftCorr);
1080 : //
1081 0 : TObjArray * arrayITSP= AliTPCcalibDButil::SmoothRelKalman(arrayITS,mstatITS, 0, 5.);
1082 0 : TObjArray * arrayITSM= AliTPCcalibDButil::SmoothRelKalman(arrayITS,mstatITS, 1, 5.);
1083 0 : TObjArray * arrayITSB= AliTPCcalibDButil::SmoothRelKalman(arrayITSP,arrayITSM);
1084 0 : TObjArray * arrayTOFP= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,mstatTOF, 0, 5.);
1085 0 : TObjArray * arrayTOFM= AliTPCcalibDButil::SmoothRelKalman(arrayTOF,mstatTOF, 1, 5.);
1086 0 : TObjArray * arrayTOFB= AliTPCcalibDButil::SmoothRelKalman(arrayTOFP,arrayTOFM);
1087 :
1088 : TObjArray * arrayTRDP= 0x0;
1089 : TObjArray * arrayTRDM= 0x0;
1090 : TObjArray * arrayTRDB= 0x0;
1091 0 : arrayTRDP= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,mstatTRD, 0, 5.);
1092 0 : arrayTRDM= AliTPCcalibDButil::SmoothRelKalman(arrayTRD,mstatTRD, 1, 5.);
1093 0 : arrayTRDB= AliTPCcalibDButil::SmoothRelKalman(arrayTRDP,arrayTRDM);
1094 : //
1095 : //
1096 0 : Int_t entries=TMath::Max(arrayITS->GetEntriesFast(),arrayTOF->GetEntriesFast());
1097 0 : TObjArray *arrays[12]={arrayITS, arrayITSP, arrayITSM, arrayITSB,
1098 : arrayTRD, arrayTRDP, arrayTRDM, arrayTRDB,
1099 : arrayTOF, arrayTOFP, arrayTOFM, arrayTOFB};
1100 0 : TString grnames[12]={"ALIGN_ITS", "ALIGN_ITSP", "ALIGN_ITSM", "ALIGN_ITSB",
1101 0 : "ALIGN_TRD", "ALIGN_TRDP", "ALIGN_TRDM","ALIGN_TRDB",
1102 0 : "ALIGN_TOF", "ALIGN_TOFP", "ALIGN_TOFM","ALIGN_TOFB"};
1103 0 : TString grpar[9]={"DELTAPSI", "DELTATHETA", "DELTAPHI",
1104 0 : "DELTAX", "DELTAY", "DELTAZ",
1105 0 : "DRIFTVD", "T0", "VDGY"};
1106 :
1107 :
1108 0 : TVectorD vX(entries);
1109 0 : TVectorD vY(entries);
1110 0 : TVectorD vEx(entries);
1111 0 : TVectorD vEy(entries);
1112 : TObjArray *arr=0;
1113 0 : for (Int_t iarray=0; iarray<12; iarray++){
1114 0 : arr = arrays[iarray];
1115 0 : if (arr==0) continue;
1116 0 : for (Int_t ipar=0; ipar<9; ipar++){
1117 : Int_t counter=0;
1118 0 : for (Int_t itime=0; itime<arr->GetEntriesFast(); itime++){
1119 0 : AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) arr->UncheckedAt(itime);
1120 0 : if (!kalman) continue;
1121 0 : vX[counter]=kalman->GetTimeStamp();
1122 0 : vY[counter]=(*(kalman->GetState()))[ipar];
1123 0 : if (ipar==6) vY[counter]=1./(*(kalman->GetState()))[ipar]-1;
1124 0 : vEx[counter]=0;
1125 0 : vEy[counter]=TMath::Sqrt((*(kalman->GetStateCov()))(ipar,ipar));
1126 0 : counter++;
1127 0 : }
1128 :
1129 0 : TGraphErrors * graph=new TGraphErrors(counter, vX.GetMatrixArray(),
1130 0 : vY.GetMatrixArray(),
1131 0 : vEx.GetMatrixArray(),
1132 0 : vEy.GetMatrixArray());
1133 0 : TString grName=grnames[iarray];
1134 0 : grName+="_TPC_";
1135 0 : grName+=grpar[ipar];
1136 0 : graph->SetName(grName.Data());
1137 0 : vdriftArray->AddLast(graph);
1138 0 : }
1139 0 : delete arrays[iarray];
1140 : }
1141 0 : }
1142 :
1143 : //_____________________________________________________________________________
1144 : void AliTPCPreprocessorOffline::AddLaserGraphs( TObjArray * vdriftArray, AliTPCcalibTime *timeDrift){
1145 : //
1146 : // add graphs for laser
1147 : //
1148 : const Double_t delayL0L1 = 0.071; //this is hack for 1/2 weeks
1149 : //THnSparse *hisN=0;
1150 0 : TGraphErrors *grLaser[6]={0,0,0,0,0,0};
1151 : //hisN = timeDrift->GetHistVdriftLaserA(0);
1152 0 : if (timeDrift->GetHistVdriftLaserA(0)){
1153 0 : grLaser[0]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserA(0),0,2,5,delayL0L1);
1154 0 : grLaser[0]->SetName("GRAPH_MEAN_DELAY_LASER_ALL_A");
1155 0 : vdriftArray->AddLast(grLaser[0]);
1156 0 : }
1157 0 : if (timeDrift->GetHistVdriftLaserA(1)){
1158 0 : grLaser[1]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserA(1),0,2,5);
1159 0 : grLaser[1]->SetName("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1160 0 : vdriftArray->AddLast(grLaser[1]);
1161 0 : }
1162 0 : if (timeDrift->GetHistVdriftLaserA(2)){
1163 0 : grLaser[2]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserA(2),0,2,5);
1164 0 : grLaser[2]->SetName("GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_A");
1165 0 : vdriftArray->AddLast(grLaser[2]);
1166 0 : }
1167 0 : if (timeDrift->GetHistVdriftLaserC(0)){
1168 0 : grLaser[3]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserC(0),0,2,5,delayL0L1);
1169 0 : grLaser[3]->SetName("GRAPH_MEAN_DELAY_LASER_ALL_C");
1170 0 : vdriftArray->AddLast(grLaser[3]);
1171 0 : }
1172 0 : if (timeDrift->GetHistVdriftLaserC(1)){
1173 0 : grLaser[4]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserC(1),0,2,5);
1174 0 : grLaser[4]->SetName("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1175 0 : vdriftArray->AddLast(grLaser[4]);
1176 0 : }
1177 0 : if (timeDrift->GetHistVdriftLaserC(2)){
1178 0 : grLaser[5]=MakeGraphFilter0(timeDrift->GetHistVdriftLaserC(2),0,2,5);
1179 0 : grLaser[5]->SetName("GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_C");
1180 0 : vdriftArray->AddLast(grLaser[5]);
1181 0 : }
1182 0 : for (Int_t i=0; i<6;i++){
1183 0 : if (grLaser[i]) {
1184 0 : SetDefaultGraphDrift(grLaser[i], 1,(i+20));
1185 0 : grLaser[i]->GetYaxis()->SetTitle("Laser Correction");
1186 0 : }
1187 : }
1188 0 : }
1189 :
1190 : //_____________________________________________________________________________
1191 : TGraphErrors * AliTPCPreprocessorOffline::MakeGraphFilter0(THnSparse *hisN, Int_t itime, Int_t ival, Int_t minEntries, Double_t offset){
1192 : //
1193 : // Make graph with mean values and rms
1194 : //
1195 0 : hisN->GetAxis(itime)->SetRange(0,100000000);
1196 0 : hisN->GetAxis(ival)->SetRange(0,100000000);
1197 0 : TH1 * hisT = hisN->Projection(itime);
1198 0 : TH1 * hisV = hisN->Projection(ival);
1199 : //
1200 0 : Int_t firstBinA = hisT->FindFirstBinAbove(2);
1201 0 : Int_t lastBinA = hisT->FindLastBinAbove(2);
1202 0 : Int_t firstBinV = hisV->FindFirstBinAbove(0);
1203 0 : Int_t lastBinV = hisV->FindLastBinAbove(0);
1204 0 : hisN->GetAxis(itime)->SetRange(firstBinA,lastBinA);
1205 0 : hisN->GetAxis(ival)->SetRange(firstBinV,lastBinV);
1206 : Int_t entries=0;
1207 0 : for (Int_t ibin=firstBinA; ibin<=lastBinA; ibin++){
1208 0 : Double_t cont = hisT->GetBinContent(ibin);
1209 0 : if (cont<minEntries) continue;
1210 0 : entries++;
1211 0 : }
1212 0 : TVectorD vecTime(entries);
1213 0 : TVectorD vecMean0(entries);
1214 0 : TVectorD vecRMS0(entries);
1215 0 : TVectorD vecMean1(entries);
1216 0 : TVectorD vecRMS1(entries);
1217 : entries=0;
1218 0 : for (Int_t ibin=firstBinA; ibin<=lastBinA; ibin++){
1219 0 : Double_t cont = hisT->GetBinContent(ibin);
1220 0 : if (cont<minEntries) continue;
1221 : //hisN->GetAxis(itime)->SetRange(ibin-1,ibin+1);
1222 0 : Int_t minBin = ibin-1;
1223 0 : Int_t maxBin = ibin+1;
1224 0 : if(minBin <= 0) minBin = 1;
1225 0 : if(maxBin >= hisN->GetAxis(itime)->GetNbins()) maxBin = hisN->GetAxis(itime)->GetNbins()-1;
1226 0 : hisN->GetAxis(itime)->SetRange(minBin,maxBin);
1227 :
1228 0 : Double_t time = hisT->GetBinCenter(ibin);
1229 0 : TH1 * his = hisN->Projection(ival);
1230 0 : Double_t nentries0= his->GetBinContent(his->FindBin(0));
1231 0 : if (cont-nentries0<minEntries) continue;
1232 : //
1233 0 : his->SetBinContent(his->FindBin(0),0);
1234 0 : vecTime[entries]=time;
1235 0 : vecMean0[entries]=his->GetMean()+offset;
1236 0 : vecMean1[entries]=his->GetMeanError();
1237 0 : vecRMS0[entries] =his->GetRMS();
1238 0 : vecRMS1[entries] =his->GetRMSError();
1239 0 : delete his;
1240 0 : entries++;
1241 0 : }
1242 0 : delete hisT;
1243 0 : delete hisV;
1244 0 : TGraphErrors * graph = new TGraphErrors(entries,vecTime.GetMatrixArray(), vecMean0.GetMatrixArray(), 0, vecMean1.GetMatrixArray());
1245 :
1246 : return graph;
1247 0 : }
1248 :
1249 : //_____________________________________________________________________________
1250 : void AliTPCPreprocessorOffline::SetDefaultGraphDrift(TGraph *graph, Int_t color, Int_t style){
1251 : //
1252 : // Set default style for QA views
1253 : //
1254 0 : graph->GetXaxis()->SetTimeDisplay(kTRUE);
1255 0 : graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
1256 0 : graph->SetMaximum( 0.025);
1257 0 : graph->SetMinimum(-0.025);
1258 0 : graph->GetXaxis()->SetTitle("Time");
1259 0 : graph->GetYaxis()->SetTitle("v_{dcorr}");
1260 : //
1261 0 : graph->GetYaxis()->SetLabelSize(0.03);
1262 0 : graph->GetXaxis()->SetLabelSize(0.03);
1263 : //
1264 0 : graph->GetXaxis()->SetNdivisions(10,5,0);
1265 0 : graph->GetYaxis()->SetNdivisions(10,5,0);
1266 : //
1267 0 : graph->GetXaxis()->SetLabelOffset(0.02);
1268 0 : graph->GetYaxis()->SetLabelOffset(0.005);
1269 : //
1270 0 : graph->GetXaxis()->SetTitleOffset(1.3);
1271 0 : graph->GetYaxis()->SetTitleOffset(1.2);
1272 : //
1273 0 : graph->SetMarkerColor(color);
1274 0 : graph->SetLineColor(color);
1275 0 : graph->SetMarkerStyle(style);
1276 0 : }
1277 :
1278 : //_____________________________________________________________________________
1279 : void AliTPCPreprocessorOffline::SetPadStyle(TPad *pad, Float_t mx0, Float_t mx1, Float_t my0, Float_t my1){
1280 : //
1281 : // Set default pad style for QA
1282 : //
1283 0 : pad->SetTicks(1,1);
1284 0 : pad->SetMargin(mx0,mx1,my0,my1);
1285 0 : }
1286 :
1287 : //_____________________________________________________________________________
1288 : void AliTPCPreprocessorOffline::MakeDefaultPlots(TObjArray * const arr, TObjArray * /*picArray*/){
1289 : //
1290 : // 0. make a default QA plots
1291 : // 1. Store them in the array
1292 : //
1293 : //
1294 : Float_t mx0=0.12, mx1=0.1, my0=0.15, my1=0.1;
1295 : //
1296 0 : TGraphErrors* laserA =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1297 0 : TGraphErrors* laserC =(TGraphErrors*)arr->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1298 0 : TGraphErrors* cosmic =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1299 0 : TGraphErrors* cross =(TGraphErrors*)arr->FindObject("TGRAPHERRORS_VDRIFT_CROSS_ALL");
1300 0 : TGraphErrors* itstpcP =(TGraphErrors*)arr->FindObject("ALIGN_ITSP_TPC_DRIFTVD");
1301 0 : TGraphErrors* itstpcM =(TGraphErrors*)arr->FindObject("ALIGN_ITSM_TPC_DRIFTVD");
1302 0 : TGraphErrors* itstpcB =(TGraphErrors*)arr->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
1303 : //
1304 0 : if (laserA) SetDefaultGraphDrift(laserA,2,25);
1305 0 : if (laserC) SetDefaultGraphDrift(laserC,4,26);
1306 0 : if (cosmic) SetDefaultGraphDrift(cosmic,3,27);
1307 0 : if (cross) SetDefaultGraphDrift(cross,4,28);
1308 0 : if (itstpcP) SetDefaultGraphDrift(itstpcP,2,29);
1309 0 : if (itstpcM) SetDefaultGraphDrift(itstpcM,4,30);
1310 0 : if (itstpcB) SetDefaultGraphDrift(itstpcB,1,31);
1311 : //
1312 : //
1313 : TPad *pad=0;
1314 : //
1315 : // Laser-Laser
1316 : //
1317 0 : if (laserA&&laserC){
1318 0 : pad = new TCanvas("TPCLaserVDrift","TPCLaserVDrift");
1319 0 : laserA->Draw("alp");
1320 0 : SetPadStyle(pad,mx0,mx1,my0,my1);
1321 0 : laserA->Draw("apl");
1322 0 : laserC->Draw("p");
1323 0 : TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1324 0 : legend->AddEntry(laserA,"Laser A side");
1325 0 : legend->AddEntry(laserC,"Laser C side");
1326 0 : legend->Draw();
1327 : //picArray->AddLast(pad);
1328 0 : }
1329 :
1330 0 : if (itstpcP&&itstpcM&&itstpcB){
1331 0 : pad = new TCanvas("ITSTPC","ITSTPC");
1332 0 : itstpcP->Draw("alp");
1333 0 : SetPadStyle(pad,mx0,mx1,my0,my1);
1334 0 : itstpcP->Draw("alp");
1335 0 : gPad->Clear();
1336 0 : itstpcM->Draw("apl");
1337 0 : itstpcP->Draw("p");
1338 0 : itstpcB->Draw("p");
1339 0 : TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1340 0 : legend->AddEntry(itstpcP,"ITS-TPC smooth plus");
1341 0 : legend->AddEntry(itstpcM,"ITS-TPC smooth minus");
1342 0 : legend->AddEntry(itstpcB,"ITS-TPC smooth ");
1343 0 : legend->Draw();
1344 : //picArray->AddLast(pad);
1345 0 : }
1346 :
1347 0 : if (itstpcB&&laserA&&itstpcP&&itstpcM){
1348 0 : pad = new TCanvas("ITSTPC_LASER","ITSTPC_LASER");
1349 0 : SetPadStyle(pad,mx0,mx1,my0,my1);
1350 0 : laserA->Draw("alp");
1351 0 : itstpcP->Draw("p");
1352 0 : itstpcM->Draw("p");
1353 0 : itstpcB->Draw("p");
1354 0 : TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1355 0 : legend->AddEntry(laserA,"TPC laser");
1356 0 : legend->AddEntry(itstpcP,"ITS-TPC smooth plus");
1357 0 : legend->AddEntry(itstpcM,"ITS-TPC smooth minus");
1358 0 : legend->AddEntry(itstpcB,"ITS-TPC smooth ");
1359 0 : legend->Draw();
1360 : //picArray->AddLast(pad);
1361 0 : }
1362 :
1363 0 : if (itstpcP&&cross){
1364 0 : pad = new TCanvas("ITSTPC_CROSS","ITSTPC_CROSS");
1365 0 : SetPadStyle(pad,mx0,mx1,my0,my1);
1366 0 : itstpcP->Draw("alp");
1367 0 : pad->Clear();
1368 0 : cross->Draw("ap");
1369 0 : itstpcP->Draw("p");
1370 : //
1371 0 : TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1372 :
1373 0 : legend->AddEntry(cross,"TPC cross tracks");
1374 0 : legend->AddEntry(itstpcB,"ITS-TPC smooth");
1375 0 : legend->Draw();
1376 : //picArray->AddLast(pad);
1377 0 : }
1378 0 : if (itstpcP&&cosmic){
1379 0 : pad = new TCanvas("ITSTPC_COSMIC","ITSTPC_COSMIC");
1380 0 : SetPadStyle(pad,mx0,mx1,my0,my1);
1381 0 : itstpcP->Draw("alp");
1382 0 : pad->Clear();
1383 0 : cosmic->Draw("ap");
1384 0 : itstpcP->Draw("p");
1385 : //
1386 0 : TLegend *legend = new TLegend(mx0+0.01,1-my1-0.2, 0.5, 1-my1-0.01, "Drift velocity correction");
1387 :
1388 0 : legend->AddEntry(cosmic,"TPC cross tracks0 up-down");
1389 0 : legend->AddEntry(itstpcB,"ITS-TPC smooth");
1390 0 : legend->Draw();
1391 : //picArray->AddLast(pad);
1392 0 : }
1393 0 : }
1394 :
1395 : //_____________________________________________________________________________
1396 : void AliTPCPreprocessorOffline::CalibTimeGain(const Char_t* fileName, Int_t startRunNumber, Int_t endRunNumber, AliCDBStorage* fullStorage, AliCDBStorage* residualStorage){
1397 : //
1398 : // Update OCDB gain
1399 : // fullStorage is where the full calibration object should go
1400 : // residualStorage is where the residual calibration for QA purposes should go
1401 : //
1402 0 : if (fullStorage==0) {
1403 0 : const TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1404 0 : fullStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
1405 0 : }
1406 0 : if ( (fGainCalibrationType==kResidualGainQA || fGainCalibrationType==kCombinedGainCalib) && residualStorage==0x0) {
1407 0 : const TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/residualOCDB";
1408 0 : residualStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
1409 0 : }
1410 :
1411 : // dump info
1412 0 : AliInfoF("Analyse gain from file %s", fileName);
1413 :
1414 : //
1415 : // 1. Read gain values
1416 : //
1417 0 : ReadGainGlobal(fileName);
1418 :
1419 : //
1420 : // 2. Extract calibration values
1421 : //
1422 0 : AnalyzeGain(startRunNumber,endRunNumber, 1000,1.43);
1423 0 : AnalyzeAttachment(startRunNumber,endRunNumber);
1424 0 : AnalyzePadRegionGain();
1425 0 : AnalyzeGainMultiplicity();
1426 0 : AnalyzeGainChamberByChamber();
1427 : //
1428 0 : AnalyzeGainDipAngle(0); // short pads
1429 0 : AnalyzeGainDipAngle(1); // medium pads
1430 0 : AnalyzeGainDipAngle(2); // long pads
1431 0 : AnalyzeGainDipAngle(3); // absolute calibration on full track
1432 :
1433 : //
1434 : // 2.a produce combined calibration if requested
1435 : //
1436 :
1437 : Bool_t combinedGainSuccessful=kTRUE;
1438 0 : if (fGainCalibrationType==kCombinedGainCalib) {
1439 0 : combinedGainSuccessful=ProduceCombinedGainCalibration();
1440 0 : AliInfoF("Result of combined gain calibration: %d", combinedGainSuccessful);
1441 0 : }
1442 :
1443 : //
1444 : // 3. Make control plots
1445 : //
1446 0 : MakeQAPlot(1.43);
1447 :
1448 : //
1449 : // 4. validate OCDB entries
1450 : //
1451 0 : if(fSwitchOnValidation==kTRUE &&
1452 0 : (fGainCalibrationType==kFullGainCalib || fGainCalibrationType==kCombinedGainCalib) &&
1453 0 : (!combinedGainSuccessful || !ValidateTimeGain()) ) {
1454 0 : AliWarning("TPC time gain OCDB parameters out of range!");
1455 0 : return;
1456 : }
1457 :
1458 : //
1459 : // 5. Update OCDB
1460 : //
1461 0 : if (fGainCalibrationType != kNoGainCalib ) {
1462 0 : UpdateOCDBGain( startRunNumber, endRunNumber, fullStorage, residualStorage);
1463 0 : }
1464 0 : }
1465 :
1466 : //_____________________________________________________________________________
1467 : void AliTPCPreprocessorOffline::ReadGainGlobal(const Char_t* fileName){
1468 : //
1469 : // read calibration entries from file
1470 : //
1471 0 : TFile fcalib(fileName);
1472 0 : TObject* obj = dynamic_cast<TObject*>(fcalib.Get("TPCCalib"));
1473 0 : TObjArray * array = dynamic_cast<TObjArray*>(obj);
1474 0 : TDirectory * dir = dynamic_cast<TDirectory*>(obj);
1475 0 : if (dir) {
1476 0 : fGainMIP = dynamic_cast<AliTPCcalibTimeGain *>(dir->Get("calibTimeGain"));
1477 0 : fGainCosmic = dynamic_cast<AliTPCcalibTimeGain *>(dir->Get("calibTimeGainCosmic"));
1478 0 : fGainMult = dynamic_cast<AliTPCcalibGainMult *>(dir->Get("calibGainMult"));
1479 0 : }
1480 0 : else if (array){
1481 0 : fGainMIP = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
1482 0 : fGainCosmic = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGainCosmic");
1483 0 : fGainMult = ( AliTPCcalibGainMult *)array->FindObject("calibGainMult");
1484 0 : }else{
1485 0 : fGainMIP = ( AliTPCcalibTimeGain *)fcalib.Get("calibTimeGain");
1486 0 : fGainCosmic = ( AliTPCcalibTimeGain *)fcalib.Get("calibTimeGainCosmic");
1487 0 : fGainMult = ( AliTPCcalibGainMult *)fcalib.Get("calibGainMult");
1488 : }
1489 0 : if (!fGainMult){
1490 0 : TFile calibMultFile("TPCMultObjects.root");
1491 0 : fGainMult = ( AliTPCcalibGainMult *)calibMultFile.Get("calibGainMult");
1492 0 : }
1493 : TH1 * hisT=0;
1494 : Int_t firstBinA =0, lastBinA=0;
1495 :
1496 0 : if (fGainCosmic){
1497 0 : hisT= fGainCosmic->GetHistGainTime()->Projection(1);
1498 0 : firstBinA = hisT->FindFirstBinAbove(2);
1499 0 : lastBinA = hisT->FindLastBinAbove(2);
1500 0 : fGainCosmic->GetHistGainTime()->GetAxis(1)->SetRange(firstBinA,lastBinA);
1501 0 : delete hisT;
1502 : }
1503 :
1504 0 : if (fGainMIP){
1505 0 : hisT= fGainMIP->GetHistGainTime()->Projection(1);
1506 0 : firstBinA = hisT->FindFirstBinAbove(2);
1507 0 : lastBinA = hisT->FindLastBinAbove(2);
1508 0 : fGainMIP->GetHistGainTime()->GetAxis(1)->SetRange(firstBinA,lastBinA);
1509 0 : delete hisT;
1510 : }
1511 :
1512 0 : }
1513 :
1514 : //_____________________________________________________________________________
1515 : Bool_t AliTPCPreprocessorOffline::AnalyzeGain(Int_t startRunNumber, Int_t endRunNumber, Int_t minEntriesGaussFit, Float_t FPtoMIPratio){
1516 : //
1517 : // Analyze gain - produce the calibration graphs
1518 : // Uses the MIP pions integrated over all chambers, eta and multiplicity
1519 : // This provides the absolute normalization to channel 50 in the end
1520 : //
1521 : // TODO: o What happens in case the MIPs are not flat in eta, and/or multiplicity?
1522 : // In CPass0 it will not be flat.
1523 : // Can we fill the histograms flat in multiplicity?
1524 : // o Does the phi distribution play a role (inactive sector parts)?
1525 : // o In principle one could normalise over eta, by first projecting in slices
1526 : // of eta, normalize each slice and then sum
1527 : //
1528 :
1529 : // 1.) try to create MIP spline
1530 0 : if (fGainMIP)
1531 : {
1532 0 : fGainMIP->GetHistGainTime()->GetAxis(5)->SetRangeUser(startRunNumber, endRunNumber);
1533 0 : fGainMIP->GetHistGainTime()->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
1534 0 : fGainMIP->GetHistGainTime()->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
1535 :
1536 : // ---| QA histogram |---
1537 0 : if (fArrQAhist) {
1538 0 : TH2D *hQA = fGainMIP->GetHistGainTime()->Projection(0,1);
1539 0 : hQA->SetName("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL_QA");
1540 0 : hQA->SetTitle("MIP calibration collisions; time;d#it{E}/d#it{x} (arb. unit)");
1541 0 : hQA->GetXaxis()->SetRangeUser(hQA->FindFirstBinAbove(1), hQA->FindLastBinAbove(1));
1542 0 : fArrQAhist->Add(hQA);
1543 0 : }
1544 :
1545 : //
1546 0 : fGraphMIP = AliTPCcalibBase::FitSlices(fGainMIP->GetHistGainTime(),0,1,minEntriesGaussFit,10,0.1,0.7);
1547 0 : if (fGraphMIP->GetN()==0) fGraphMIP = 0x0;
1548 0 : if (fGraphMIP) fFitMIP = AliTPCcalibTimeGain::MakeSplineFit(fGraphMIP);
1549 0 : if (fGraphMIP) fGraphMIP->SetName("TGRAPHERRORS_MEAN_GAIN_BEAM_ALL");// set proper names according to naming convention
1550 0 : fGainArray->AddAt(fFitMIP,0);
1551 0 : }
1552 :
1553 : // 2.) try to create Cosmic spline
1554 0 : if (fGainCosmic)
1555 : {
1556 0 : fGainCosmic->GetHistGainTime()->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
1557 0 : fGainCosmic->GetHistGainTime()->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
1558 :
1559 : // ---| QA histogram |---
1560 0 : if (fArrQAhist) {
1561 0 : TH2D *hQA = fGainMIP->GetHistGainTime()->Projection(0,1);
1562 0 : hQA->SetName("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL_QA");
1563 0 : hQA->SetTitle("MIP calibration cosmics; time;d#it{E}/d#it{x} (arb. unit)");
1564 0 : hQA->GetXaxis()->SetRangeUser(hQA->FindFirstBinAbove(1), hQA->FindLastBinAbove(1));
1565 0 : fArrQAhist->Add(hQA);
1566 0 : }
1567 :
1568 : //
1569 0 : fGraphCosmic = AliTPCcalibBase::FitSlices(fGainCosmic->GetHistGainTime(),0,1,minEntriesGaussFit,10);
1570 0 : if (fGraphCosmic->GetN()==0) fGraphCosmic = 0x0;
1571 : //
1572 0 : if (fGraphCosmic) {
1573 0 : for(Int_t i=0; i < fGraphCosmic->GetN(); i++) {
1574 0 : fGraphCosmic->GetY()[i] /= FPtoMIPratio;
1575 0 : fGraphCosmic->GetEY()[i] /= FPtoMIPratio;
1576 : }
1577 0 : }
1578 : //
1579 0 : if (fGraphCosmic) fFitCosmic = AliTPCcalibTimeGain::MakeSplineFit(fGraphCosmic);
1580 0 : if (fGraphCosmic) fGraphCosmic->SetName("TGRAPHERRORS_MEAN_GAIN_COSMIC_ALL"); // set proper names according to naming convention
1581 0 : fGainArray->AddAt(fFitCosmic,1);
1582 0 : }
1583 : // with naming convention and backward compatibility
1584 0 : fGainArray->AddAt(fGraphMIP,2);
1585 0 : fGainArray->AddAt(fGraphCosmic,3);
1586 : //
1587 : // 3.) Add HV and PT correction parameterization which was used
1588 : //
1589 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
1590 0 : if (param->GetGainSlopesHV()) fGainArray->AddLast(param->GetGainSlopesHV());
1591 0 : if (param->GetGainSlopesPT()) fGainArray->AddLast(param->GetGainSlopesPT());
1592 : //
1593 : // 4.) Add the RecoParam and ClusterParam - for compatibility checks -deffrent sets of paramters can invalidate calibration
1594 : //
1595 0 : AliTPCClusterParam *clParam = AliTPCcalibDB::Instance()->GetClusterParam();
1596 0 : TObjArray *recoParams = new TObjArray(4) ;
1597 0 : for (Int_t i=0;i<4;i++) recoParams->AddAt(AliTPCcalibDB::Instance()->GetRecoParam(i),i);
1598 0 : fGainArray->AddLast(clParam);
1599 0 : fGainArray->AddLast(recoParams);
1600 : //
1601 0 : cout << "fGraphCosmic: " << fGraphCosmic << " fGraphMIP " << fGraphMIP << endl;
1602 0 : return kTRUE;
1603 :
1604 0 : }
1605 :
1606 : //_____________________________________________________________________________
1607 : Bool_t AliTPCPreprocessorOffline::AnalyzeAttachment(Int_t startRunNumber, Int_t endRunNumber, Int_t minEntriesFit) {
1608 : //
1609 : // determine slope as a function of mean driftlength
1610 : //
1611 0 : if(!fGainMIP) return kFALSE;
1612 :
1613 0 : fGainMIP->GetHistGainTime()->GetAxis(5)->SetRangeUser(startRunNumber, endRunNumber);
1614 : //
1615 0 : fGainMIP->GetHistGainTime()->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
1616 0 : fGainMIP->GetHistGainTime()->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
1617 : //
1618 0 : fGainMIP->GetHistGainTime()->GetAxis(3)->SetRangeUser(125,250);// only full tracking region (driftlength)
1619 0 : fGainMIP->GetHistGainTime()->GetAxis(0)->SetRangeUser(1.5,3.5);// only full tracking region (driftlength)
1620 : //
1621 0 : TH3D * hist = fGainMIP->GetHistGainTime()->Projection(1, 0, 3);
1622 : //
1623 0 : Double_t *xvec = new Double_t[hist->GetNbinsX()];
1624 0 : Double_t *yvec = new Double_t[hist->GetNbinsX()];
1625 0 : Double_t *xerr = new Double_t[hist->GetNbinsX()];
1626 0 : Double_t *yerr = new Double_t[hist->GetNbinsX()];
1627 : Int_t counter = 0;
1628 : //
1629 0 : for(Int_t i=1; i < hist->GetNbinsX(); i++) {
1630 : Int_t nsum=0;
1631 : Int_t imin = i;
1632 : Int_t imax = i;
1633 0 : for (Int_t idelta=0; idelta<5; idelta++){
1634 : //
1635 0 : imin = TMath::Max(i-idelta,1);
1636 0 : imax = TMath::Min(i+idelta,hist->GetNbinsX());
1637 0 : nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
1638 : //if (nsum==0) break;
1639 0 : if (nsum>minEntriesFit) break;
1640 : }
1641 0 : if (nsum<minEntriesFit) continue;
1642 : //
1643 0 : fGainMIP->GetHistGainTime()->GetAxis(1)->SetRangeUser(hist->GetXaxis()->GetBinCenter(imin-1),hist->GetXaxis()->GetBinCenter(imax+1)); // define time range
1644 0 : TH2D * histZdep = fGainMIP->GetHistGainTime()->Projection(0,3);
1645 0 : TObjArray arr;
1646 0 : histZdep->FitSlicesY(0,0,-1,0,"QNR",&arr);
1647 0 : TH1D * driftDep = (TH1D*)arr.At(1);
1648 0 : delete histZdep;
1649 : //TGraphErrors * driftDep = AliTPCcalibBase::FitSlices(fGainMIP->GetHistGainTime(),0,3,100,1,0.,1);
1650 : /*if (driftDep->GetN() < 4) {
1651 : delete driftDep;
1652 : continue;
1653 : }*/
1654 : //
1655 : //TObjArray arr;
1656 : //
1657 0 : TF1 pol1("polynom1","pol1",125,240);
1658 : //driftDep->Fit(&pol1,"QNRROB=0.8");
1659 0 : driftDep->Fit(&pol1,"QNR");
1660 0 : xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin-1)+hist->GetXaxis()->GetBinCenter(imax+1));
1661 0 : yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
1662 0 : xerr[counter] = hist->GetXaxis()->GetBinCenter(imax+1)-hist->GetXaxis()->GetBinCenter(imin-1);
1663 0 : yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
1664 0 : counter++;
1665 : //
1666 : //delete driftDep;
1667 0 : }
1668 : //
1669 0 : fGraphAttachmentMIP = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
1670 0 : if (fGraphAttachmentMIP) fGraphAttachmentMIP->SetName("TGRAPHERRORS_MEAN_ATTACHMENT_BEAM_ALL");// set proper names according to naming convention
1671 0 : fGainArray->AddLast(fGraphAttachmentMIP);
1672 : //
1673 0 : delete [] xvec;
1674 0 : delete [] yvec;
1675 0 : delete [] xerr;
1676 0 : delete [] yerr;
1677 0 : delete hist;
1678 : //
1679 0 : if (counter < 1) return kFALSE;
1680 0 : return kTRUE;
1681 :
1682 0 : }
1683 :
1684 : //_____________________________________________________________________________
1685 : Bool_t AliTPCPreprocessorOffline::AnalyzePadRegionGain(){
1686 : //
1687 : // Analyze gain for different pad regions - produce the calibration graphs 0,1,2
1688 : // Uses the MIP pions integrated over all chambers, eta and multiplicity
1689 : // TODO: o What happens in case the MIPs are not equally distributed over
1690 : // eta and multiplicity? In CPass0 neither will be flat.
1691 : // Can we fill the histograms flat in eta (and multiplicity)?
1692 : // o In principle one could normalise over eta, by first projecting in slices
1693 : // of eta, normalize each slice and then sum
1694 : //
1695 0 : if (!fGainMult) return kFALSE;
1696 :
1697 0 : THnSparseF *histPadEqual=fGainMult->GetHistPadEqual();
1698 : // histPadEqual->GetAxis(4)->SetRangeUser(.35,.65);
1699 : // TH2D * histQmaxTmp = (TH2D*) histPadEqual->Projection(0,2);
1700 : // TH2D * histQtotTmp = (TH2D*) histPadEqual->Projection(1,2);
1701 : // histQmaxTmp->SetDirectory(0);
1702 : // histQtotTmp->SetDirectory(0);
1703 : // histPadEqual->GetAxis(4)->SetRangeUser(-.65,-.35);
1704 : // TH2D * histQmax = (TH2D*) histPadEqual->Projection(0,2);
1705 : // TH2D * histQtot = (TH2D*) histPadEqual->Projection(1,2);
1706 : // histQmax->Add(histQmaxTmp);
1707 : // histQtot->Add(histQtotTmp);
1708 : // delete histQmaxTmp;
1709 : // delete histQtotTmp;
1710 0 : TH2 * histQmax = AliTPCcalibBase::NormalizedProjection(histPadEqual,0,2,4);
1711 0 : TH2 * histQtot = AliTPCcalibBase::NormalizedProjection(histPadEqual,1,2,4);
1712 :
1713 : // === old part
1714 : // TObjArray arr;
1715 : // histQmax->FitSlicesY(0,0,-1,0,"QNR",&arr);
1716 : // Double_t xMax[3] = {0,1,2};
1717 : // Double_t yMax[3] = {((TH1D*)arr.At(1))->GetBinContent(1),
1718 : // ((TH1D*)arr.At(1))->GetBinContent(2),
1719 : // ((TH1D*)arr.At(1))->GetBinContent(3)};
1720 : // Double_t yMaxErr[3] = {((TH1D*)arr.At(1))->GetBinError(1),
1721 : // ((TH1D*)arr.At(1))->GetBinError(2),
1722 : // ((TH1D*)arr.At(1))->GetBinError(3)};
1723 : // //
1724 : // histQtot->FitSlicesY(0,0,-1,0,"QNR",&arr);
1725 : // Double_t xTot[3] = {0,1,2};
1726 : // Double_t yTot[3] = {((TH1D*)arr.At(1))->GetBinContent(1),
1727 : // ((TH1D*)arr.At(1))->GetBinContent(2),
1728 : // ((TH1D*)arr.At(1))->GetBinContent(3)};
1729 : // Double_t yTotErr[3] = {((TH1D*)arr.At(1))->GetBinError(1),
1730 : // ((TH1D*)arr.At(1))->GetBinError(2),
1731 : // ((TH1D*)arr.At(1))->GetBinError(3)};
1732 : // Double_t truncationFactor=AliTPCcalibGainMult::GetTruncatedMeanPosition(yTot[0],yTot[1],yTot[2],1000);
1733 : // for (Int_t i=0;i<3; i++){
1734 : // yMax[i]*=truncationFactor;
1735 : // yTot[i]*=truncationFactor;
1736 : // }
1737 : //
1738 : // TGraphErrors * fitPadRegionQmax = new TGraphErrors(3, xMax, yMax, 0, yMaxErr);
1739 : // TGraphErrors * fitPadRegionQtot = new TGraphErrors(3, xTot, yTot, 0, yTotErr);
1740 : // ^^^ End old part
1741 :
1742 : // histQmax->GetXaxis()->SetRange(1,3);
1743 : // histQtot->GetXaxis()->SetRange(1,3);
1744 0 : TGraphErrors *fitPadRegionQmax = AliTPCcalibBase::FitSlices(histQmax,200,1,.15,.85);
1745 0 : TGraphErrors *fitPadRegionQtot = AliTPCcalibBase::FitSlices(histQtot,200,1,.15,.85);
1746 0 : histQmax->GetXaxis()->SetRange(0,-1);
1747 0 : histQtot->GetXaxis()->SetRange(0,-1);
1748 0 : fitPadRegionQmax->RemovePoint(3);
1749 0 : fitPadRegionQtot->RemovePoint(3);
1750 :
1751 0 : Double_t *yMax =fitPadRegionQmax->GetY();
1752 0 : Double_t *yTot =fitPadRegionQtot->GetY();
1753 0 : Double_t *eyMax=fitPadRegionQmax->GetEY();
1754 0 : Double_t *eyTot=fitPadRegionQtot->GetEY();
1755 0 : Double_t truncationFactorMax=AliTPCcalibGainMult::GetTruncatedMeanPosition(yMax[0],yMax[1],yMax[2],1000);
1756 0 : Double_t truncationFactorTot=AliTPCcalibGainMult::GetTruncatedMeanPosition(yTot[0],yTot[1],yTot[2],1000);
1757 :
1758 0 : if (truncationFactorMax>0) {
1759 0 : for (Int_t i=0;i<3; i++){
1760 0 : yMax[i] /=truncationFactorMax;
1761 0 : eyMax[i]/=truncationFactorMax;
1762 : }
1763 0 : }
1764 : else {
1765 0 : AliError("truncationFactorMax<=0, could not normalize");
1766 : }
1767 :
1768 0 : if (truncationFactorTot>0) {
1769 0 : for (Int_t i=0;i<3; i++){
1770 0 : yTot[i] /=truncationFactorTot;
1771 0 : eyTot[i]/=truncationFactorTot;
1772 : }
1773 0 : }
1774 : else {
1775 0 : AliError("truncationFactorTot<=0, could not normalize");
1776 : }
1777 :
1778 : //
1779 0 : fitPadRegionQtot->SetName("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");// set proper names according to naming convention
1780 0 : fitPadRegionQmax->SetName("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");// set proper names according to naming convention
1781 : //
1782 0 : fGainArray->AddLast(fitPadRegionQtot);
1783 0 : fGainArray->AddLast(fitPadRegionQmax);
1784 :
1785 : // ---| QA histograms |---
1786 0 : if (fArrQAhist) {
1787 : // --- Qmax ---
1788 0 : histQmax->SetName("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL_QA");
1789 0 : histQmax->SetTitle("Pad region calibration Q_{max}; pad region;d#it{E}/d#it{x}_{Qmax} (arb. unit)");
1790 0 : fArrQAhist->Add(histQmax);
1791 :
1792 : // --- Qtot ---
1793 0 : histQtot->SetName("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL_QA");
1794 0 : histQtot->SetTitle("Pad region calibration Q_{tot}; pad region;d#it{E}/d#it{x}_{Qtot} (arb. unit)");
1795 0 : fArrQAhist->Add(histQtot);
1796 :
1797 0 : } else {
1798 0 : delete histQmax;
1799 0 : delete histQtot;
1800 : }
1801 :
1802 : return kTRUE;
1803 0 : }
1804 :
1805 : //_____________________________________________________________________________
1806 : Bool_t AliTPCPreprocessorOffline::AnalyzeGainDipAngle(Int_t padRegion) {
1807 : //
1808 : // Analyze gain as a function of multiplicity and produce calibration graphs
1809 : // padRegion -- 0: short, 1: medium, 2: long, 3: absolute calibration of full track
1810 : // Uses the MIP pions integrated over all chambers and multiplicity
1811 : // Noramlisation is done to the mean of the slices fit (i.e. tan(lambda)=0.5)
1812 : //
1813 : // TODO: What happens in case the MIPs are not equally distributed over
1814 : // multiplicity? In CPass0 it will not be flat.
1815 : // Can we fill the histograms flat in multiplicity?
1816 : //
1817 0 : Int_t kMarkers[10]={25,24,20,21,22};
1818 0 : Int_t kColors[10]={1,2,4,3,6};
1819 0 : if (!fGainMult) return kFALSE;
1820 0 : if (!(fGainMult->GetHistTopology())) return kFALSE;
1821 : const Double_t kMinStat=100;
1822 : //
1823 : // "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength"
1824 0 : TObjArray arrMax;
1825 0 : TObjArray arrTot;
1826 : //
1827 : TH2D * histQmax = 0x0;
1828 : TH2D * histQtot = 0x0;
1829 0 : fGainMult->GetHistPadEqual()->GetAxis(4)->SetRangeUser(-0.85,0.85);
1830 0 : fGainMult->GetHistTopology()->GetAxis(2)->SetRangeUser(-0.85,0.85);
1831 0 : if (padRegion < 3) {
1832 0 : fGainMult->GetHistPadEqual()->GetAxis(2)->SetRangeUser(padRegion,padRegion); // short,medium,long
1833 0 : histQmax = (TH2D*) fGainMult->GetHistPadEqual()->Projection(0,4);
1834 0 : histQtot = (TH2D*) fGainMult->GetHistPadEqual()->Projection(1,4);
1835 0 : } else {
1836 0 : fGainMult->GetHistTopology()->GetAxis(1)->SetRangeUser(1,1); //Qmax
1837 0 : histQmax = (TH2D*) fGainMult->GetHistTopology()->Projection(0,2);
1838 0 : histQmax->SetName("fGainMult_GetHistPadEqual_11");
1839 0 : fGainMult->GetHistTopology()->GetAxis(1)->SetRangeUser(0,0); //Qtot
1840 0 : histQtot = (TH2D*) fGainMult->GetHistTopology()->Projection(0,2);
1841 0 : histQtot->SetName("fGainMult_GetHistPadEqual_00");
1842 : }
1843 : //
1844 :
1845 0 : if (histQmax->GetEntries()<=kMinStat || histQtot->GetEntries()<=kMinStat) {
1846 0 : AliError(Form("hisQtot.GetEntries()=%f",histQtot->GetEntries()));
1847 0 : AliError(Form("hisQmax.GetEntries()=%f",histQmax->GetEntries()));
1848 0 : return kFALSE;
1849 : }
1850 :
1851 : // TGraphErrors * graphMax = TStatToolkit::MakeStat1D( histQmax,0,0.8,4,kMarkers[padRegion],kColors[padRegion]);
1852 : // TGraphErrors * graphTot = TStatToolkit::MakeStat1D( histQtot,0,0.8,4,kMarkers[padRegion],kColors[padRegion]);
1853 0 : TGraphErrors * graphMax = TStatToolkit::MakeStat1D( histQmax,0,0.9,6,kMarkers[padRegion],kColors[padRegion]);
1854 0 : TGraphErrors * graphTot = TStatToolkit::MakeStat1D( histQtot,0,0.9,6,kMarkers[padRegion],kColors[padRegion]);
1855 :
1856 : //
1857 0 : const char* names[4]={"SHORT","MEDIUM","LONG","ABSOLUTE"};
1858 : //
1859 0 : const Double_t meanMax = TMath::Mean(graphMax->GetN(), graphMax->GetY());
1860 0 : const Double_t meanTot = TMath::Mean(graphTot->GetN(), graphTot->GetY());
1861 :
1862 : // ---| QA histograms |---
1863 0 : if (fArrQAhist) {
1864 : // --- Qmax ---
1865 0 : histQmax->SetName(Form("TGRAPHERRORS_QMAX_DIPANGLE_%s_BEAM_ALL_QA",names[padRegion]));
1866 0 : histQmax->SetTitle(Form("tan(#lambda) calibration Q_{max} %s; tan(#lambda);d#it{E}/d#it{x}_{Qmax} (arb. unit)",names[padRegion]));
1867 0 : fArrQAhist->Add(histQmax);
1868 :
1869 : // --- Qtot ---
1870 0 : histQtot->SetName(Form("TGRAPHERRORS_QTOT_DIPANGLE_%s_BEAM_ALL_QA",names[padRegion]));
1871 0 : histQtot->SetTitle(Form("tan(#lambda) calibration Q_{tot} %s; tan(#lambda);d#it{E}/d#it{x}_{Qtot} (arb. unit)",names[padRegion]));
1872 0 : fArrQAhist->Add(histQtot);
1873 :
1874 : // ---| scale to mean multiplicity |---
1875 0 : if(fNormaliseQA && meanMax>0) {
1876 0 : TAxis *a=histQmax->GetYaxis();
1877 0 : a->SetLimits(a->GetXmin()/meanMax, a->GetXmax()/meanMax);
1878 0 : }
1879 0 : if(fNormaliseQA && meanTot>0) {
1880 0 : TAxis *a=histQtot->GetYaxis();
1881 0 : a->SetLimits(a->GetXmin()/meanTot, a->GetXmax()/meanTot);
1882 0 : }
1883 :
1884 : } else {
1885 0 : delete histQmax;
1886 0 : delete histQtot;
1887 : }
1888 :
1889 : //
1890 0 : if (meanMax<=0 || meanTot<=0){
1891 0 : AliError(Form("meanMax=%f",meanMax));
1892 0 : AliError(Form("meanTot=%f",meanTot));
1893 0 : return kFALSE;
1894 : }
1895 : //
1896 0 : for (Int_t ipoint=0; ipoint<graphMax->GetN(); ipoint++) {
1897 0 : graphMax->GetY()[ipoint]/=meanMax;
1898 0 : graphMax->GetEY()[ipoint]/=meanMax;
1899 : }
1900 0 : for (Int_t ipoint=0; ipoint<graphTot->GetN(); ipoint++) {
1901 0 : graphTot->GetY()[ipoint]/=meanTot;
1902 0 : graphTot->GetEY()[ipoint]/=meanTot;
1903 : }
1904 : //
1905 0 : graphMax->SetNameTitle(Form("TGRAPHERRORS_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]),
1906 0 : Form("TGRAPHERRORS_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]));
1907 0 : graphTot->SetNameTitle(Form("TGRAPHERRORS_QTOT_DIPANGLE_%s_BEAM_ALL",names[padRegion]),
1908 0 : Form("TGRAPHERRORS_QTOT_DIPANGLE_%s_BEAM_ALL",names[padRegion]));
1909 : //
1910 0 : fGainArray->AddLast(graphMax);
1911 0 : fGainArray->AddLast(graphTot);
1912 : //
1913 : // Normalization to 1 (mean of the graph.fY --> 1)
1914 : //
1915 0 : TF1 * funMax= new TF1("","1++abs(x)++abs(x*x)");
1916 0 : TF1 * funTot= new TF1("","1++abs(x)++abs(x*x)");
1917 0 : graphMax->Fit(funMax,"w","rob=0.9",-0.8,0.8);
1918 0 : graphTot->Fit(funTot,"w","rob=0.9",-0.8,0.8);
1919 0 : funMax->SetNameTitle(Form("TF1_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]),
1920 0 : Form("TF1_QMAX_DIPANGLE_%s_BEAM_ALL",names[padRegion]));
1921 0 : funTot->SetNameTitle(Form("TF1_QTOT_DIPANGLE_%s_BEAM_ALL",names[padRegion]),
1922 0 : Form("TF1_QTOT_DIPANGLE_%s_BEAM_ALL",names[padRegion]));
1923 :
1924 : //
1925 0 : fGainArray->AddLast(funMax);
1926 0 : fGainArray->AddLast(funTot);
1927 : //
1928 : return kTRUE;
1929 0 : }
1930 :
1931 : //_____________________________________________________________________________
1932 : Bool_t AliTPCPreprocessorOffline::AnalyzeGainMultiplicity() {
1933 : //
1934 : // Analyze gain as a function of multiplicity and produce calibration graphs
1935 : // Uses the MIP pions integrated over all chambers and eta
1936 : // Noramlisation is done to the mean of the slices fit (i.e. sector average)
1937 : //
1938 : // TODO: o What happens in case the MIPs are not flat in eta?
1939 : // In CPass0 it will not be flat.
1940 : // Can we fill the histograms flat in multiplicity?
1941 : // o Is it better to use the median in case one chamber fit fails?
1942 : // o The gain mult does not have tan(lambda) as dim
1943 : //
1944 0 : if (!fGainMult) return kFALSE;
1945 0 : fGainMult->GetHistGainMult()->GetAxis(3)->SetRangeUser(3,3);
1946 0 : TH2D * histMultMax = fGainMult->GetHistGainMult()->Projection(0,4);
1947 0 : TH2D * histMultTot = fGainMult->GetHistGainMult()->Projection(1,4);
1948 0 : histMultMax->RebinX(4);
1949 0 : histMultTot->RebinX(4);
1950 : //
1951 0 : TObjArray arrMax;
1952 0 : TObjArray arrTot;
1953 0 : TF1 fitGaus("fitGaus","gaus(0)",histMultMax->GetYaxis()->GetXmin(),histMultMax->GetYaxis()->GetXmax());
1954 0 : fitGaus.SetParameters(histMultMax->GetEntries()/10., histMultMax->GetMean(2), TMath::Sqrt(TMath::Abs(histMultMax->GetMean(2))));
1955 0 : histMultMax->FitSlicesY(&fitGaus,0,-1,1,"QNRB",&arrMax);
1956 0 : fitGaus.SetParameters(histMultTot->GetEntries()/10., histMultTot->GetMean(2), TMath::Sqrt(TMath::Abs(histMultTot->GetMean(2))));
1957 0 : histMultTot->FitSlicesY(&fitGaus,0,-1,1,"QNRB",&arrTot);
1958 : //
1959 0 : TH1D * meanMax = (TH1D*)arrMax.At(1);
1960 0 : TH1D * meanTot = (TH1D*)arrTot.At(1);
1961 0 : Float_t meanMult = histMultMax->GetMean();
1962 0 : const Double_t qMaxCont=meanMax->GetBinContent(meanMax->FindBin(meanMult));
1963 0 : const Double_t qTotCont=meanTot->GetBinContent(meanTot->FindBin(meanMult));
1964 :
1965 : // ---| QA histograms |---
1966 0 : if (fArrQAhist) {
1967 : // --- Qmax ---
1968 0 : histMultMax->SetName("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL_QA");
1969 0 : histMultMax->SetTitle("Multiplicity correction Q_{max};#ESD tracks;d#it{E}/d#it{x}_{Qmax} (arb. unit)");
1970 0 : fArrQAhist->Add(histMultMax);
1971 :
1972 : // --- Qtot ---
1973 0 : histMultTot->SetName("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL_QA");
1974 0 : histMultTot->SetTitle("Multiplicity correction Q_{tot};#ESD tracks;d#it{E}/d#it{x}_{Qtot} (arb. unit)");
1975 0 : fArrQAhist->Add(histMultTot);
1976 :
1977 : // ---| scale to mean multiplicity |---
1978 0 : if(fNormaliseQA && qMaxCont>0) {
1979 0 : TAxis *a=histMultMax->GetYaxis();
1980 0 : a->SetLimits(a->GetXmin()/qMaxCont, a->GetXmax()/qMaxCont);
1981 0 : }
1982 0 : if(fNormaliseQA && qTotCont>0) {
1983 0 : TAxis *a=histMultTot->GetYaxis();
1984 0 : a->SetLimits(a->GetXmin()/qTotCont, a->GetXmax()/qTotCont);
1985 0 : }
1986 :
1987 : } else {
1988 0 : delete histMultMax;
1989 0 : delete histMultTot;
1990 : }
1991 :
1992 : // ---| scale to mean multiplicity |---
1993 0 : if(qMaxCont) {
1994 0 : meanMax->Scale(1./qMaxCont);
1995 : }
1996 : else {
1997 0 : return kFALSE;
1998 : }
1999 0 : if(qTotCont) {
2000 0 : meanTot->Scale(1./qTotCont);
2001 : }
2002 : else {
2003 0 : return kFALSE;
2004 : }
2005 0 : Float_t xMultMax[50];
2006 0 : Float_t yMultMax[50];
2007 0 : Float_t yMultErrMax[50];
2008 0 : Float_t xMultTot[50];
2009 0 : Float_t yMultTot[50];
2010 0 : Float_t yMultErrTot[50];
2011 : //
2012 : Int_t nCountMax = 0;
2013 0 : for(Int_t iBin = 1; iBin < meanMax->GetXaxis()->GetNbins(); iBin++) {
2014 0 : Float_t yValMax = meanMax->GetBinContent(iBin);
2015 0 : if (yValMax < 0.7) continue;
2016 0 : if (yValMax > 1.3) continue;
2017 0 : if (meanMax->GetBinError(iBin)/yValMax > 0.01) continue;
2018 0 : xMultMax[nCountMax] = meanMax->GetXaxis()->GetBinCenter(iBin);
2019 0 : yMultMax[nCountMax] = yValMax;
2020 0 : yMultErrMax[nCountMax] = meanMax->GetBinError(iBin);
2021 0 : nCountMax++;
2022 0 : }
2023 : //
2024 0 : if (nCountMax < 10) return kFALSE;
2025 0 : TGraphErrors * fitMultMax = new TGraphErrors(nCountMax, xMultMax, yMultMax, 0, yMultErrMax);
2026 0 : fitMultMax->SetName("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2027 : //
2028 : Int_t nCountTot = 0;
2029 0 : for(Int_t iBin = 1; iBin < meanTot->GetXaxis()->GetNbins(); iBin++) {
2030 0 : Float_t yValTot = meanTot->GetBinContent(iBin);
2031 0 : if (yValTot < 0.7) continue;
2032 0 : if (yValTot > 1.3) continue;
2033 0 : if (meanTot->GetBinError(iBin)/yValTot > 0.1) continue;
2034 0 : xMultTot[nCountTot] = meanTot->GetXaxis()->GetBinCenter(iBin);
2035 0 : yMultTot[nCountTot] = yValTot;
2036 0 : yMultErrTot[nCountTot] = meanTot->GetBinError(iBin);
2037 0 : nCountTot++;
2038 0 : }
2039 : //
2040 0 : if (nCountTot < 10) return kFALSE;
2041 0 : TGraphErrors * fitMultTot = new TGraphErrors(nCountTot, xMultTot, yMultTot, 0, yMultErrTot);
2042 0 : fitMultTot->SetName("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2043 : //
2044 0 : fGainArray->AddLast(fitMultMax);
2045 0 : fGainArray->AddLast(fitMultTot);
2046 : //
2047 : return kTRUE;
2048 :
2049 0 : }
2050 :
2051 : //_____________________________________________________________________________
2052 : Bool_t AliTPCPreprocessorOffline::AnalyzeGainChamberByChamber(){
2053 : //
2054 : // get chamber by chamber gain
2055 : //
2056 0 : if (!fGainMult) return kFALSE;
2057 0 : TGraphErrors *grShort = fGainMult->GetGainPerChamberRobust(0, kFALSE, fArrQAhist, fNormaliseQA);
2058 0 : TGraphErrors *grMedium = fGainMult->GetGainPerChamberRobust(1, kFALSE, fArrQAhist, fNormaliseQA);
2059 0 : TGraphErrors *grLong = fGainMult->GetGainPerChamberRobust(2, kFALSE, fArrQAhist, fNormaliseQA);
2060 0 : if (grShort==0x0 || grMedium==0x0 || grLong==0x0) {
2061 0 : delete grShort;
2062 0 : delete grMedium;
2063 0 : delete grLong;
2064 0 : return kFALSE;
2065 : }
2066 :
2067 0 : fGainArray->AddLast(grShort);
2068 0 : fGainArray->AddLast(grMedium);
2069 0 : fGainArray->AddLast(grLong);
2070 :
2071 0 : return kTRUE;
2072 0 : }
2073 :
2074 : //_____________________________________________________________________________
2075 : void AliTPCPreprocessorOffline::UpdateOCDBGain(Int_t startRunNumber, Int_t endRunNumber, AliCDBStorage *fullStorage, AliCDBStorage* residualStorage/*=0x0*/)
2076 : {
2077 : //
2078 : // Update OCDB entry
2079 : //
2080 0 : AliCDBMetaData *metaData= new AliCDBMetaData();
2081 0 : metaData->SetObjectClassName("TObjArray");
2082 0 : metaData->SetResponsible("Jens Wiechula");
2083 0 : metaData->SetBeamPeriod(1);
2084 0 : metaData->SetAliRootVersion("05-24-00"); //root version
2085 0 : metaData->SetComment("Calibration of the time dependence of the gain due to pressure and temperature changes.");
2086 0 : AliCDBId id1("TPC/Calib/TimeGain", startRunNumber, endRunNumber);
2087 0 : if (fGainCalibrationType==kFullGainCalib) {
2088 0 : AliInfoF("Writing gain calibration object to the full storage: %s", fullStorage->GetURI().Data());
2089 0 : fullStorage->Put(fGainArray, id1, metaData);
2090 : }
2091 0 : else if (fGainCalibrationType==kResidualGainQA) {
2092 0 : AliInfoF("Writing gain calibration object to the residual storage: %s", residualStorage->GetURI().Data());
2093 0 : residualStorage->Put(fGainArray, id1, metaData);
2094 : }
2095 0 : else if (fGainCalibrationType==kCombinedGainCalib) {
2096 0 : if (residualStorage){
2097 0 : AliInfoF("Writing residual gain calibration object to the residual storage: %s", residualStorage->GetURI().Data());
2098 0 : residualStorage->Put(fGainArray, id1, metaData);
2099 : }
2100 : else {
2101 0 : AliError("No residual storage set, but calibration type is combined + residual QA");
2102 : }
2103 :
2104 : // write the combined calibration after the residual calibration to make sure this is the
2105 : // latest object in case fullStorage and residualStorage are identical
2106 0 : AliInfoF("Writing combined gain calibration object to the full storage: %s", fullStorage->GetURI().Data());
2107 0 : fullStorage->Put(fGainArrayCombined, id1, metaData);
2108 : }
2109 : else {
2110 0 : AliFatalF("Unsupported gain calibration type: %d", Int_t(fGainCalibrationType));
2111 : }
2112 0 : }
2113 :
2114 : //_____________________________________________________________________________
2115 : void AliTPCPreprocessorOffline::MakeQAPlot(Float_t FPtoMIPratio) {
2116 : //
2117 : // Make QA plot to visualize results
2118 : //
2119 : //
2120 : //
2121 0 : if (fGraphCosmic) {
2122 0 : TCanvas * canvasCosmic = new TCanvas("gain Cosmic", "time dependent gain QA histogram cosmic");
2123 0 : canvasCosmic->cd();
2124 0 : TH2D * gainHistoCosmic = fGainCosmic->GetHistGainTime()->Projection(0,1);
2125 0 : gainHistoCosmic->SetDirectory(0);
2126 0 : gainHistoCosmic->SetName("GainHistoCosmic");
2127 0 : gainHistoCosmic->GetXaxis()->SetTimeDisplay(kTRUE);
2128 0 : gainHistoCosmic->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
2129 0 : gainHistoCosmic->Draw("colz");
2130 0 : fGraphCosmic->SetMarkerStyle(25);
2131 0 : fGraphCosmic->Draw("lp");
2132 0 : fGraphCosmic->SetMarkerStyle(25);
2133 0 : TGraph * grfFitCosmic = fFitCosmic->MakeGraph(fGraphCosmic->GetX()[0],fGraphCosmic->GetX()[fGraphCosmic->GetN()-1],50000,0);
2134 0 : if (grfFitCosmic) {
2135 0 : for(Int_t i=0; i < grfFitCosmic->GetN(); i++) {
2136 0 : grfFitCosmic->GetY()[i] *= FPtoMIPratio;
2137 : }
2138 0 : for(Int_t i=0; i < fGraphCosmic->GetN(); i++) {
2139 0 : fGraphCosmic->GetY()[i] *= FPtoMIPratio;
2140 : }
2141 0 : }
2142 0 : fGraphCosmic->Draw("lp");
2143 0 : if (grfFitCosmic) {
2144 0 : grfFitCosmic->SetLineColor(2);
2145 0 : grfFitCosmic->Draw("lu");
2146 0 : }
2147 0 : fGainArray->AddLast(gainHistoCosmic);
2148 : //fGainArray->AddLast(canvasCosmic->Clone());
2149 0 : delete canvasCosmic;
2150 0 : }
2151 0 : if (fFitMIP) {
2152 0 : TCanvas * canvasMIP = new TCanvas("gain MIP", "time dependent gain QA histogram MIP");
2153 0 : canvasMIP->cd();
2154 0 : TH2D * gainHistoMIP = fGainMIP->GetHistGainTime()->Projection(0,1);
2155 0 : gainHistoMIP->SetName("GainHistoCosmic");
2156 0 : gainHistoMIP->SetDirectory(0);
2157 0 : gainHistoMIP->GetXaxis()->SetTimeDisplay(kTRUE);
2158 0 : gainHistoMIP->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
2159 0 : gainHistoMIP->Draw("colz");
2160 0 : fGraphMIP->SetMarkerStyle(25);
2161 0 : fGraphMIP->Draw("lp");
2162 0 : TGraph * grfFitMIP = fFitMIP->MakeGraph(fGraphMIP->GetX()[0],fGraphMIP->GetX()[fGraphMIP->GetN()-1],50000,0);
2163 0 : grfFitMIP->SetLineColor(2);
2164 0 : grfFitMIP->Draw("lu");
2165 0 : fGainArray->AddLast(gainHistoMIP);
2166 : //fGainArray->AddLast(canvasMIP->Clone());
2167 0 : delete canvasMIP;
2168 0 : delete grfFitMIP;
2169 0 : }
2170 0 : }
2171 :
2172 : //_____________________________________________________________________________
2173 : void AliTPCPreprocessorOffline::MakeFitTime(){
2174 : //
2175 : // make aligment fit - store results in the file
2176 : //
2177 : const Int_t kMinEntries=1000;
2178 0 : MakeChainTime();
2179 0 : MakePrimitivesTime();
2180 0 : if (!fAlignTree) return;
2181 0 : if (fAlignTree->GetEntries()<kMinEntries) return;
2182 0 : fAlignTree->SetAlias("ptype","type");
2183 0 : fAlignTree->SetAlias("hasITS","(1+0)");
2184 0 : fAlignTree->SetAlias("dITS","1-2*(refX<40)");
2185 0 : fAlignTree->SetAlias("isITS","refX>10");
2186 0 : fAlignTree->SetAlias("isVertex","refX<10");
2187 : //
2188 : Int_t npointsMax=30000000;
2189 0 : Double_t chi2=0;
2190 0 : Int_t npoints=0;
2191 0 : TVectorD param;
2192 0 : TMatrixD covar;
2193 :
2194 0 : TString fstringFast="";
2195 0 : fstringFast+="FExBTwistX++";
2196 0 : fstringFast+="FExBTwistY++";
2197 0 : fstringFast+="FAlignRot0D++";
2198 0 : fstringFast+="FAlignTrans0D++";
2199 0 : fstringFast+="FAlignTrans1D++";
2200 : //
2201 0 : fstringFast+="hasITS*FAlignTrans0++";
2202 0 : fstringFast+="hasITS*FAlignTrans1++";
2203 0 : fstringFast+="hasITS*FAlignRot0++";
2204 0 : fstringFast+="hasITS*FAlignRot1++";
2205 0 : fstringFast+="hasITS*FAlignRot2++";
2206 : //
2207 0 : fstringFast+="dITS*FAlignTrans0++";
2208 0 : fstringFast+="dITS*FAlignTrans1++";
2209 0 : fstringFast+="dITS*FAlignRot0++";
2210 0 : fstringFast+="dITS*FAlignRot1++";
2211 0 : fstringFast+="dITS*FAlignRot2++";
2212 :
2213 0 : TCut cutFit="entries>10&&abs(mean)>0.00001&&rms>0";
2214 0 : fAlignTree->SetAlias("err","rms");
2215 :
2216 0 : TString *strDeltaITS = TStatToolkit::FitPlaneConstrain(fAlignTree,"mean:err", fstringFast.Data(),cutFit, chi2,npoints,param,covar,-1,0, npointsMax, 1);
2217 0 : TObjArray* tokArr = strDeltaITS->Tokenize("++");
2218 0 : static bool verboseOutput = !(getenv("HLT_ONLINE_MODE") && strcmp(getenv("HLT_ONLINE_MODE"), "on") == 0);
2219 0 : if (verboseOutput) tokArr->Print();
2220 0 : delete tokArr;
2221 0 : fAlignTree->SetAlias("fitYFast",strDeltaITS->Data());
2222 0 : delete strDeltaITS;
2223 :
2224 : //
2225 0 : TVectorD paramC= param;
2226 0 : TMatrixD covarC= covar;
2227 0 : TStatToolkit::Constrain1D(fstringFast,"Trans0D",paramC,covarC,0, 0.1);
2228 0 : TStatToolkit::Constrain1D(fstringFast,"Trans1D",paramC,covarC,0, 0.1);
2229 0 : TStatToolkit::Constrain1D(fstringFast,"TwistX",paramC,covarC,0, 0.1);
2230 0 : TStatToolkit::Constrain1D(fstringFast,"TwistY",paramC,covarC,0, 0.1);
2231 0 : TString strFitConst=TStatToolkit::MakeFitString(fstringFast, paramC,covar);
2232 0 : fAlignTree->SetAlias("fitYFastC",strFitConst.Data());
2233 0 : CreateAlignTime(fstringFast,paramC);
2234 :
2235 0 : }
2236 :
2237 : //_____________________________________________________________________________
2238 : void AliTPCPreprocessorOffline::MakeChainTime(){
2239 : //
2240 : //
2241 : //
2242 0 : TFile f("CalibObjects.root");
2243 :
2244 : // const char *cdtype[7]={"ITS","TRD","Vertex","TOF","TPC","TPC0","TPC1"};
2245 : //const char *cptype[5]={"dy","dz","dsnp","dtheta","d1pt"};
2246 0 : const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
2247 : Int_t run=0;
2248 : AliTPCcalibTime *calibTime = 0;
2249 0 : TObject* obj = dynamic_cast<TObject*>(f.Get("TPCCalib"));
2250 0 : TObjArray * array = dynamic_cast<TObjArray*>(obj);
2251 0 : TDirectory * dir = dynamic_cast<TDirectory*>(obj);
2252 0 : if (dir) {
2253 0 : calibTime = dynamic_cast<AliTPCcalibTime*>(dir->Get("calibTime"));
2254 0 : }
2255 0 : else if (array){
2256 0 : calibTime = (AliTPCcalibTime *)array->FindObject("calibTime");
2257 0 : } else {
2258 0 : calibTime = (AliTPCcalibTime*)f.Get("calibTime");
2259 : }
2260 0 : if (!calibTime) return;
2261 0 : AliTPCCorrectionFit::CreateAlignMaps(AliTracker::GetBz(), run);
2262 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("meanITSVertex.root");
2263 : //
2264 : Int_t ihis=0;
2265 0 : THnSparse *his = calibTime->GetResHistoTPCITS(ihis);
2266 0 : if (his){
2267 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2268 0 : his->GetAxis(2)->SetRange(0,1000000);
2269 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2270 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,3);
2271 : }
2272 : ihis=1;
2273 0 : his = calibTime->GetResHistoTPCITS(ihis);
2274 0 : if (his){
2275 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2276 0 : his->GetAxis(2)->SetRange(0,1000000);
2277 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2278 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,3);
2279 : }
2280 : ihis=2;
2281 0 : his = calibTime->GetResHistoTPCITS(ihis);
2282 0 : if (his){
2283 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2284 0 : his->GetAxis(2)->SetRange(0,1000000);
2285 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2286 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("ITS%s",hname[ihis]),run,85.,ihis,3);
2287 : }
2288 : ihis=0;
2289 0 : his = calibTime->GetResHistoTPCvertex(ihis);
2290 0 : if (his){
2291 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2292 0 : his->GetAxis(2)->SetRange(0,1000000);
2293 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2294 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,3);
2295 : }
2296 : ihis=2;
2297 0 : his = calibTime->GetResHistoTPCvertex(ihis);
2298 0 : if (his){
2299 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2300 0 : his->GetAxis(2)->SetRange(0,1000000);
2301 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2302 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,3);
2303 :
2304 : }
2305 : ihis=1;
2306 0 : his = calibTime->GetResHistoTPCvertex(ihis);
2307 0 : if (his){
2308 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2309 0 : his->GetAxis(2)->SetRange(0,1000000);
2310 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2311 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("Vertex%s",hname[ihis]),run,0.,ihis,3);
2312 :
2313 : }
2314 : ihis=0;
2315 0 : his = calibTime->GetResHistoTPCTOF(ihis);
2316 0 : if (his){
2317 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2318 0 : his->GetAxis(2)->SetRange(0,1000000);
2319 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2320 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TOF%s",hname[ihis]),run,0.,ihis,3);
2321 :
2322 : }
2323 : ihis=0;
2324 0 : his = calibTime->GetResHistoTPCTRD(ihis);
2325 0 : if (his){
2326 0 : his->GetAxis(1)->SetRangeUser(-1.1,1.1);
2327 0 : his->GetAxis(2)->SetRange(0,1000000);
2328 0 : his->GetAxis(3)->SetRangeUser(-0.35,0.35);
2329 0 : AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TRD%s",hname[ihis]),run,0.,ihis,3);
2330 :
2331 : }
2332 0 : if (dir || (!dir && !array)) { // the object is taken from a directory or a file
2333 0 : delete calibTime;
2334 : }
2335 0 : delete pcstream;
2336 0 : }
2337 :
2338 : //_____________________________________________________________________________
2339 : Double_t AliTPCPreprocessorOffline::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){
2340 : //
2341 : //
2342 : //
2343 0 : Double_t sector = 9*phi/TMath::Pi();
2344 0 : if (sector<0) sector+=18;
2345 0 : Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
2346 0 : Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
2347 0 : if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
2348 0 : if (ptype==2) return (y245-y85)/(245.-85.);
2349 0 : return 0;
2350 0 : }
2351 :
2352 : //_____________________________________________________________________________
2353 : Double_t AliTPCPreprocessorOffline::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
2354 : //
2355 : // Fit the distortion along the line with the parabolic model
2356 : // Parameters:
2357 : // phi0 - phi at the entrance of the TPC
2358 : // snp - local inclination angle at the entrance of the TPC
2359 : // refX - ref X where the distortion is evanluated
2360 : // theta
2361 : //
2362 0 : static TLinearFitter fitter(3,"pol2");
2363 0 : fitter.ClearPoints();
2364 0 : if (nsteps<3) nsteps=3;
2365 0 : Double_t deltaX=(245-85)/(nsteps);
2366 0 : for (Int_t istep=0; istep<(nsteps+1); istep++){
2367 : //
2368 0 : Double_t localX =85.+deltaX*istep;
2369 0 : Double_t localPhi=phi0+deltaX*snp*istep;
2370 0 : Double_t sector = 9*localPhi/TMath::Pi();
2371 0 : if (sector<0) sector+=18;
2372 0 : Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
2373 0 : Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
2374 0 : Double_t x[1]={localX-dlocalX};
2375 0 : fitter.AddPoint(x,y);
2376 0 : }
2377 0 : fitter.Eval();
2378 : Double_t par[3];
2379 0 : par[0]=fitter.GetParameter(0);
2380 0 : par[1]=fitter.GetParameter(1);
2381 0 : par[2]=fitter.GetParameter(2);
2382 :
2383 0 : if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
2384 0 : if (ptype==2) return par[1]+2*par[2]*refX;
2385 0 : if (ptype==4) return par[2];
2386 0 : return 0;
2387 0 : }
2388 :
2389 : //_____________________________________________________________________________
2390 : void AliTPCPreprocessorOffline::MakePrimitivesTime(){
2391 : //
2392 : // Create primitive transformation to fit
2393 : //
2394 0 : fAlignTree=new TChain("fit","fit");
2395 0 : fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdy");
2396 0 : fAlignTree->AddFile("meanITSVertex.root",10000000,"ITSdsnp");
2397 0 : fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdy");
2398 0 : fAlignTree->AddFile("meanITSVertex.root",10000000,"Vertexdsnp");
2399 : //
2400 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
2401 0 : Double_t bzField=AliTrackerBase::GetBz();
2402 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
2403 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
2404 0 : Double_t wtP = -10.0 * (bzField) * vdrift / ezField ;
2405 0 : AliTPCExBTwist *fitExBTwistX= new AliTPCExBTwist;
2406 0 : AliTPCExBTwist *fitExBTwistY= new AliTPCExBTwist;
2407 0 : AliTPCCalibGlobalMisalignment *trans0 =new AliTPCCalibGlobalMisalignment;
2408 0 : AliTPCCalibGlobalMisalignment *trans1 =new AliTPCCalibGlobalMisalignment;
2409 0 : AliTPCCalibGlobalMisalignment *trans0D =new AliTPCCalibGlobalMisalignment;
2410 0 : AliTPCCalibGlobalMisalignment *trans1D =new AliTPCCalibGlobalMisalignment;
2411 0 : AliTPCCalibGlobalMisalignment *rot0 =new AliTPCCalibGlobalMisalignment;
2412 0 : AliTPCCalibGlobalMisalignment *rot1 =new AliTPCCalibGlobalMisalignment;
2413 0 : AliTPCCalibGlobalMisalignment *rot2 =new AliTPCCalibGlobalMisalignment;
2414 0 : AliTPCCalibGlobalMisalignment *rot3 =new AliTPCCalibGlobalMisalignment;
2415 : //
2416 : //
2417 0 : fitExBTwistX->SetXTwist(0.001);
2418 0 : fitExBTwistX->SetOmegaTauT1T2(wtP,1,1);
2419 : //
2420 0 : fitExBTwistY->SetYTwist(0.001);
2421 0 : fitExBTwistY->SetOmegaTauT1T2(wtP,1,1);
2422 : //
2423 0 : TGeoHMatrix *matrixRot = new TGeoHMatrix;
2424 0 : TGeoHMatrix *matrixX = new TGeoHMatrix;
2425 0 : TGeoHMatrix *matrixY = new TGeoHMatrix;
2426 0 : matrixX->SetDx(0.1);
2427 0 : matrixY->SetDy(0.1);
2428 0 : Double_t rotAngles0[9]={0};
2429 0 : Double_t rotAngles1[9]={0};
2430 0 : Double_t rotAngles2[9]={0};
2431 : //
2432 0 : Double_t rotAngles3[9]={0};
2433 :
2434 0 : rotAngles0[0]=1; rotAngles0[4]=1; rotAngles0[8]=1;
2435 0 : rotAngles1[0]=1; rotAngles1[4]=1; rotAngles1[8]=1;
2436 0 : rotAngles2[0]=1; rotAngles2[4]=1; rotAngles2[8]=1;
2437 0 : rotAngles3[0]=1; rotAngles3[4]=1; rotAngles3[8]=1;
2438 :
2439 0 : rotAngles0[1]=-0.001;rotAngles0[3]=0.001;
2440 0 : rotAngles1[5]=-0.001;rotAngles1[7]=0.001;
2441 0 : rotAngles2[2]=0.001;rotAngles2[6]=-0.001;
2442 0 : rotAngles3[1]=0.001;rotAngles3[3]=-0.001;
2443 0 : matrixRot->SetRotation(rotAngles0);
2444 0 : rot0->SetAlignGlobal(matrixRot);
2445 0 : matrixRot->SetRotation(rotAngles1);
2446 0 : rot1->SetAlignGlobal(matrixRot);
2447 0 : matrixRot->SetRotation(rotAngles2);
2448 0 : rot2->SetAlignGlobal(matrixRot);
2449 0 : matrixRot->SetRotation(rotAngles3);
2450 0 : rot3->SetAlignGlobalDelta(matrixRot);
2451 : //
2452 0 : trans0->SetAlignGlobal(matrixX);
2453 0 : trans1->SetAlignGlobal(matrixY);
2454 0 : trans0D->SetAlignGlobalDelta(matrixX);
2455 0 : trans1D->SetAlignGlobalDelta(matrixY);
2456 0 : fitExBTwistX->Init();
2457 0 : fitExBTwistY->Init();
2458 : //
2459 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistX->Clone()),100);
2460 0 : fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwistY->Clone()),101);
2461 : //
2462 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot0->Clone()),102);
2463 0 : fitExBTwistY->AddVisualCorrection((AliTPCExBTwist*)(rot1->Clone()),103);
2464 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot2->Clone()),104);
2465 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(rot3->Clone()),105);
2466 :
2467 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0->Clone()),106);
2468 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1->Clone()),107);
2469 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans0D->Clone()),108);
2470 0 : fitExBTwistX->AddVisualCorrection((AliTPCExBTwist*)(trans1D->Clone()),109);
2471 : //
2472 0 : fAlignTree->SetAlias("FExBTwistX", "AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,100,ptype)+0");
2473 0 : fAlignTree->SetAlias("FExBTwistY","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,101,ptype)+0");
2474 0 : fAlignTree->SetAlias("FAlignRot0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,102,ptype)+0");
2475 0 : fAlignTree->SetAlias("FAlignRot0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,105,ptype)+0");
2476 0 : fAlignTree->SetAlias("FAlignRot1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,103,ptype)+0");
2477 0 : fAlignTree->SetAlias("FAlignRot2","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,104,ptype)+0");
2478 0 : fAlignTree->SetAlias("FAlignTrans0","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,106,ptype)+0");
2479 0 : fAlignTree->SetAlias("FAlignTrans1","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,107,ptype)+0");
2480 0 : fAlignTree->SetAlias("FAlignTrans0D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,108,ptype)+0");
2481 0 : fAlignTree->SetAlias("FAlignTrans1D","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,109,ptype)+0");
2482 : //
2483 : // test fast function
2484 : //
2485 : // fAlignTree->Draw("FExBTwistX:ExBTwistX","isITS&&ptype==0&&abs(snp)<0.05","");
2486 : // fAlignTree->Draw("FExBTwistY:ExBTwistY","isITS&&ptype==0&&abs(snp)<0.05","");
2487 : // fAlignTree->Draw("FAlignRot0:alignRot0","isITS&&ptype==0&&abs(snp)<0.05","");
2488 : // fAlignTree->Draw("FAlignRot1:alignRot1","isITS&&ptype==0&&abs(snp)<0.05","");
2489 : // fAlignTree->Draw("FAlignRot2:alignRot2","isITS&&ptype==0&&abs(snp)<0.05","");
2490 : // //
2491 : // fAlignTree->Draw("FAlignTrans0:alignTrans0","isITS&&ptype==0&&abs(snp)<0.05","");
2492 : // fAlignTree->Draw("FAlignTrans1:alignTrans1","isITS&&ptype==0&&abs(snp)<0.05","");
2493 :
2494 0 : }
2495 :
2496 : //_____________________________________________________________________________
2497 : void AliTPCPreprocessorOffline::CreateAlignTime(TString fstring, TVectorD paramC){
2498 : //
2499 : //
2500 : //
2501 : //
2502 0 : TGeoHMatrix *matrixDelta = new TGeoHMatrix;
2503 0 : TGeoHMatrix *matrixGlobal = new TGeoHMatrix;
2504 0 : Double_t rAngles[9];
2505 : Int_t index=0;
2506 : //
2507 0 : index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans0D");
2508 0 : if (index>=0) matrixDelta->SetDx(paramC[index+1]*0.1);
2509 0 : index=TStatToolkit::GetFitIndex(fstring,"FAlignTrans1D");
2510 0 : if (index>=0) matrixDelta->SetDy(paramC[index+1]*0.1);
2511 0 : rAngles[0]=1; rAngles[4]=1; rAngles[8]=1;
2512 0 : index=TStatToolkit::GetFitIndex(fstring,"FAlignRot0D");
2513 0 : rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001;
2514 0 : rAngles[5]=0; rAngles[7] =0;
2515 0 : rAngles[2]=0; rAngles[6] =0;
2516 0 : matrixDelta->SetRotation(rAngles);
2517 : //
2518 : //
2519 : //
2520 0 : index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans0");
2521 0 : if (index>=0) matrixGlobal->SetDx(paramC[index+1]*0.1);
2522 0 : index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignTrans1");
2523 0 : if (index>=0) matrixGlobal->SetDy(paramC[index+1]*0.1);
2524 0 : rAngles[0]=1; rAngles[4]=1; rAngles[8]=1;
2525 0 : index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot0");
2526 0 : rAngles[1]=-paramC[index+1]*0.001; rAngles[3]=paramC[index+1]*0.001;
2527 0 : index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot1");
2528 0 : rAngles[5]=-paramC[index+1]*0.001; rAngles[7]=paramC[index+1]*0.001;
2529 0 : index=TStatToolkit::GetFitIndex(fstring,"hasITS*FAlignRot2");
2530 0 : rAngles[2]=paramC[index+1]*0.001; rAngles[6] =-paramC[index+1]*0.001;
2531 0 : matrixGlobal->SetRotation(rAngles);
2532 : //
2533 : AliTPCCalibGlobalMisalignment *fitAlignTime =0;
2534 0 : fitAlignTime =new AliTPCCalibGlobalMisalignment;
2535 0 : fitAlignTime->SetName("FitAlignTime");
2536 0 : fitAlignTime->SetTitle("FitAlignTime");
2537 0 : fitAlignTime->SetAlignGlobalDelta(matrixDelta);
2538 0 : fitAlignTime->SetAlignGlobal(matrixGlobal);
2539 : //
2540 0 : AliTPCExBTwist * fitExBTwist= new AliTPCExBTwist;
2541 0 : Int_t indexX=TStatToolkit::GetFitIndex(fstring,"ExBTwistX");
2542 0 : Int_t indexY=TStatToolkit::GetFitIndex(fstring,"ExBTwistY");
2543 0 : fitExBTwist->SetXTwist(0.001*paramC[indexX+1]); // 1 mrad twist in x
2544 0 : fitExBTwist->SetYTwist(0.001*paramC[indexY+1]); // 1 mrad twist in x
2545 0 : fitExBTwist->SetName("FitExBTwistTime");
2546 0 : fitExBTwist->SetTitle("FitExBTwistTime");
2547 0 : AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
2548 0 : Double_t bzField=AliTrackerBase::GetBz();
2549 0 : Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us] // From dataBase: to be updated: per second (ideally)
2550 :
2551 : Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
2552 0 : Double_t wt = -10.0 * (bzField) * vdrift / ezField ;
2553 : //
2554 0 : fitExBTwist->SetOmegaTauT1T2(wt,1,1);
2555 0 : fitExBTwist->Init();
2556 :
2557 0 : AliTPCComposedCorrection *corrTime = new AliTPCComposedCorrection;
2558 0 : TObjArray *arr = new TObjArray;
2559 0 : corrTime->SetCorrections(arr);
2560 :
2561 0 : corrTime->GetCorrections()->Add(fitExBTwist);
2562 0 : corrTime->GetCorrections()->Add(fitAlignTime);
2563 0 : corrTime->SetName("FitCorrectionTime");
2564 0 : corrTime->SetTitle("FitCorrectionTime");
2565 :
2566 0 : fitExBTwist->AddVisualCorrection((AliTPCExBTwist*)(fitExBTwist->Clone()),1001);
2567 0 : fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(fitAlignTime->Clone()),1002);
2568 0 : fitAlignTime->AddVisualCorrection((AliTPCExBTwist*)(corrTime->Clone()),1003);
2569 :
2570 :
2571 0 : fAlignTree->SetAlias("ExBTwistTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1001,ptype)+0");
2572 0 : fAlignTree->SetAlias("AlignTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1002,ptype)+0");
2573 0 : fAlignTree->SetAlias("FitCorrectionTime","AliTPCPreprocessorOffline::EvalAt(phi,refX,theta,1003,ptype)+0");
2574 :
2575 :
2576 0 : TFile *f = new TFile("fitITSVertex.root","update");
2577 0 : corrTime->Write("FitCorrectionTime");
2578 0 : f->Close();
2579 0 : }
2580 :
2581 : //_____________________________________________________________________________
2582 : Int_t AliTPCPreprocessorOffline::GetStatus()
2583 : {
2584 : //get the calibration status
2585 : // 0 means OK
2586 : // positive numbers invalidate for unknown reasons.
2587 : // negative numbers invalidate with a known reason (e.g. low statistics).
2588 : // the returned integer has one bit set for every component that failed.
2589 :
2590 0 : Bool_t enoughStatistics = (fNtracksVdrift>fMinTracksVdrift && fNeventsVdrift>fMinEventsVdrift);
2591 :
2592 0 : if (!enoughStatistics)
2593 : {
2594 0 : fCalibrationStatus=-TMath::Abs(fCalibrationStatus);
2595 0 : }
2596 :
2597 0 : return fCalibrationStatus;
2598 : }
2599 :
2600 : //_____________________________________________________________________________
2601 : void AliTPCPreprocessorOffline::FillQA(Bool_t qa, Bool_t norm/*=kTRUE*/)
2602 : {
2603 : // setup QA histogram array
2604 0 : if (qa && !fArrQAhist) {
2605 0 : fArrQAhist = new TObjArray;
2606 0 : fArrQAhist->SetOwner();
2607 0 : } else {
2608 0 : delete fArrQAhist;
2609 0 : fArrQAhist = 0x0;
2610 : }
2611 :
2612 0 : fNormaliseQA=norm;
2613 0 : }
2614 :
2615 : //_____________________________________________________________________________
2616 : void AliTPCPreprocessorOffline::MakeQAPlotsGain(TString outputDirectory/*=""*/, TString fileTypes/*="png"*/)
2617 : {
2618 : // Draw QA histograms, one per file
2619 : // if outputDirectory is non empty, the QA canvases will be written to the output directory
2620 : // fileTypes can contain output formats, comma separated. Default is one png per canvas.
2621 : // also recognized jpg, gif, root.
2622 : // in case of root, all canvases are written to one root file
2623 0 : if (!fArrQAhist) return;
2624 :
2625 0 : TDirectory *dir=gDirectory;
2626 : TFile *f=0x0;
2627 0 : if (fileTypes.Contains("root")) {
2628 0 : f=new TFile(TString::Format("%s/GainQA.root", outputDirectory.Data()),"recreate");
2629 0 : }
2630 :
2631 : const Int_t ntypes=5;
2632 0 : const TString ftypes[ntypes]={"png","jpg","gif","pdf","eps"};
2633 :
2634 0 : for (Int_t ihist=0; ihist<fArrQAhist->GetEntriesFast(); ++ihist) {
2635 0 : dir->cd();
2636 0 : TH2 *h = static_cast<TH2*>(fArrQAhist->UncheckedAt(ihist));
2637 0 : if (!h) continue;
2638 0 : TString histName = h->GetName();
2639 0 : TString canvName=histName;
2640 0 : canvName.Prepend("c_");
2641 0 : TCanvas *c = static_cast<TCanvas*>(gROOT->GetListOfCanvases()->FindObject(canvName));
2642 0 : if (!c) {
2643 0 : c=new TCanvas(canvName, canvName, 700,500);
2644 0 : }
2645 0 : c->Clear();
2646 0 : c->SetLogz();
2647 :
2648 0 : h->Draw("colz");
2649 :
2650 : // ---| check for derived histogram and draw it on top |---
2651 0 : histName.ReplaceAll("_QA","");
2652 0 : TObject *o = fGainArray->FindObject(histName);
2653 0 : if (o) {
2654 0 : o->Draw("same");
2655 : }
2656 :
2657 : // ---| save output |---
2658 0 : if (!outputDirectory.IsNull()) {
2659 : // --- loop over file types ---
2660 0 : for (Int_t itype=0; itype<ntypes; ++itype) {
2661 0 : if (fileTypes.Contains(ftypes[itype])) {
2662 0 : c->SaveAs(TString::Format("%s/%s.%s", outputDirectory.Data(), c->GetName(), ftypes[itype].Data()));
2663 0 : }
2664 : }
2665 :
2666 : // --- save to file if requested ---
2667 0 : if (f) {
2668 0 : f->cd();
2669 0 : c->Write();
2670 0 : dir->cd();
2671 : }
2672 : }
2673 0 : }
2674 :
2675 0 : delete f;
2676 0 : }
2677 :
2678 :
2679 : void AliTPCPreprocessorOffline::ScaleY(TGraphErrors *graph, Double_t normval)
2680 : {
2681 0 : for (Int_t ipoint=0; ipoint<graph->GetN(); ipoint++) {
2682 0 : graph->GetY()[ipoint]*=normval;
2683 0 : graph->GetEY()[ipoint]*=normval;
2684 : }
2685 0 : }
2686 :
2687 : Bool_t AliTPCPreprocessorOffline::NormaliseYToMean(TGraphErrors *graph)
2688 : {
2689 0 : const Double_t mean = TMath::Mean(graph->GetN(), graph->GetY());
2690 :
2691 0 : if (mean<=0){
2692 0 : AliError(Form("mean=%f",mean));
2693 0 : return kFALSE;
2694 : }
2695 :
2696 0 : ScaleY(graph, 1./mean);
2697 0 : }
2698 :
2699 : Bool_t AliTPCPreprocessorOffline::NormaliseYToWeightedMeandEdx(TGraphErrors *graph)
2700 : {
2701 0 : const Double_t *y=graph->GetY();
2702 0 : Double_t scaleFactor=(63.*y[0] + 64.*y[1] + 32*y[2])/159.;
2703 :
2704 0 : if (scaleFactor<=0){
2705 0 : AliError(Form("scaleFactor=%f",scaleFactor));
2706 0 : return kFALSE;
2707 : }
2708 :
2709 0 : ScaleY(graph, 1./scaleFactor);
2710 0 : }
2711 :
2712 : Bool_t AliTPCPreprocessorOffline::NormaliseYToTruncateddEdx(TGraphErrors *graph)
2713 : {
2714 0 : const Double_t *y=graph->GetY();
2715 0 : Double_t truncationFactor=AliTPCcalibGainMult::GetTruncatedMeanPosition(y[0],y[1],y[2],1000);
2716 :
2717 0 : if (truncationFactor<=0){
2718 0 : AliError(Form("truncationFactor=%f",truncationFactor));
2719 0 : return kFALSE;
2720 : }
2721 :
2722 0 : ScaleY(graph, 1./truncationFactor);
2723 0 : }
2724 :
2725 : /*
2726 : Short sequence to acces the calbration entry:
2727 : TFile *f = TFile::Open("CalibObjects.root");
2728 : AliTPCcalibGainMult * fGainMult = (AliTPCcalibGainMult *)f->Get("TPCCalib/calibGainMult");
2729 :
2730 :
2731 : */
|
