]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDCalibraFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix Coverity defects
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraFit.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18/////////////////////////////////////////////////////////////////////////////////
19//
20// AliTRDCalibraFit
21//
22// This class is for the TRD calibration of the relative gain factor, the drift velocity,
23// the time 0 and the pad response function. It fits the histos.
24// The 2D histograms or vectors (first converted in 1D histos) will be fitted
25// if they have enough entries, otherwise the (default) value of the choosen database
26// will be put. For the relative gain calibration the resulted factors will be globally
27// normalized to the gain factors of the choosen database. It unables to precise
28// previous calibration procedure.
29// The function SetDebug enables the user to see:
30// _fDebug = 0: nothing, only the values are written in the tree if wanted
31// _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32// _fDebug = 2: a comparaison of the coefficients found and the default values
33// in the choosen database.
34// fCoef , histogram of the coefs as function of the calibration group number
35// fDelta , histogram of the relative difference of the coef with the default
36// value in the database as function of the calibration group number
37// fError , dirstribution of this relative difference
38// _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39// pad row and col number
40// _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41// also the comparaison histograms of the 1 for this detector
42//
43//
44// Author:
45// R. Bailhache (R.Bailhache@gsi.de)
46//
47//////////////////////////////////////////////////////////////////////////////////////
48
49#include <TLine.h>
50#include <TH1I.h>
51#include <TStyle.h>
52#include <TProfile2D.h>
53#include <TCanvas.h>
54#include <TGraphErrors.h>
55#include <TObjArray.h>
56#include <TH1.h>
57#include <TH1F.h>
58#include <TF1.h>
59#include <TAxis.h>
60#include <TMath.h>
61#include <TDirectory.h>
62#include <TTreeStream.h>
63#include <TLinearFitter.h>
64#include <TVectorD.h>
65#include <TROOT.h>
66#include <TString.h>
67
68#include "AliLog.h"
69#include "AliMathBase.h"
70
71#include "AliTRDCalibraFit.h"
72#include "AliTRDCalibraMode.h"
73#include "AliTRDCalibraVector.h"
74#include "AliTRDCalibraVdriftLinearFit.h"
75#include "AliTRDcalibDB.h"
76#include "AliTRDgeometry.h"
77#include "AliTRDpadPlane.h"
78#include "AliTRDgeometry.h"
79#include "AliTRDCommonParam.h"
80#include "./Cal/AliTRDCalROC.h"
81#include "./Cal/AliTRDCalPad.h"
82#include "./Cal/AliTRDCalDet.h"
83
84
85ClassImp(AliTRDCalibraFit)
86
87AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
89
90//_____________singleton implementation_________________________________________________
91AliTRDCalibraFit *AliTRDCalibraFit::Instance()
92{
93 //
94 // Singleton implementation
95 //
96
97 if (fgTerminated != kFALSE) {
98 return 0;
99 }
100
101 if (fgInstance == 0) {
102 fgInstance = new AliTRDCalibraFit();
103 }
104
105 return fgInstance;
106
107}
108//______________________________________________________________________________________
109void AliTRDCalibraFit::Terminate()
110{
111 //
112 // Singleton implementation
113 // Deletes the instance of this class
114 //
115
116 fgTerminated = kTRUE;
117
118 if (fgInstance != 0) {
119 delete fgInstance;
120 fgInstance = 0;
121 }
122
123}
124//______________________________________________________________________________________
125AliTRDCalibraFit::AliTRDCalibraFit()
126 :TObject()
127 ,fGeo(0)
128 ,fNumberOfBinsExpected(0)
129 ,fMethod(0)
130 ,fBeginFitCharge(3.5)
131 ,fFitPHPeriode(1)
132 ,fTakeTheMaxPH(kTRUE)
133 ,fT0Shift0(0.124797)
134 ,fT0Shift1(0.267451)
135 ,fRangeFitPRF(1.0)
136 ,fAccCDB(kFALSE)
137 ,fMinEntries(800)
138 ,fRebin(1)
139 ,fNumberFit(0)
140 ,fNumberFitSuccess(0)
141 ,fNumberEnt(0)
142 ,fStatisticMean(0.0)
143 ,fDebugStreamer(0x0)
144 ,fDebugLevel(0)
145 ,fFitVoir(0)
146 ,fMagneticField(0.5)
147 ,fCalibraMode(new AliTRDCalibraMode())
148 ,fCurrentCoefE(0.0)
149 ,fCurrentCoefE2(0.0)
150 ,fDect1(0)
151 ,fDect2(0)
152 ,fScaleFitFactor(0.0)
153 ,fEntriesCurrent(0)
154 ,fCountDet(0)
155 ,fCount(0)
156 ,fNbDet(0)
157 ,fCalDet(0x0)
158 ,fCalROC(0x0)
159 ,fCalDet2(0x0)
160 ,fCalROC2(0x0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
163 ,fVectorFit(0)
164 ,fVectorFit2(0)
165{
166 //
167 // Default constructor
168 //
169
170 fGeo = new AliTRDgeometry();
171
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
176 }
177 for (Int_t i = 0; i < 3; i++) {
178 fPhd[i] = 0.0;
179 fDet[i] = 0;
180 }
181
182}
183//______________________________________________________________________________________
184AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
185:TObject(c)
186,fGeo(0)
187,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
188,fMethod(c.fMethod)
189,fBeginFitCharge(c.fBeginFitCharge)
190,fFitPHPeriode(c.fFitPHPeriode)
191,fTakeTheMaxPH(c.fTakeTheMaxPH)
192,fT0Shift0(c.fT0Shift0)
193,fT0Shift1(c.fT0Shift1)
194,fRangeFitPRF(c.fRangeFitPRF)
195,fAccCDB(c.fAccCDB)
196,fMinEntries(c.fMinEntries)
197,fRebin(c.fRebin)
198,fNumberFit(c.fNumberFit)
199,fNumberFitSuccess(c.fNumberFitSuccess)
200,fNumberEnt(c.fNumberEnt)
201,fStatisticMean(c.fStatisticMean)
202,fDebugStreamer(0x0)
203,fDebugLevel(c.fDebugLevel)
204,fFitVoir(c.fFitVoir)
205,fMagneticField(c.fMagneticField)
206,fCalibraMode(0x0)
207,fCurrentCoefE(c.fCurrentCoefE)
208,fCurrentCoefE2(c.fCurrentCoefE2)
209,fDect1(c.fDect1)
210,fDect2(c.fDect2)
211,fScaleFitFactor(c.fScaleFitFactor)
212,fEntriesCurrent(c.fEntriesCurrent)
213,fCountDet(c.fCountDet)
214,fCount(c.fCount)
215,fNbDet(c.fNbDet)
216,fCalDet(0x0)
217,fCalROC(0x0)
218,fCalDet2(0x0)
219,fCalROC2(0x0)
220,fCurrentCoefDetector(0x0)
221,fCurrentCoefDetector2(0x0)
222,fVectorFit(0)
223,fVectorFit2(0)
224{
225 //
226 // Copy constructor
227 //
228
229 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
230
231 //Current variables initialised
232 for (Int_t k = 0; k < 2; k++) {
233 fCurrentCoef[k] = 0.0;
234 fCurrentCoef2[k] = 0.0;
235 }
236 for (Int_t i = 0; i < 3; i++) {
237 fPhd[i] = 0.0;
238 fDet[i] = 0;
239 }
240 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
241 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
242
243 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
244 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
245
246 fVectorFit.SetName(c.fVectorFit.GetName());
247 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
248 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
249 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
250 Int_t ntotal = 1;
251 if (GetStack(detector) == 2) {
252 ntotal = 1728;
253 }
254 else {
255 ntotal = 2304;
256 }
257 Float_t *coef = new Float_t[ntotal];
258 for (Int_t i = 0; i < ntotal; i++) {
259 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
260 }
261 fitInfo->SetCoef(coef);
262 fitInfo->SetDetector(detector);
263 fVectorFit.Add((TObject *) fitInfo);
264 }
265 fVectorFit.SetName(c.fVectorFit.GetName());
266 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
267 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
268 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
269 Int_t ntotal = 1;
270 if (GetStack(detector) == 2) {
271 ntotal = 1728;
272 }
273 else {
274 ntotal = 2304;
275 }
276 Float_t *coef = new Float_t[ntotal];
277 for (Int_t i = 0; i < ntotal; i++) {
278 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
279 }
280 fitInfo->SetCoef(coef);
281 fitInfo->SetDetector(detector);
282 fVectorFit2.Add((TObject *) fitInfo);
283 }
284 if (fGeo) {
285 delete fGeo;
286 }
287 fGeo = new AliTRDgeometry();
288
289}
290//____________________________________________________________________________________
291AliTRDCalibraFit::~AliTRDCalibraFit()
292{
293 //
294 // AliTRDCalibraFit destructor
295 //
296 if ( fDebugStreamer ) delete fDebugStreamer;
297 if ( fCalDet ) delete fCalDet;
298 if ( fCalDet2 ) delete fCalDet2;
299 if ( fCalROC ) delete fCalROC;
300 if ( fCalROC2 ) delete fCalROC2;
301 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
302 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
303 fVectorFit.Delete();
304 fVectorFit2.Delete();
305 if (fGeo) {
306 delete fGeo;
307 }
308
309}
310//_____________________________________________________________________________
311void AliTRDCalibraFit::Destroy()
312{
313 //
314 // Delete instance
315 //
316
317 if (fgInstance) {
318 delete fgInstance;
319 fgInstance = 0x0;
320 }
321
322}
323//_____________________________________________________________________________
324void AliTRDCalibraFit::DestroyDebugStreamer()
325{
326 //
327 // Delete DebugStreamer
328 //
329
330 if ( fDebugStreamer ) delete fDebugStreamer;
331 fDebugStreamer = 0x0;
332
333}
334//__________________________________________________________________________________
335void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
336{
337 //
338 // From the drift velocity and t0
339 // return the position of the peak and maximum negative slope
340 //
341
342 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
343 Double_t widbins = 0.1; // 0.1 mus
344
345 //peak and maxnegslope in mus
346 Double_t begind = t0*widbins + fT0Shift0;
347 Double_t peakd = t0*widbins + fT0Shift1;
348 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
349
350 // peak and maxnegslope in timebin
351 begin = TMath::Nint(begind*widbins);
352 peak = TMath::Nint(peakd*widbins);
353 end = TMath::Nint(maxnegslope*widbins);
354
355}
356//____________Functions fit Online CH2d________________________________________
357Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
358{
359 //
360 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
361 // calibration group normalized the resulted coefficients (to 1 normally)
362 //
363
364 // Set the calibration mode
365 //const char *name = ch->GetTitle();
366 TString name = ch->GetTitle();
367 if(!SetModeCalibration(name,0)) return kFALSE;
368
369 // Number of Ybins (detectors or groups of pads)
370 Int_t nbins = ch->GetNbinsX();// charge
371 Int_t nybins = ch->GetNbinsY();// groups number
372 if (!InitFit(nybins,0)) {
373 return kFALSE;
374 }
375 if (!InitFitCH()) {
376 return kFALSE;
377 }
378 fStatisticMean = 0.0;
379 fNumberFit = 0;
380 fNumberFitSuccess = 0;
381 fNumberEnt = 0;
382 // Init fCountDet and fCount
383 InitfCountDetAndfCount(0);
384 // Beginning of the loop betwwen dect1 and dect2
385 for (Int_t idect = fDect1; idect < fDect2; idect++) {
386 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
387 UpdatefCountDetAndfCount(idect,0);
388 ReconstructFitRowMinRowMax(idect, 0);
389 // Take the histo
390 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
391 projch->SetDirectory(0);
392 // Number of entries for this calibration group
393 Double_t nentries = 0.0;
394 Double_t mean = 0.0;
395 for (Int_t k = 0; k < nbins; k++) {
396 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
397 nentries += ch->GetBinContent(binnb);
398 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
399 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
400 }
401 projch->SetEntries(nentries);
402 //printf("The number of entries for the group %d is %f\n",idect,nentries);
403 if (nentries > 0) {
404 fNumberEnt++;
405 mean /= nentries;
406 }
407 // Rebin and statistic stuff
408 if (fRebin > 1) {
409 projch = ReBin((TH1I *) projch);
410 }
411 // This detector has not enough statistics or was off
412 if (nentries <= fMinEntries) {
413 NotEnoughStatisticCH(idect);
414 if (fDebugLevel != 1) {
415 delete projch;
416 }
417 continue;
418 }
419 // Statistics of the group fitted
420 fStatisticMean += nentries;
421 fNumberFit++;
422 //Method choosen
423 switch(fMethod)
424 {
425 case 0: FitMeanW((TH1 *) projch, nentries); break;
426 case 1: FitMean((TH1 *) projch, nentries, mean); break;
427 case 2: FitCH((TH1 *) projch, mean); break;
428 case 3: FitBisCH((TH1 *) projch, mean); break;
429 default: return kFALSE;
430 }
431 // Fill Infos Fit
432 FillInfosFitCH(idect);
433 // Memory!!!
434 if (fDebugLevel != 1) {
435 delete projch;
436 }
437 } // Boucle object
438 // Normierungcharge
439 if (fDebugLevel != 1) {
440 NormierungCharge();
441 }
442 // Mean Statistic
443 if (fNumberFit > 0) {
444 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
445 fStatisticMean = fStatisticMean / fNumberFit;
446 }
447 else {
448 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
449 }
450 delete fDebugStreamer;
451 fDebugStreamer = 0x0;
452
453 return kTRUE;
454}
455//____________Functions fit Online CH2d________________________________________
456Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
457{
458 //
459 // Reconstruct a 1D histo from the vectorCH for each calibration group,
460 // fit the histo, normalized the resulted coefficients (to 1 normally)
461 //
462
463 // Set the calibraMode
464 //const char *name = calvect->GetNameCH();
465 TString name = calvect->GetNameCH();
466 if(!SetModeCalibration(name,0)) return kFALSE;
467
468 // Number of Xbins (detectors or groups of pads)
469 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
470 return kFALSE;
471 }
472 if (!InitFitCH()) {
473 return kFALSE;
474 }
475 fStatisticMean = 0.0;
476 fNumberFit = 0;
477 fNumberFitSuccess = 0;
478 fNumberEnt = 0;
479 // Init fCountDet and fCount
480 InitfCountDetAndfCount(0);
481 // Beginning of the loop between dect1 and dect2
482 for (Int_t idect = fDect1; idect < fDect2; idect++) {
483 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
484 UpdatefCountDetAndfCount(idect,0);
485 ReconstructFitRowMinRowMax(idect,0);
486 // Take the histo
487 Double_t nentries = 0.0;
488 Double_t mean = 0.0;
489 if(!calvect->GetCHEntries(fCountDet)) {
490 NotEnoughStatisticCH(idect);
491 continue;
492 }
493
494 TString tname("CH");
495 tname += idect;
496 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
497 projch->SetDirectory(0);
498 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
499 nentries += projch->GetBinContent(k+1);
500 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
501 }
502 if (nentries > 0) {
503 fNumberEnt++;
504 mean /= nentries;
505 }
506 //printf("The number of entries for the group %d is %f\n",idect,nentries);
507 // Rebin
508 if (fRebin > 1) {
509 projch = ReBin((TH1F *) projch);
510 }
511 // This detector has not enough statistics or was not found in VectorCH
512 if (nentries <= fMinEntries) {
513 NotEnoughStatisticCH(idect);
514 continue;
515 }
516 // Statistic of the histos fitted
517 fStatisticMean += nentries;
518 fNumberFit++;
519 //Method choosen
520 switch(fMethod)
521 {
522 case 0: FitMeanW((TH1 *) projch, nentries); break;
523 case 1: FitMean((TH1 *) projch, nentries, mean); break;
524 case 2: FitCH((TH1 *) projch, mean); break;
525 case 3: FitBisCH((TH1 *) projch, mean); break;
526 default: return kFALSE;
527 }
528 // Fill Infos Fit
529 FillInfosFitCH(idect);
530 } // Boucle object
531 // Normierungcharge
532 if (fDebugLevel != 1) {
533 NormierungCharge();
534 }
535 // Mean Statistics
536 if (fNumberFit > 0) {
537 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
538 fStatisticMean = fStatisticMean / fNumberFit;
539 }
540 else {
541 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
542 }
543 delete fDebugStreamer;
544 fDebugStreamer = 0x0;
545 return kTRUE;
546}
547//________________functions fit Online PH2d____________________________________
548Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
549{
550 //
551 // Take the 1D profiles (average pulse height), projections of the 2D PH
552 // on the Xaxis, for each calibration group
553 // Reconstruct a drift velocity
554 // A first calibration of T0 is also made using the same method
555 //
556
557 // Set the calibration mode
558 //const char *name = ph->GetTitle();
559 TString name = ph->GetTitle();
560 if(!SetModeCalibration(name,1)) return kFALSE;
561
562 //printf("Mode calibration set\n");
563
564 // Number of Xbins (detectors or groups of pads)
565 Int_t nbins = ph->GetNbinsX();// time
566 Int_t nybins = ph->GetNbinsY();// calibration group
567 if (!InitFit(nybins,1)) {
568 return kFALSE;
569 }
570
571 //printf("Init fit\n");
572
573 if (!InitFitPH()) {
574 return kFALSE;
575 }
576
577 //printf("Init fit PH\n");
578
579 fStatisticMean = 0.0;
580 fNumberFit = 0;
581 fNumberFitSuccess = 0;
582 fNumberEnt = 0;
583 // Init fCountDet and fCount
584 InitfCountDetAndfCount(1);
585 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
586
587 // Beginning of the loop
588 for (Int_t idect = fDect1; idect < fDect2; idect++) {
589 //printf("idect = %d\n",idect);
590 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
591 UpdatefCountDetAndfCount(idect,1);
592 ReconstructFitRowMinRowMax(idect,1);
593 // Take the histo
594 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
595 projph->SetDirectory(0);
596 // Number of entries for this calibration group
597 Double_t nentries = 0;
598 for (Int_t k = 0; k < nbins; k++) {
599 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
600 nentries += ph->GetBinEntries(binnb);
601 }
602 if (nentries > 0) {
603 fNumberEnt++;
604 }
605 //printf("The number of entries for the group %d is %f\n",idect,nentries);
606 // This detector has not enough statistics or was off
607 if (nentries <= fMinEntries) {
608 //printf("Not enough statistic!\n");
609 NotEnoughStatisticPH(idect,nentries);
610 if (fDebugLevel != 1) {
611 delete projph;
612 }
613 continue;
614 }
615 // Statistics of the histos fitted
616 fNumberFit++;
617 fStatisticMean += nentries;
618 // Calcul of "real" coef
619 CalculVdriftCoefMean();
620 CalculT0CoefMean();
621 //Method choosen
622 //printf("Method\n");
623 switch(fMethod)
624 {
625 case 0: FitLagrangePoly((TH1 *) projph); break;
626 case 1: FitPente((TH1 *) projph); break;
627 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
628 default: return kFALSE;
629 }
630 // Fill the tree if end of a detector or only the pointer to the branch!!!
631 FillInfosFitPH(idect,nentries);
632 // Memory!!!
633 if (fDebugLevel != 1) {
634 delete projph;
635 }
636 } // Boucle object
637 // Mean Statistic
638 if (fNumberFit > 0) {
639 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
640 fStatisticMean = fStatisticMean / fNumberFit;
641 }
642 else {
643 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
644 }
645 delete fDebugStreamer;
646 fDebugStreamer = 0x0;
647 return kTRUE;
648}
649//____________Functions fit Online PH2d________________________________________
650Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
651{
652 //
653 // Reconstruct the average pulse height from the vectorPH for each
654 // calibration group
655 // Reconstruct a drift velocity
656 // A first calibration of T0 is also made using the same method (slope method)
657 //
658
659 // Set the calibration mode
660 //const char *name = calvect->GetNamePH();
661 TString name = calvect->GetNamePH();
662 if(!SetModeCalibration(name,1)) return kFALSE;
663
664 // Number of Xbins (detectors or groups of pads)
665 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
666 return kFALSE;
667 }
668 if (!InitFitPH()) {
669 return kFALSE;
670 }
671 fStatisticMean = 0.0;
672 fNumberFit = 0;
673 fNumberFitSuccess = 0;
674 fNumberEnt = 0;
675 // Init fCountDet and fCount
676 InitfCountDetAndfCount(1);
677 // Beginning of the loop
678 for (Int_t idect = fDect1; idect < fDect2; idect++) {
679 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
680 UpdatefCountDetAndfCount(idect,1);
681 ReconstructFitRowMinRowMax(idect,1);
682 // Take the histo
683 fEntriesCurrent = 0;
684 if(!calvect->GetPHEntries(fCountDet)) {
685 NotEnoughStatisticPH(idect,fEntriesCurrent);
686 continue;
687 }
688 TString tname("PH");
689 tname += idect;
690 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
691 projph->SetDirectory(0);
692 if(fEntriesCurrent > 0) fNumberEnt++;
693 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
694 // This detector has not enough statistics or was off
695 if (fEntriesCurrent <= fMinEntries) {
696 //printf("Not enough stat!\n");
697 NotEnoughStatisticPH(idect,fEntriesCurrent);
698 continue;
699 }
700 // Statistic of the histos fitted
701 fNumberFit++;
702 fStatisticMean += fEntriesCurrent;
703 // Calcul of "real" coef
704 CalculVdriftCoefMean();
705 CalculT0CoefMean();
706 //Method choosen
707 switch(fMethod)
708 {
709 case 0: FitLagrangePoly((TH1 *) projph); break;
710 case 1: FitPente((TH1 *) projph); break;
711 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
712 default: return kFALSE;
713 }
714 // Fill the tree if end of a detector or only the pointer to the branch!!!
715 FillInfosFitPH(idect,fEntriesCurrent);
716 } // Boucle object
717
718 // Mean Statistic
719 if (fNumberFit > 0) {
720 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
721 fStatisticMean = fStatisticMean / fNumberFit;
722 }
723 else {
724 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
725 }
726 delete fDebugStreamer;
727 fDebugStreamer = 0x0;
728 return kTRUE;
729}
730//____________Functions fit Online PRF2d_______________________________________
731Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
732{
733 //
734 // Take the 1D profiles (pad response function), projections of the 2D PRF
735 // on the Xaxis, for each calibration group
736 // Fit with a gaussian to reconstruct the sigma of the pad response function
737 //
738
739 // Set the calibration mode
740 //const char *name = prf->GetTitle();
741 TString name = prf->GetTitle();
742 if(!SetModeCalibration(name,2)) return kFALSE;
743
744 // Number of Ybins (detectors or groups of pads)
745 Int_t nybins = prf->GetNbinsY();// calibration groups
746 Int_t nbins = prf->GetNbinsX();// bins
747 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
748 if((nbg > 0) || (nbg == -1)) return kFALSE;
749 if (!InitFit(nybins,2)) {
750 return kFALSE;
751 }
752 if (!InitFitPRF()) {
753 return kFALSE;
754 }
755 fStatisticMean = 0.0;
756 fNumberFit = 0;
757 fNumberFitSuccess = 0;
758 fNumberEnt = 0;
759 // Init fCountDet and fCount
760 InitfCountDetAndfCount(2);
761 // Beginning of the loop
762 for (Int_t idect = fDect1; idect < fDect2; idect++) {
763 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
764 UpdatefCountDetAndfCount(idect,2);
765 ReconstructFitRowMinRowMax(idect,2);
766 // Take the histo
767 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
768 projprf->SetDirectory(0);
769 // Number of entries for this calibration group
770 Double_t nentries = 0;
771 for (Int_t k = 0; k < nbins; k++) {
772 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
773 nentries += prf->GetBinEntries(binnb);
774 }
775 if(nentries > 0) fNumberEnt++;
776 // This detector has not enough statistics or was off
777 if (nentries <= fMinEntries) {
778 NotEnoughStatisticPRF(idect);
779 if (fDebugLevel != 1) {
780 delete projprf;
781 }
782 continue;
783 }
784 // Statistics of the histos fitted
785 fNumberFit++;
786 fStatisticMean += nentries;
787 // Calcul of "real" coef
788 CalculPRFCoefMean();
789 //Method choosen
790 switch(fMethod)
791 {
792 case 0: FitPRF((TH1 *) projprf); break;
793 case 1: RmsPRF((TH1 *) projprf); break;
794 default: return kFALSE;
795 }
796 // Fill the tree if end of a detector or only the pointer to the branch!!!
797 FillInfosFitPRF(idect);
798 // Memory!!!
799 if (fDebugLevel != 1) {
800 delete projprf;
801 }
802 } // Boucle object
803 // Mean Statistic
804 if (fNumberFit > 0) {
805 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
806 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
807 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
808 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
809 fStatisticMean = fStatisticMean / fNumberFit;
810 }
811 else {
812 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
813 }
814 delete fDebugStreamer;
815 fDebugStreamer = 0x0;
816 return kTRUE;
817}
818//____________Functions fit Online PRF2d_______________________________________
819Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
820{
821 //
822 // Take the 1D profiles (pad response function), projections of the 2D PRF
823 // on the Xaxis, for each calibration group
824 // Fit with a gaussian to reconstruct the sigma of the pad response function
825 //
826
827 // Set the calibration mode
828 //const char *name = prf->GetTitle();
829 TString name = prf->GetTitle();
830 if(!SetModeCalibration(name,2)) return kFALSE;
831
832 // Number of Ybins (detectors or groups of pads)
833 TAxis *xprf = prf->GetXaxis();
834 TAxis *yprf = prf->GetYaxis();
835 Int_t nybins = yprf->GetNbins();// calibration groups
836 Int_t nbins = xprf->GetNbins();// bins
837 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
838 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
839 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
840 if(nbg == -1) return kFALSE;
841 if(nbg > 0) fMethod = 1;
842 else fMethod = 0;
843 if (!InitFit(nybins,2)) {
844 return kFALSE;
845 }
846 if (!InitFitPRF()) {
847 return kFALSE;
848 }
849 fStatisticMean = 0.0;
850 fNumberFit = 0;
851 fNumberFitSuccess = 0;
852 fNumberEnt = 0;
853 // Init fCountDet and fCount
854 InitfCountDetAndfCount(2);
855 // Beginning of the loop
856 for (Int_t idect = fDect1; idect < fDect2; idect++) {
857 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
858 UpdatefCountDetAndfCount(idect,2);
859 ReconstructFitRowMinRowMax(idect,2);
860 // Build the array of entries and sum
861 TArrayD arraye = TArrayD(nbins);
862 TArrayD arraym = TArrayD(nbins);
863 TArrayD arrayme = TArrayD(nbins);
864 Double_t nentries = 0;
865 //printf("nbins %d\n",nbins);
866 for (Int_t k = 0; k < nbins; k++) {
867 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
868 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
869 Double_t mean = (Double_t)prf->GetBinContent(binnb);
870 Double_t error = (Double_t)prf->GetBinError(binnb);
871 //printf("for %d we have %f\n",k,entries);
872 nentries += entries;
873 arraye.AddAt(entries,k);
874 arraym.AddAt(mean,k);
875 arrayme.AddAt(error,k);
876 }
877 if(nentries > 0) fNumberEnt++;
878 //printf("The number of entries for the group %d is %f\n",idect,nentries);
879 // This detector has not enough statistics or was off
880 if (nentries <= fMinEntries) {
881 NotEnoughStatisticPRF(idect);
882 continue;
883 }
884 // Statistics of the histos fitted
885 fNumberFit++;
886 fStatisticMean += nentries;
887 // Calcul of "real" coef
888 CalculPRFCoefMean();
889 //Method choosen
890 switch(fMethod)
891 {
892 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
893 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
894 default: return kFALSE;
895 }
896 // Fill the tree if end of a detector or only the pointer to the branch!!!
897 FillInfosFitPRF(idect);
898 } // Boucle object
899 // Mean Statistic
900 if (fNumberFit > 0) {
901 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
902 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
903 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
904 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
905 fStatisticMean = fStatisticMean / fNumberFit;
906 }
907 else {
908 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
909 }
910 delete fDebugStreamer;
911 fDebugStreamer = 0x0;
912 return kTRUE;
913}
914//____________Functions fit Online PRF2d_______________________________________
915Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
916{
917 //
918 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
919 // each calibration group
920 // Fit with a gaussian to reconstruct the sigma of the pad response function
921 //
922
923 // Set the calibra mode
924 //const char *name = calvect->GetNamePRF();
925 TString name = calvect->GetNamePRF();
926 if(!SetModeCalibration(name,2)) return kFALSE;
927 //printf("test0 %s\n",name);
928
929 // Number of Xbins (detectors or groups of pads)
930 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
931 //printf("test1\n");
932 return kFALSE;
933 }
934 if (!InitFitPRF()) {
935 ///printf("test2\n");
936 return kFALSE;
937 }
938 fStatisticMean = 0.0;
939 fNumberFit = 0;
940 fNumberFitSuccess = 0;
941 fNumberEnt = 0;
942 // Init fCountDet and fCount
943 InitfCountDetAndfCount(2);
944 // Beginning of the loop
945 for (Int_t idect = fDect1; idect < fDect2; idect++) {
946 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
947 UpdatefCountDetAndfCount(idect,2);
948 ReconstructFitRowMinRowMax(idect,2);
949 // Take the histo
950 fEntriesCurrent = 0;
951 if(!calvect->GetPRFEntries(fCountDet)) {
952 NotEnoughStatisticPRF(idect);
953 continue;
954 }
955 TString tname("PRF");
956 tname += idect;
957 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
958 projprf->SetDirectory(0);
959 if(fEntriesCurrent > 0) fNumberEnt++;
960 // This detector has not enough statistics or was off
961 if (fEntriesCurrent <= fMinEntries) {
962 NotEnoughStatisticPRF(idect);
963 continue;
964 }
965 // Statistic of the histos fitted
966 fNumberFit++;
967 fStatisticMean += fEntriesCurrent;
968 // Calcul of "real" coef
969 CalculPRFCoefMean();
970 //Method choosen
971 switch(fMethod)
972 {
973 case 1: FitPRF((TH1 *) projprf); break;
974 case 2: RmsPRF((TH1 *) projprf); break;
975 default: return kFALSE;
976 }
977 // Fill the tree if end of a detector or only the pointer to the branch!!!
978 FillInfosFitPRF(idect);
979 } // Boucle object
980 // Mean Statistics
981 if (fNumberFit > 0) {
982 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
983 }
984 else {
985 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
986 }
987 delete fDebugStreamer;
988 fDebugStreamer = 0x0;
989 return kTRUE;
990}
991//____________Functions fit Online PRF2d_______________________________________
992Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
993{
994 //
995 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
996 // each calibration group
997 // Fit with a gaussian to reconstruct the sigma of the pad response function
998 //
999
1000 // Set the calibra mode
1001 //const char *name = calvect->GetNamePRF();
1002 TString name = calvect->GetNamePRF();
1003 if(!SetModeCalibration(name,2)) return kFALSE;
1004 //printf("test0 %s\n",name);
1005 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1006 //printf("test1 %d\n",nbg);
1007 if(nbg == -1) return kFALSE;
1008 if(nbg > 0) fMethod = 1;
1009 else fMethod = 0;
1010 // Number of Xbins (detectors or groups of pads)
1011 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1012 //printf("test2\n");
1013 return kFALSE;
1014 }
1015 if (!InitFitPRF()) {
1016 //printf("test3\n");
1017 return kFALSE;
1018 }
1019 fStatisticMean = 0.0;
1020 fNumberFit = 0;
1021 fNumberFitSuccess = 0;
1022 fNumberEnt = 0;
1023 // Variables
1024 Int_t nbins = 0;
1025 Double_t *arrayx = 0;
1026 Double_t *arraye = 0;
1027 Double_t *arraym = 0;
1028 Double_t *arrayme = 0;
1029 Float_t lowedge = 0.0;
1030 Float_t upedge = 0.0;
1031 // Init fCountDet and fCount
1032 InitfCountDetAndfCount(2);
1033 // Beginning of the loop
1034 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1035 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1036 UpdatefCountDetAndfCount(idect,2);
1037 ReconstructFitRowMinRowMax(idect,2);
1038 // Take the histo
1039 fEntriesCurrent = 0;
1040 if(!calvect->GetPRFEntries(fCountDet)) {
1041 NotEnoughStatisticPRF(idect);
1042 continue;
1043 }
1044 TString tname("PRF");
1045 tname += idect;
1046 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1047 nbins = projprftree->GetN();
1048 arrayx = (Double_t *)projprftree->GetX();
1049 arraye = (Double_t *)projprftree->GetEX();
1050 arraym = (Double_t *)projprftree->GetY();
1051 arrayme = (Double_t *)projprftree->GetEY();
1052 Float_t step = arrayx[1]-arrayx[0];
1053 lowedge = arrayx[0] - step/2.0;
1054 upedge = arrayx[(nbins-1)] + step/2.0;
1055 //printf("nbins est %d\n",nbins);
1056 for(Int_t k = 0; k < nbins; k++){
1057 fEntriesCurrent += (Int_t)arraye[k];
1058 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1059 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1060 }
1061 if(fEntriesCurrent > 0) fNumberEnt++;
1062 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1063 // This detector has not enough statistics or was off
1064 if (fEntriesCurrent <= fMinEntries) {
1065 NotEnoughStatisticPRF(idect);
1066 continue;
1067 }
1068 // Statistic of the histos fitted
1069 fNumberFit++;
1070 fStatisticMean += fEntriesCurrent;
1071 // Calcul of "real" coef
1072 CalculPRFCoefMean();
1073 //Method choosen
1074 switch(fMethod)
1075 {
1076 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1077 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1078 default: return kFALSE;
1079 }
1080 // Fill the tree if end of a detector or only the pointer to the branch!!!
1081 FillInfosFitPRF(idect);
1082 } // Boucle object
1083 // Mean Statistics
1084 if (fNumberFit > 0) {
1085 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1086 }
1087 else {
1088 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1089 }
1090 delete fDebugStreamer;
1091 fDebugStreamer = 0x0;
1092 return kTRUE;
1093}
1094//____________Functions fit Online CH2d________________________________________
1095Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1096{
1097 //
1098 // The linear method
1099 //
1100
1101 fStatisticMean = 0.0;
1102 fNumberFit = 0;
1103 fNumberFitSuccess = 0;
1104 fNumberEnt = 0;
1105 if(!InitFitLinearFitter()) return kFALSE;
1106
1107
1108 for(Int_t idet = 0; idet < 540; idet++){
1109
1110
1111 //printf("detector number %d\n",idet);
1112
1113 // Take the result
1114 TVectorD param(2);
1115 TVectorD error(3);
1116 fEntriesCurrent = 0;
1117 fCountDet = idet;
1118 Bool_t here = calivdli->GetParam(idet,&param);
1119 Bool_t heree = calivdli->GetError(idet,&error);
1120 //printf("here %d and heree %d\n",here, heree);
1121 if(heree) {
1122 fEntriesCurrent = (Int_t) error[2];
1123 fNumberEnt++;
1124 }
1125 //printf("Number of entries %d\n",fEntriesCurrent);
1126 // Nothing found or not enough statistic
1127 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1128 NotEnoughStatisticLinearFitter();
1129 continue;
1130 }
1131 //param.Print();
1132 //error.Print();
1133 //Statistics
1134 fNumberFit++;
1135 fStatisticMean += fEntriesCurrent;
1136
1137 // Check the fit
1138 if((-(param[1])) <= 0.0) {
1139 NotEnoughStatisticLinearFitter();
1140 continue;
1141 }
1142
1143 // CalculDatabaseVdriftandTan
1144 CalculVdriftLorentzCoef();
1145
1146 // Statistics
1147 fNumberFitSuccess ++;
1148
1149 // Put the fCurrentCoef
1150 fCurrentCoef[0] = -param[1];
1151 // here the database must be the one of the reconstruction for the lorentz angle....
1152 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1153 fCurrentCoefE = error[1];
1154 fCurrentCoefE2 = error[0];
1155 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1156 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1157 }
1158
1159 // Fill
1160 FillInfosFitLinearFitter();
1161
1162
1163 }
1164 // Mean Statistics
1165 if (fNumberFit > 0) {
1166 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1167 }
1168 else {
1169 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1170 }
1171 delete fDebugStreamer;
1172 fDebugStreamer = 0x0;
1173 return kTRUE;
1174
1175}
1176//____________Functions for seeing if the pad is really okey___________________
1177//_____________________________________________________________________________
1178Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1179{
1180 //
1181 // Get numberofgroupsprf
1182 //
1183
1184 // Some patterns
1185 const Char_t *pattern0 = "Ngp0";
1186 const Char_t *pattern1 = "Ngp1";
1187 const Char_t *pattern2 = "Ngp2";
1188 const Char_t *pattern3 = "Ngp3";
1189 const Char_t *pattern4 = "Ngp4";
1190 const Char_t *pattern5 = "Ngp5";
1191 const Char_t *pattern6 = "Ngp6";
1192
1193 // Nrphi mode
1194 if (strstr(nametitle.Data(),pattern0)) {
1195 return 0;
1196 }
1197 if (strstr(nametitle.Data(),pattern1)) {
1198 return 1;
1199 }
1200 if (strstr(nametitle.Data(),pattern2)) {
1201 return 2;
1202 }
1203 if (strstr(nametitle.Data(),pattern3)) {
1204 return 3;
1205 }
1206 if (strstr(nametitle.Data(),pattern4)) {
1207 return 4;
1208 }
1209 if (strstr(nametitle.Data(),pattern5)) {
1210 return 5;
1211 }
1212 if (strstr(nametitle.Data(),pattern6)){
1213 return 6;
1214 }
1215 else return -1;
1216
1217
1218}
1219//_____________________________________________________________________________
1220Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1221{
1222 //
1223 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1224 // corresponding to the given name
1225 //
1226
1227 if(!SetNzFromTObject(name,i)) return kFALSE;
1228 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1229
1230 return kTRUE;
1231
1232}
1233//_____________________________________________________________________________
1234Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1235{
1236 //
1237 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1238 // corresponding to the given TObject
1239 //
1240
1241 // Some patterns
1242 const Char_t *patternrphi0 = "Nrphi0";
1243 const Char_t *patternrphi1 = "Nrphi1";
1244 const Char_t *patternrphi2 = "Nrphi2";
1245 const Char_t *patternrphi3 = "Nrphi3";
1246 const Char_t *patternrphi4 = "Nrphi4";
1247 const Char_t *patternrphi5 = "Nrphi5";
1248 const Char_t *patternrphi6 = "Nrphi6";
1249
1250
1251 const Char_t *patternrphi10 = "Nrphi10";
1252 const Char_t *patternrphi100 = "Nrphi100";
1253 const Char_t *patternz10 = "Nz10";
1254 const Char_t *patternz100 = "Nz100";
1255
1256 // Nrphi mode
1257 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1258 fCalibraMode->SetAllTogether(i);
1259 fNbDet = 540;
1260 if (fDebugLevel > 1) {
1261 AliInfo(Form("fNbDet %d and 100",fNbDet));
1262 }
1263 return kTRUE;
1264 }
1265 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1266 fCalibraMode->SetPerSuperModule(i);
1267 fNbDet = 30;
1268 if (fDebugLevel > 1) {
1269 AliInfo(Form("fNDet %d and 100",fNbDet));
1270 }
1271 return kTRUE;
1272 }
1273
1274 if (strstr(name.Data(),patternrphi0)) {
1275 fCalibraMode->SetNrphi(i ,0);
1276 if (fDebugLevel > 1) {
1277 AliInfo(Form("fNbDet %d and 0",fNbDet));
1278 }
1279 return kTRUE;
1280 }
1281 if (strstr(name.Data(),patternrphi1)) {
1282 fCalibraMode->SetNrphi(i, 1);
1283 if (fDebugLevel > 1) {
1284 AliInfo(Form("fNbDet %d and 1",fNbDet));
1285 }
1286 return kTRUE;
1287 }
1288 if (strstr(name.Data(),patternrphi2)) {
1289 fCalibraMode->SetNrphi(i, 2);
1290 if (fDebugLevel > 1) {
1291 AliInfo(Form("fNbDet %d and 2",fNbDet));
1292 }
1293 return kTRUE;
1294 }
1295 if (strstr(name.Data(),patternrphi3)) {
1296 fCalibraMode->SetNrphi(i, 3);
1297 if (fDebugLevel > 1) {
1298 AliInfo(Form("fNbDet %d and 3",fNbDet));
1299 }
1300 return kTRUE;
1301 }
1302 if (strstr(name.Data(),patternrphi4)) {
1303 fCalibraMode->SetNrphi(i, 4);
1304 if (fDebugLevel > 1) {
1305 AliInfo(Form("fNbDet %d and 4",fNbDet));
1306 }
1307 return kTRUE;
1308 }
1309 if (strstr(name.Data(),patternrphi5)) {
1310 fCalibraMode->SetNrphi(i, 5);
1311 if (fDebugLevel > 1) {
1312 AliInfo(Form("fNbDet %d and 5",fNbDet));
1313 }
1314 return kTRUE;
1315 }
1316 if (strstr(name.Data(),patternrphi6)) {
1317 fCalibraMode->SetNrphi(i, 6);
1318 if (fDebugLevel > 1) {
1319 AliInfo(Form("fNbDet %d and 6",fNbDet));
1320 }
1321 return kTRUE;
1322 }
1323
1324 if (fDebugLevel > 1) {
1325 AliInfo(Form("fNbDet %d and rest",fNbDet));
1326 }
1327 fCalibraMode->SetNrphi(i ,0);
1328 return kFALSE;
1329
1330}
1331//_____________________________________________________________________________
1332Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1333{
1334 //
1335 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1336 // corresponding to the given TObject
1337 //
1338
1339 // Some patterns
1340 const Char_t *patternz0 = "Nz0";
1341 const Char_t *patternz1 = "Nz1";
1342 const Char_t *patternz2 = "Nz2";
1343 const Char_t *patternz3 = "Nz3";
1344 const Char_t *patternz4 = "Nz4";
1345
1346 const Char_t *patternrphi10 = "Nrphi10";
1347 const Char_t *patternrphi100 = "Nrphi100";
1348 const Char_t *patternz10 = "Nz10";
1349 const Char_t *patternz100 = "Nz100";
1350
1351 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1352 fCalibraMode->SetAllTogether(i);
1353 fNbDet = 540;
1354 if (fDebugLevel > 1) {
1355 AliInfo(Form("fNbDet %d and 100",fNbDet));
1356 }
1357 return kTRUE;
1358 }
1359 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1360 fCalibraMode->SetPerSuperModule(i);
1361 fNbDet = 30;
1362 if (fDebugLevel > 1) {
1363 AliInfo(Form("fNbDet %d and 10",fNbDet));
1364 }
1365 return kTRUE;
1366 }
1367 if (strstr(name.Data(),patternz0)) {
1368 fCalibraMode->SetNz(i, 0);
1369 if (fDebugLevel > 1) {
1370 AliInfo(Form("fNbDet %d and 0",fNbDet));
1371 }
1372 return kTRUE;
1373 }
1374 if (strstr(name.Data(),patternz1)) {
1375 fCalibraMode->SetNz(i ,1);
1376 if (fDebugLevel > 1) {
1377 AliInfo(Form("fNbDet %d and 1",fNbDet));
1378 }
1379 return kTRUE;
1380 }
1381 if (strstr(name.Data(),patternz2)) {
1382 fCalibraMode->SetNz(i ,2);
1383 if (fDebugLevel > 1) {
1384 AliInfo(Form("fNbDet %d and 2",fNbDet));
1385 }
1386 return kTRUE;
1387 }
1388 if (strstr(name.Data(),patternz3)) {
1389 fCalibraMode->SetNz(i ,3);
1390 if (fDebugLevel > 1) {
1391 AliInfo(Form("fNbDet %d and 3",fNbDet));
1392 }
1393 return kTRUE;
1394 }
1395 if (strstr(name.Data(),patternz4)) {
1396 fCalibraMode->SetNz(i ,4);
1397 if (fDebugLevel > 1) {
1398 AliInfo(Form("fNbDet %d and 4",fNbDet));
1399 }
1400 return kTRUE;
1401 }
1402
1403 if (fDebugLevel > 1) {
1404 AliInfo(Form("fNbDet %d and rest",fNbDet));
1405 }
1406 fCalibraMode->SetNz(i ,0);
1407 return kFALSE;
1408}
1409//______________________________________________________________________
1410void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1411 //
1412 // Remove the results too far from the mean value and rms
1413 // type: 0 gain, 1 vdrift
1414 // perdetector
1415 //
1416
1417 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1418 if(loop != 540) {
1419 AliInfo("The Vector Fit is not complete!");
1420 return;
1421 }
1422 Int_t detector = -1;
1423 Int_t sector = -1;
1424 Float_t value = 0.0;
1425
1426 /////////////////////////////////
1427 // Calculate the mean values
1428 ////////////////////////////////
1429 // Initialisation
1430 ////////////////////////
1431 Double_t meanAll = 0.0;
1432 Double_t rmsAll = 0.0;
1433 Int_t countAll = 0;
1434 ////////////
1435 // compute
1436 ////////////
1437 for (Int_t k = 0; k < loop; k++) {
1438 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1439 sector = GetSector(detector);
1440 if(perdetector){
1441 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1442 if(value > 0.0) {
1443 rmsAll += value*value;
1444 meanAll += value;
1445 countAll++;
1446 }
1447 }
1448 else {
1449 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1450 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1451 for (Int_t row = 0; row < rowMax; row++) {
1452 for (Int_t col = 0; col < colMax; col++) {
1453 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1454 if(value > 0.0) {
1455 rmsAll += value*value;
1456 meanAll += value;
1457 countAll++;
1458 }
1459
1460 } // Col
1461 } // Row
1462 }
1463 }
1464 if(countAll > 0) {
1465 meanAll = meanAll/countAll;
1466 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1467 }
1468 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1469 /////////////////////////////////////////////////
1470 // Remove outliers
1471 ////////////////////////////////////////////////
1472 Double_t defaultvalue = -1.0;
1473 if(type==1) defaultvalue = -1.5;
1474 for (Int_t k = 0; k < loop; k++) {
1475 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1476 sector = GetSector(detector);
1477 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1478 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1479 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1480
1481 // remove the results too far away
1482 for (Int_t row = 0; row < rowMax; row++) {
1483 for (Int_t col = 0; col < colMax; col++) {
1484 value = coef[(Int_t)(col*rowMax+row)];
1485 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1486 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1487 }
1488 } // Col
1489 } // Row
1490 }
1491}
1492//______________________________________________________________________
1493void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1494 //
1495 // Remove the results too far from the mean and rms
1496 // perdetector
1497 //
1498
1499 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1500 if(loop != 540) {
1501 AliInfo("The Vector Fit is not complete!");
1502 return;
1503 }
1504 Int_t detector = -1;
1505 Int_t sector = -1;
1506 Float_t value = 0.0;
1507
1508 /////////////////////////////////
1509 // Calculate the mean values
1510 ////////////////////////////////
1511 // Initialisation
1512 ////////////////////////
1513 Double_t meanAll = 0.0;
1514 Double_t rmsAll = 0.0;
1515 Int_t countAll = 0;
1516 /////////////
1517 // compute
1518 ////////////
1519 for (Int_t k = 0; k < loop; k++) {
1520 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1521 sector = GetSector(detector);
1522 if(perdetector){
1523 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1524 if(value < 70.0) {
1525 meanAll += value;
1526 rmsAll += value*value;
1527 countAll++;
1528 }
1529 }
1530 else {
1531 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1532 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1533 for (Int_t row = 0; row < rowMax; row++) {
1534 for (Int_t col = 0; col < colMax; col++) {
1535 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1536 if(value < 70.0) {
1537 rmsAll += value*value;
1538 meanAll += value;
1539 countAll++;
1540 }
1541 } // Col
1542 } // Row
1543 }
1544 }
1545 if(countAll > 0) {
1546 meanAll = meanAll/countAll;
1547 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1548 }
1549 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1550 /////////////////////////////////////////////////
1551 // Remove outliers
1552 ////////////////////////////////////////////////
1553 for (Int_t k = 0; k < loop; k++) {
1554 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1555 sector = GetSector(detector);
1556 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1557 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1558 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1559
1560 // remove the results too far away
1561 for (Int_t row = 0; row < rowMax; row++) {
1562 for (Int_t col = 0; col < colMax; col++) {
1563 value = coef[(Int_t)(col*rowMax+row)];
1564 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) coef[(Int_t)(col*rowMax+row)] = 100.0;
1565 } // Col
1566 } // Row
1567 }
1568}
1569//______________________________________________________________________
1570void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1571 //
1572 // ofwhat is equaled to 0: mean value of all passing detectors
1573 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1574 //
1575
1576 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1577 if(loop != 540) {
1578 AliInfo("The Vector Fit is not complete!");
1579 return;
1580 }
1581 Int_t detector = -1;
1582 Int_t sector = -1;
1583 Float_t value = 0.0;
1584
1585 /////////////////////////////////
1586 // Calculate the mean values
1587 ////////////////////////////////
1588 // Initialisation
1589 ////////////////////////
1590 Double_t meanAll = 0.0;
1591 Double_t meanSupermodule[18];
1592 Double_t meanDetector[540];
1593 Double_t rmsAll = 0.0;
1594 Double_t rmsSupermodule[18];
1595 Double_t rmsDetector[540];
1596 Int_t countAll = 0;
1597 Int_t countSupermodule[18];
1598 Int_t countDetector[540];
1599 for(Int_t sm = 0; sm < 18; sm++){
1600 rmsSupermodule[sm] = 0.0;
1601 meanSupermodule[sm] = 0.0;
1602 countSupermodule[sm] = 0;
1603 }
1604 for(Int_t det = 0; det < 540; det++){
1605 rmsDetector[det] = 0.0;
1606 meanDetector[det] = 0.0;
1607 countDetector[det] = 0;
1608 }
1609 ////////////
1610 // compute
1611 ////////////
1612 for (Int_t k = 0; k < loop; k++) {
1613 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1614 sector = GetSector(detector);
1615 if(perdetector){
1616 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1617 if(value > 0.0) {
1618 rmsDetector[detector] += value*value;
1619 meanDetector[detector] += value;
1620 countDetector[detector]++;
1621 rmsSupermodule[sector] += value*value;
1622 meanSupermodule[sector] += value;
1623 countSupermodule[sector]++;
1624 rmsAll += value*value;
1625 meanAll += value;
1626 countAll++;
1627 }
1628 }
1629 else {
1630 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1631 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1632 for (Int_t row = 0; row < rowMax; row++) {
1633 for (Int_t col = 0; col < colMax; col++) {
1634 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1635 if(value > 0.0) {
1636 rmsDetector[detector] += value*value;
1637 meanDetector[detector] += value;
1638 countDetector[detector]++;
1639 rmsSupermodule[sector] += value*value;
1640 meanSupermodule[sector] += value;
1641 countSupermodule[sector]++;
1642 rmsAll += value*value;
1643 meanAll += value;
1644 countAll++;
1645 }
1646
1647 } // Col
1648 } // Row
1649 }
1650 }
1651 if(countAll > 0) {
1652 meanAll = meanAll/countAll;
1653 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1654 }
1655 for(Int_t sm = 0; sm < 18; sm++){
1656 if(countSupermodule[sm] > 0) {
1657 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1658 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1659 }
1660 }
1661 for(Int_t det = 0; det < 540; det++){
1662 if(countDetector[det] > 0) {
1663 meanDetector[det] = meanDetector[det]/countDetector[det];
1664 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1665 }
1666 }
1667 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1668 ///////////////////////////////////////////////
1669 // Put the mean value for the no-fitted
1670 /////////////////////////////////////////////
1671 for (Int_t k = 0; k < loop; k++) {
1672 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1673 sector = GetSector(detector);
1674 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1675 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1676 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1677
1678 for (Int_t row = 0; row < rowMax; row++) {
1679 for (Int_t col = 0; col < colMax; col++) {
1680 value = coef[(Int_t)(col*rowMax+row)];
1681 if(value < 0.0) {
1682 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1683 if(ofwhat == 1){
1684 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1685 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1686 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1687 }
1688 }
1689 // Debug
1690 if(fDebugLevel > 1){
1691
1692 if ( !fDebugStreamer ) {
1693 //debug stream
1694 TDirectory *backup = gDirectory;
1695 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1696 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1697 }
1698
1699 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1700
1701 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1702 "detector="<<detector<<
1703 "sector="<<sector<<
1704 "row="<<row<<
1705 "col="<<col<<
1706 "before="<<value<<
1707 "after="<<coefnow<<
1708 "\n";
1709 }
1710 } // Col
1711 } // Row
1712 }
1713}
1714//______________________________________________________________________
1715void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1716 //
1717 // ofwhat is equaled to 0: mean value of all passing detectors
1718 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1719 //
1720
1721 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1722 if(loop != 540) {
1723 AliInfo("The Vector Fit is not complete!");
1724 return;
1725 }
1726 Int_t detector = -1;
1727 Int_t sector = -1;
1728 Float_t value = 0.0;
1729
1730 /////////////////////////////////
1731 // Calculate the mean values
1732 ////////////////////////////////
1733 // Initialisation
1734 ////////////////////////
1735 Double_t meanAll = 0.0;
1736 Double_t rmsAll = 0.0;
1737 Double_t meanSupermodule[18];
1738 Double_t rmsSupermodule[18];
1739 Double_t meanDetector[540];
1740 Double_t rmsDetector[540];
1741 Int_t countAll = 0;
1742 Int_t countSupermodule[18];
1743 Int_t countDetector[540];
1744 for(Int_t sm = 0; sm < 18; sm++){
1745 rmsSupermodule[sm] = 0.0;
1746 meanSupermodule[sm] = 0.0;
1747 countSupermodule[sm] = 0;
1748 }
1749 for(Int_t det = 0; det < 540; det++){
1750 rmsDetector[det] = 0.0;
1751 meanDetector[det] = 0.0;
1752 countDetector[det] = 0;
1753 }
1754 // compute
1755 ////////////
1756 for (Int_t k = 0; k < loop; k++) {
1757 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1758 sector = GetSector(detector);
1759 if(perdetector){
1760 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1761 if(value < 70.0) {
1762 rmsDetector[detector] += value*value;
1763 meanDetector[detector] += value;
1764 countDetector[detector]++;
1765 rmsSupermodule[sector] += value*value;
1766 meanSupermodule[sector] += value;
1767 countSupermodule[sector]++;
1768 meanAll += value;
1769 rmsAll += value*value;
1770 countAll++;
1771 }
1772 }
1773 else {
1774 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1775 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1776 for (Int_t row = 0; row < rowMax; row++) {
1777 for (Int_t col = 0; col < colMax; col++) {
1778 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1779 if(value < 70.0) {
1780 rmsDetector[detector] += value*value;
1781 meanDetector[detector] += value;
1782 countDetector[detector]++;
1783 rmsSupermodule[sector] += value*value;
1784 meanSupermodule[sector] += value;
1785 countSupermodule[sector]++;
1786 rmsAll += value*value;
1787 meanAll += value;
1788 countAll++;
1789 }
1790
1791 } // Col
1792 } // Row
1793 }
1794 }
1795 if(countAll > 0) {
1796 meanAll = meanAll/countAll;
1797 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1798 }
1799 for(Int_t sm = 0; sm < 18; sm++){
1800 if(countSupermodule[sm] > 0) {
1801 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1802 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1803 }
1804 }
1805 for(Int_t det = 0; det < 540; det++){
1806 if(countDetector[det] > 0) {
1807 meanDetector[det] = meanDetector[det]/countDetector[det];
1808 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1809 }
1810 }
1811 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1812 ////////////////////////////////////////////
1813 // Put the mean value for the no-fitted
1814 /////////////////////////////////////////////
1815 for (Int_t k = 0; k < loop; k++) {
1816 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1817 sector = GetSector(detector);
1818 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1819 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1820 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1821
1822 for (Int_t row = 0; row < rowMax; row++) {
1823 for (Int_t col = 0; col < colMax; col++) {
1824 value = coef[(Int_t)(col*rowMax+row)];
1825 if(value > 70.0) {
1826 if((ofwhat == 0) && (meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1827 if(ofwhat == 1){
1828 if((meanDetector[detector] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1829 else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1830 else if((meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1831 }
1832 }
1833 // Debug
1834 if(fDebugLevel > 1){
1835
1836 if ( !fDebugStreamer ) {
1837 //debug stream
1838 TDirectory *backup = gDirectory;
1839 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1840 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1841 }
1842
1843 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1844
1845 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1846 "detector="<<detector<<
1847 "sector="<<sector<<
1848 "row="<<row<<
1849 "col="<<col<<
1850 "before="<<value<<
1851 "after="<<coefnow<<
1852 "\n";
1853 }
1854 } // Col
1855 } // Row
1856 }
1857
1858}
1859//_____________________________________________________________________________
1860AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1861{
1862 //
1863 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1864 // It takes the mean value of the coefficients per detector
1865 // This object has to be written in the database
1866 //
1867
1868 // Create the DetObject
1869 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1870
1871 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1872 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1873 Int_t detector = -1;
1874 Float_t value = 0.0;
1875
1876 //
1877 for (Int_t k = 0; k < loop; k++) {
1878 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1879 Float_t mean = 0.0;
1880 if(perdetector){
1881 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1882 }
1883 else {
1884 Int_t count = 0;
1885 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1886 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1887 for (Int_t row = 0; row < rowMax; row++) {
1888 for (Int_t col = 0; col < colMax; col++) {
1889 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1890 mean += TMath::Abs(value);
1891 count++;
1892 } // Col
1893 } // Row
1894 if(count > 0) mean = mean/count;
1895 }
1896 object->SetValue(detector,mean);
1897 }
1898
1899 return object;
1900}
1901//_____________________________________________________________________________
1902AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1903{
1904 //
1905 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1906 // It takes the mean value of the coefficients per detector
1907 // This object has to be written in the database
1908 //
1909
1910 // Create the DetObject
1911 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1912
1913
1914 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1915 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1916 Int_t detector = -1;
1917 Float_t value = 0.0;
1918
1919 for (Int_t k = 0; k < loop; k++) {
1920 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1921 Float_t mean = 0.0;
1922 if(perdetector){
1923 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1924 if(!meanOtherBefore){
1925 if(value > 0) value = value*scaleFitFactor;
1926 }
1927 else value = value*scaleFitFactor;
1928 mean = TMath::Abs(value);
1929 }
1930 else{
1931 Int_t count = 0;
1932 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1933 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1934 for (Int_t row = 0; row < rowMax; row++) {
1935 for (Int_t col = 0; col < colMax; col++) {
1936 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1937 if(!meanOtherBefore) {
1938 if(value > 0) value = value*scaleFitFactor;
1939 }
1940 else value = value*scaleFitFactor;
1941 mean += TMath::Abs(value);
1942 count++;
1943 } // Col
1944 } // Row
1945 if(count > 0) mean = mean/count;
1946 }
1947 object->SetValue(detector,mean);
1948 }
1949
1950 return object;
1951}
1952//_____________________________________________________________________________
1953AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1954{
1955 //
1956 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1957 // It takes the min value of the coefficients per detector
1958 // This object has to be written in the database
1959 //
1960
1961 // Create the DetObject
1962 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1963
1964 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1965 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1966 Int_t detector = -1;
1967 Float_t value = 0.0;
1968
1969 for (Int_t k = 0; k < loop; k++) {
1970 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1971 Float_t min = 100.0;
1972 if(perdetector){
1973 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1974 // check successful
1975 if(value > 70.0) value = value-100.0;
1976 //
1977 min = value;
1978 }
1979 else{
1980 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1981 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1982 for (Int_t row = 0; row < rowMax; row++) {
1983 for (Int_t col = 0; col < colMax; col++) {
1984 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1985 // check successful
1986 if(value > 70.0) value = value-100.0;
1987 //
1988 if(min > value) min = value;
1989 } // Col
1990 } // Row
1991 }
1992 object->SetValue(detector,min);
1993 }
1994
1995 return object;
1996
1997}
1998//_____________________________________________________________________________
1999AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2000{
2001 //
2002 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2003 // It takes the min value of the coefficients per detector
2004 // This object has to be written in the database
2005 //
2006
2007 // Create the DetObject
2008 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2009
2010
2011 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2012 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2013 Int_t detector = -1;
2014 Float_t value = 0.0;
2015
2016 for (Int_t k = 0; k < loop; k++) {
2017 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2018 /*
2019 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2020 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2021 Float_t min = 100.0;
2022 for (Int_t row = 0; row < rowMax; row++) {
2023 for (Int_t col = 0; col < colMax; col++) {
2024 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2025 mean += -TMath::Abs(value);
2026 count++;
2027 } // Col
2028 } // Row
2029 if(count > 0) mean = mean/count;
2030 */
2031 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2032 object->SetValue(detector,-TMath::Abs(value));
2033 }
2034
2035 return object;
2036
2037}
2038//_____________________________________________________________________________
2039TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2040{
2041 //
2042 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2043 // You need first to create the object for the detectors,
2044 // where the mean value is put.
2045 // This object has to be written in the database
2046 //
2047
2048 // Create the DetObject
2049 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2050
2051 if(!vectorFit){
2052 for(Int_t k = 0; k < 540; k++){
2053 AliTRDCalROC *calROC = object->GetCalROC(k);
2054 Int_t nchannels = calROC->GetNchannels();
2055 for(Int_t ch = 0; ch < nchannels; ch++){
2056 calROC->SetValue(ch,1.0);
2057 }
2058 }
2059 }
2060 else{
2061
2062 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2063 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2064 Int_t detector = -1;
2065 Float_t value = 0.0;
2066
2067 for (Int_t k = 0; k < loop; k++) {
2068 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2069 AliTRDCalROC *calROC = object->GetCalROC(detector);
2070 Float_t mean = detobject->GetValue(detector);
2071 if(TMath::Abs(mean) <= 0.0000000001) continue;
2072 Int_t rowMax = calROC->GetNrows();
2073 Int_t colMax = calROC->GetNcols();
2074 for (Int_t row = 0; row < rowMax; row++) {
2075 for (Int_t col = 0; col < colMax; col++) {
2076 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2077 if(value > 0) value = value*scaleFitFactor;
2078 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2079 } // Col
2080 } // Row
2081 }
2082 }
2083
2084 return object;
2085}
2086//_____________________________________________________________________________
2087TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2088{
2089 //
2090 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2091 // You need first to create the object for the detectors,
2092 // where the mean value is put.
2093 // This object has to be written in the database
2094 //
2095
2096 // Create the DetObject
2097 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2098
2099 if(!vectorFit){
2100 for(Int_t k = 0; k < 540; k++){
2101 AliTRDCalROC *calROC = object->GetCalROC(k);
2102 Int_t nchannels = calROC->GetNchannels();
2103 for(Int_t ch = 0; ch < nchannels; ch++){
2104 calROC->SetValue(ch,1.0);
2105 }
2106 }
2107 }
2108 else {
2109
2110 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2111 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2112 Int_t detector = -1;
2113 Float_t value = 0.0;
2114
2115 for (Int_t k = 0; k < loop; k++) {
2116 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2117 AliTRDCalROC *calROC = object->GetCalROC(detector);
2118 Float_t mean = detobject->GetValue(detector);
2119 if(mean == 0) continue;
2120 Int_t rowMax = calROC->GetNrows();
2121 Int_t colMax = calROC->GetNcols();
2122 for (Int_t row = 0; row < rowMax; row++) {
2123 for (Int_t col = 0; col < colMax; col++) {
2124 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2125 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2126 } // Col
2127 } // Row
2128 }
2129 }
2130 return object;
2131
2132}
2133//_____________________________________________________________________________
2134TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2135{
2136 //
2137 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2138 // You need first to create the object for the detectors,
2139 // where the mean value is put.
2140 // This object has to be written in the database
2141 //
2142
2143 // Create the DetObject
2144 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2145
2146 if(!vectorFit){
2147 for(Int_t k = 0; k < 540; k++){
2148 AliTRDCalROC *calROC = object->GetCalROC(k);
2149 Int_t nchannels = calROC->GetNchannels();
2150 for(Int_t ch = 0; ch < nchannels; ch++){
2151 calROC->SetValue(ch,0.0);
2152 }
2153 }
2154 }
2155 else {
2156
2157 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2158 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2159 Int_t detector = -1;
2160 Float_t value = 0.0;
2161
2162 for (Int_t k = 0; k < loop; k++) {
2163 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2164 AliTRDCalROC *calROC = object->GetCalROC(detector);
2165 Float_t min = detobject->GetValue(detector);
2166 Int_t rowMax = calROC->GetNrows();
2167 Int_t colMax = calROC->GetNcols();
2168 for (Int_t row = 0; row < rowMax; row++) {
2169 for (Int_t col = 0; col < colMax; col++) {
2170 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2171 // check successful
2172 if(value > 70.0) value = value - 100.0;
2173 //
2174 calROC->SetValue(col,row,value-min);
2175 } // Col
2176 } // Row
2177 }
2178 }
2179 return object;
2180
2181}
2182//_____________________________________________________________________________
2183TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2184{
2185 //
2186 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2187 // This object has to be written in the database
2188 //
2189
2190 // Create the DetObject
2191 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2192
2193 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2194 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2195 Int_t detector = -1;
2196 Float_t value = 0.0;
2197
2198 for (Int_t k = 0; k < loop; k++) {
2199 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2200 AliTRDCalROC *calROC = object->GetCalROC(detector);
2201 Int_t rowMax = calROC->GetNrows();
2202 Int_t colMax = calROC->GetNcols();
2203 for (Int_t row = 0; row < rowMax; row++) {
2204 for (Int_t col = 0; col < colMax; col++) {
2205 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2206 calROC->SetValue(col,row,TMath::Abs(value));
2207 } // Col
2208 } // Row
2209 }
2210
2211 return object;
2212
2213}
2214//_____________________________________________________________________________
2215AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2216{
2217 //
2218 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2219 // 0 successful fit 1 not successful fit
2220 // mean is the mean value over the successful fit
2221 // do not use it for t0: no meaning
2222 //
2223
2224 // Create the CalObject
2225 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2226 mean = 0.0;
2227 Int_t count = 0;
2228
2229 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2230 if(loop != 540) {
2231 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2232 for(Int_t k = 0; k < 540; k++){
2233 object->SetValue(k,1.0);
2234 }
2235 }
2236 Int_t detector = -1;
2237 Float_t value = 0.0;
2238
2239 for (Int_t k = 0; k < loop; k++) {
2240 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2241 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2242 if(value <= 0) object->SetValue(detector,1.0);
2243 else {
2244 object->SetValue(detector,0.0);
2245 mean += value;
2246 count++;
2247 }
2248 }
2249 if(count > 0) mean /= count;
2250 return object;
2251}
2252//_____________________________________________________________________________
2253TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2254{
2255 //
2256 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2257 // 0 not successful fit 1 successful fit
2258 // mean mean value over the successful fit
2259 //
2260
2261 // Create the CalObject
2262 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2263 mean = 0.0;
2264 Int_t count = 0;
2265
2266 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2267 if(loop != 540) {
2268 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2269 for(Int_t k = 0; k < 540; k++){
2270 AliTRDCalROC *calROC = object->GetCalROC(k);
2271 Int_t nchannels = calROC->GetNchannels();
2272 for(Int_t ch = 0; ch < nchannels; ch++){
2273 calROC->SetValue(ch,1.0);
2274 }
2275 }
2276 }
2277 Int_t detector = -1;
2278 Float_t value = 0.0;
2279
2280 for (Int_t k = 0; k < loop; k++) {
2281 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2282 AliTRDCalROC *calROC = object->GetCalROC(detector);
2283 Int_t nchannels = calROC->GetNchannels();
2284 for (Int_t ch = 0; ch < nchannels; ch++) {
2285 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2286 if(value <= 0) calROC->SetValue(ch,1.0);
2287 else {
2288 calROC->SetValue(ch,0.0);
2289 mean += value;
2290 count++;
2291 }
2292 } // channels
2293 }
2294 if(count > 0) mean /= count;
2295 return object;
2296}
2297//_____________________________________________________________________________
2298void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2299{
2300 //
2301 // Set FitPH if 1 then each detector will be fitted
2302 //
2303
2304 if (periodeFitPH > 0) {
2305 fFitPHPeriode = periodeFitPH;
2306 }
2307 else {
2308 AliInfo("periodeFitPH must be higher than 0!");
2309 }
2310
2311}
2312//_____________________________________________________________________________
2313void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2314{
2315 //
2316 // The fit of the deposited charge distribution begins at
2317 // histo->Mean()/beginFitCharge
2318 // You can here set beginFitCharge
2319 //
2320
2321 if (beginFitCharge > 0) {
2322 fBeginFitCharge = beginFitCharge;
2323 }
2324 else {
2325 AliInfo("beginFitCharge must be strict positif!");
2326 }
2327
2328}
2329
2330//_____________________________________________________________________________
2331void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2332{
2333 //
2334 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2335 // You can here set t0Shift0
2336 //
2337
2338 if (t0Shift > 0) {
2339 fT0Shift0 = t0Shift;
2340 }
2341 else {
2342 AliInfo("t0Shift0 must be strict positif!");
2343 }
2344
2345}
2346
2347//_____________________________________________________________________________
2348void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2349{
2350 //
2351 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2352 // You can here set t0Shift1
2353 //
2354
2355 if (t0Shift > 0) {
2356 fT0Shift1 = t0Shift;
2357 }
2358 else {
2359 AliInfo("t0Shift must be strict positif!");
2360 }
2361
2362}
2363
2364//_____________________________________________________________________________
2365void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2366{
2367 //
2368 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2369 // You can here set rangeFitPRF
2370 //
2371
2372 if ((rangeFitPRF > 0) &&
2373 (rangeFitPRF <= 1.5)) {
2374 fRangeFitPRF = rangeFitPRF;
2375 }
2376 else {
2377 AliInfo("rangeFitPRF must be between 0 and 1.0");
2378 }
2379
2380}
2381
2382//_____________________________________________________________________________
2383void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2384{
2385 //
2386 // Minimum entries for fitting
2387 //
2388
2389 if (minEntries > 0) {
2390 fMinEntries = minEntries;
2391 }
2392 else {
2393 AliInfo("fMinEntries must be >= 0.");
2394 }
2395
2396}
2397
2398//_____________________________________________________________________________
2399void AliTRDCalibraFit::SetRebin(Short_t rebin)
2400{
2401 //
2402 // Rebin with rebin time less bins the Ch histo
2403 // You can set here rebin that should divide the number of bins of CH histo
2404 //
2405
2406 if (rebin > 0) {
2407 fRebin = rebin;
2408 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2409 }
2410 else {
2411 AliInfo("You have to choose a positiv value!");
2412 }
2413
2414}
2415//_____________________________________________________________________________
2416Bool_t AliTRDCalibraFit::FillVectorFit()
2417{
2418 //
2419 // For the Fit functions fill the vector Fit
2420 //
2421
2422 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2423
2424 Int_t ntotal = 1;
2425 if (GetStack(fCountDet) == 2) {
2426 ntotal = 1728;
2427 }
2428 else {
2429 ntotal = 2304;
2430 }
2431
2432 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2433 Float_t *coef = new Float_t[ntotal];
2434 for (Int_t i = 0; i < ntotal; i++) {
2435 coef[i] = fCurrentCoefDetector[i];
2436 }
2437
2438 Int_t detector = fCountDet;
2439 // Set
2440 fitInfo->SetCoef(coef);
2441 fitInfo->SetDetector(detector);
2442 fVectorFit.Add((TObject *) fitInfo);
2443
2444 return kTRUE;
2445
2446}
2447//_____________________________________________________________________________
2448Bool_t AliTRDCalibraFit::FillVectorFit2()
2449{
2450 //
2451 // For the Fit functions fill the vector Fit
2452 //
2453
2454 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2455
2456 Int_t ntotal = 1;
2457 if (GetStack(fCountDet) == 2) {
2458 ntotal = 1728;
2459 }
2460 else {
2461 ntotal = 2304;
2462 }
2463
2464 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2465 Float_t *coef = new Float_t[ntotal];
2466 for (Int_t i = 0; i < ntotal; i++) {
2467 coef[i] = fCurrentCoefDetector2[i];
2468 }
2469
2470 Int_t detector = fCountDet;
2471 // Set
2472 fitInfo->SetCoef(coef);
2473 fitInfo->SetDetector(detector);
2474 fVectorFit2.Add((TObject *) fitInfo);
2475
2476 return kTRUE;
2477
2478}
2479//____________Functions for initialising the AliTRDCalibraFit in the code_________
2480Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2481{
2482 //
2483 // Init the number of expected bins and fDect1[i] fDect2[i]
2484 //
2485
2486 gStyle->SetPalette(1);
2487 gStyle->SetOptStat(1111);
2488 gStyle->SetPadBorderMode(0);
2489 gStyle->SetCanvasColor(10);
2490 gStyle->SetPadLeftMargin(0.13);
2491 gStyle->SetPadRightMargin(0.01);
2492
2493 // Mode groups of pads: the total number of bins!
2494 CalculNumberOfBinsExpected(i);
2495
2496 // Quick verification that we have the good pad calibration mode!
2497 if (fNumberOfBinsExpected != nbins) {
2498 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2499 return kFALSE;
2500 }
2501
2502 // Security for fDebug 3 and 4
2503 if ((fDebugLevel >= 3) &&
2504 ((fDet[0] > 5) ||
2505 (fDet[1] > 4) ||
2506 (fDet[2] > 17))) {
2507 AliInfo("This detector doesn't exit!");
2508 return kFALSE;
2509 }
2510
2511 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2512 CalculDect1Dect2(i);
2513
2514
2515 return kTRUE;
2516}
2517//____________Functions for initialising the AliTRDCalibraFit in the code_________
2518Bool_t AliTRDCalibraFit::InitFitCH()
2519{
2520 //
2521 // Init the fVectorFitCH for normalisation
2522 // Init the histo for debugging
2523 //
2524
2525 gDirectory = gROOT;
2526
2527 fScaleFitFactor = 0.0;
2528 fCurrentCoefDetector = new Float_t[2304];
2529 for (Int_t k = 0; k < 2304; k++) {
2530 fCurrentCoefDetector[k] = 0.0;
2531 }
2532 fVectorFit.SetName("gainfactorscoefficients");
2533
2534 // fDebug == 0 nothing
2535 // fDebug == 1 and fFitVoir no histo
2536 if (fDebugLevel == 1) {
2537 if(!CheckFitVoir()) return kFALSE;
2538 }
2539 //Get the CalDet object
2540 if(fAccCDB){
2541 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2542 if (!cal) {
2543 AliInfo("Could not get calibDB");
2544 return kFALSE;
2545 }
2546 if(fCalDet) delete fCalDet;
2547 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2548 }
2549 else{
2550 Float_t devalue = 1.0;
2551 if(fCalDet) delete fCalDet;
2552 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2553 for(Int_t k = 0; k < 540; k++){
2554 fCalDet->SetValue(k,devalue);
2555 }
2556 }
2557 return kTRUE;
2558
2559}
2560//____________Functions for initialising the AliTRDCalibraFit in the code_________
2561Bool_t AliTRDCalibraFit::InitFitPH()
2562{
2563 //
2564 // Init the arrays of results
2565 // Init the histos for debugging
2566 //
2567
2568 gDirectory = gROOT;
2569 fVectorFit.SetName("driftvelocitycoefficients");
2570 fVectorFit2.SetName("t0coefficients");
2571
2572 fCurrentCoefDetector = new Float_t[2304];
2573 for (Int_t k = 0; k < 2304; k++) {
2574 fCurrentCoefDetector[k] = 0.0;
2575 }
2576
2577 fCurrentCoefDetector2 = new Float_t[2304];
2578 for (Int_t k = 0; k < 2304; k++) {
2579 fCurrentCoefDetector2[k] = 0.0;
2580 }
2581
2582 //fDebug == 0 nothing
2583 // fDebug == 1 and fFitVoir no histo
2584 if (fDebugLevel == 1) {
2585 if(!CheckFitVoir()) return kFALSE;
2586 }
2587 //Get the CalDet object
2588 if(fAccCDB){
2589 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2590 if (!cal) {
2591 AliInfo("Could not get calibDB");
2592 return kFALSE;
2593 }
2594 if(fCalDet) delete fCalDet;
2595 if(fCalDet2) delete fCalDet2;
2596 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2597 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2598 }
2599 else{
2600 Float_t devalue = 1.5;
2601 Float_t devalue2 = 0.0;
2602 if(fCalDet) delete fCalDet;
2603 if(fCalDet2) delete fCalDet2;
2604 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2605 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2606 for(Int_t k = 0; k < 540; k++){
2607 fCalDet->SetValue(k,devalue);
2608 fCalDet2->SetValue(k,devalue2);
2609 }
2610 }
2611 return kTRUE;
2612}
2613//____________Functions for initialising the AliTRDCalibraFit in the code_________
2614Bool_t AliTRDCalibraFit::InitFitPRF()
2615{
2616 //
2617 // Init the calibration mode (Nz, Nrphi), the histograms for
2618 // debugging the fit methods if fDebug > 0,
2619 //
2620
2621 gDirectory = gROOT;
2622 fVectorFit.SetName("prfwidthcoefficients");
2623
2624 fCurrentCoefDetector = new Float_t[2304];
2625 for (Int_t k = 0; k < 2304; k++) {
2626 fCurrentCoefDetector[k] = 0.0;
2627 }
2628
2629 // fDebug == 0 nothing
2630 // fDebug == 1 and fFitVoir no histo
2631 if (fDebugLevel == 1) {
2632 if(!CheckFitVoir()) return kFALSE;
2633 }
2634 return kTRUE;
2635}
2636//____________Functions for initialising the AliTRDCalibraFit in the code_________
2637Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2638{
2639 //
2640 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2641 //
2642
2643 gDirectory = gROOT;
2644
2645 fCurrentCoefDetector = new Float_t[2304];
2646 fCurrentCoefDetector2 = new Float_t[2304];
2647 for (Int_t k = 0; k < 2304; k++) {
2648 fCurrentCoefDetector[k] = 0.0;
2649 fCurrentCoefDetector2[k] = 0.0;
2650 }
2651
2652 //printf("test0\n");
2653
2654 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2655 if (!cal) {
2656 AliInfo("Could not get calibDB");
2657 return kFALSE;
2658 }
2659
2660 //Get the CalDet object
2661 if(fAccCDB){
2662 if(fCalDet) delete fCalDet;
2663 if(fCalDet2) delete fCalDet2;
2664 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2665 //printf("test1\n");
2666 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2667 //printf("test2\n");
2668 for(Int_t k = 0; k < 540; k++){
2669 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2670 }
2671 //printf("test3\n");
2672 }
2673 else{
2674 Float_t devalue = 1.5;
2675 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2676 if(fCalDet) delete fCalDet;
2677 if(fCalDet2) delete fCalDet2;
2678 //printf("test1\n");
2679 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2680 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2681 //printf("test2\n");
2682 for(Int_t k = 0; k < 540; k++){
2683 fCalDet->SetValue(k,devalue);
2684 fCalDet2->SetValue(k,devalue2);
2685 }
2686 //printf("test3\n");
2687 }
2688 return kTRUE;
2689}
2690
2691//____________Functions for initialising the AliTRDCalibraFit in the code_________
2692void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2693{
2694 //
2695 // Init the current detector where we are fCountDet and the
2696 // next fCount for the functions Fit...
2697 //
2698
2699 // Loop on the Xbins of ch!!
2700 fCountDet = -1; // Current detector
2701 fCount = 0; // To find the next detector
2702
2703 // If fDebug >= 3
2704 if (fDebugLevel >= 3) {
2705 // Set countdet to the detector
2706 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2707 // Set counter to write at the end of the detector
2708 fCount = fDect2;
2709 // Get the right calib objects
2710 SetCalROC(i);
2711 }
2712 if(fDebugLevel == 1) {
2713 fCountDet = 0;
2714 fCalibraMode->CalculXBins(fCountDet,i);
2715 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2716 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2717 fCountDet++;
2718 fCalibraMode->CalculXBins(fCountDet,i);
2719 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2720 }
2721 }
2722 else {
2723 fCountDet++;
2724 }
2725 fCount = fCalibraMode->GetXbins(i);
2726 fCountDet--;
2727 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2728 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2729 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2730 ,(Int_t) GetStack(fCountDet)
2731 ,(Int_t) GetSector(fCountDet),i);
2732 }
2733}
2734//_______________________________________________________________________________
2735void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2736{
2737 //
2738 // Calculate the number of bins expected (calibration groups)
2739 //
2740
2741 fNumberOfBinsExpected = 0;
2742 // All
2743 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2744 fNumberOfBinsExpected = 1;
2745 return;
2746 }
2747 // Per supermodule
2748 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2749 fNumberOfBinsExpected = 18;
2750 return;
2751 }
2752 // More
2753 fCalibraMode->ModePadCalibration(2,i);
2754 fCalibraMode->ModePadFragmentation(0,2,0,i);
2755 fCalibraMode->SetDetChamb2(i);
2756 if (fDebugLevel > 1) {
2757 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2758 }
2759 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2760 fCalibraMode->ModePadCalibration(0,i);
2761 fCalibraMode->ModePadFragmentation(0,0,0,i);
2762 fCalibraMode->SetDetChamb0(i);
2763 if (fDebugLevel > 1) {
2764 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2765 }
2766 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2767
2768}
2769//_______________________________________________________________________________
2770void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2771{
2772 //
2773 // Calculate the range of fits
2774 //
2775
2776 fDect1 = -1;
2777 fDect2 = -1;
2778 if (fDebugLevel == 1) {
2779 fDect1 = fFitVoir;
2780 fDect2 = fDect1 +1;
2781 }
2782 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2783 fDect1 = 0;
2784 fDect2 = fNumberOfBinsExpected;
2785 }
2786 if (fDebugLevel >= 3) {
2787 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2788 fCalibraMode->CalculXBins(fCountDet,i);
2789 fDect1 = fCalibraMode->GetXbins(i);
2790 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2791 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2792 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2793 ,(Int_t) GetStack(fCountDet)
2794 ,(Int_t) GetSector(fCountDet),i);
2795 // Set for the next detector
2796 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2797 }
2798}
2799//_______________________________________________________________________________
2800Bool_t AliTRDCalibraFit::CheckFitVoir()
2801{
2802 //
2803 // Check if fFitVoir is in the range
2804 //
2805
2806 if (fFitVoir < fNumberOfBinsExpected) {
2807 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2808 }
2809 else {
2810 AliInfo("fFitVoir is out of range of the histo!");
2811 return kFALSE;
2812 }
2813 return kTRUE;
2814}
2815//____________Functions for initialising the AliTRDCalibraFit in the code_________
2816void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2817{
2818 //
2819 // See if we are in a new detector and update the
2820 // variables fNfragZ and fNfragRphi if yes
2821 // Will never happen for only one detector (3 and 4)
2822 // Doesn't matter for 2
2823 //
2824 if (fCount == idect) {
2825 // On en est au detector (or first detector in the group)
2826 fCountDet += 1;
2827 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2828 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2829 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2830 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2831 ,(Int_t) GetStack(fCountDet)
2832 ,(Int_t) GetSector(fCountDet),i);
2833 // Set for the next detector
2834 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2835 // calib objects
2836 SetCalROC(i);
2837 }
2838}
2839//____________Functions for initialising the AliTRDCalibraFit in the code_________
2840void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2841{
2842 //
2843 // Reconstruct the min pad row, max pad row, min pad col and
2844 // max pad col of the calibration group for the Fit functions
2845 // idect is the calibration group inside the detector
2846 //
2847 if (fDebugLevel != 1) {
2848 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2849 }
2850 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2851 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2852}
2853//____________Functions for initialising the AliTRDCalibraFit in the code_________
2854Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2855{
2856 //
2857 // For the case where there are not enough entries in the histograms
2858 // of the calibration group, the value present in the choosen database
2859 // will be put. A negativ sign enables to know that a fit was not possible.
2860 //
2861
2862 if (fDebugLevel == 1) {
2863 AliInfo("The element has not enough statistic to be fitted");
2864 }
2865 else if (fNbDet > 0){
2866 Int_t firstdetector = fCountDet;
2867 Int_t lastdetector = fCountDet+fNbDet;
2868 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2869 ,idect,firstdetector,lastdetector));
2870 // loop over detectors
2871 for(Int_t det = firstdetector; det < lastdetector; det++){
2872
2873 //Set the calibration object again
2874 fCountDet = det;
2875 SetCalROC(0);
2876
2877 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2878 // Put them at 1
2879 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2880 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2881 ,(Int_t) GetStack(fCountDet)
2882 ,(Int_t) GetSector(fCountDet),0);
2883 // Reconstruct row min row max
2884 ReconstructFitRowMinRowMax(idect,0);
2885
2886 // Calcul the coef from the database choosen for the detector
2887 CalculChargeCoefMean(kFALSE);
2888
2889 //stack 2, not stack 2
2890 Int_t factor = 0;
2891 if(GetStack(fCountDet) == 2) factor = 12;
2892 else factor = 16;
2893
2894 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2895 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2896 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2897 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2898 }
2899 }
2900
2901 //Put default value negative
2902 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2903 fCurrentCoefE = 0.0;
2904
2905 // Fill the stuff
2906 FillVectorFit();
2907 // Debug
2908 if(fDebugLevel > 1){
2909
2910 if ( !fDebugStreamer ) {
2911 //debug stream
2912 TDirectory *backup = gDirectory;
2913 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2914 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2915 }
2916
2917 Int_t detector = fCountDet;
2918 Int_t caligroup = idect;
2919 Short_t rowmin = fCalibraMode->GetRowMin(0);
2920 Short_t rowmax = fCalibraMode->GetRowMax(0);
2921 Short_t colmin = fCalibraMode->GetColMin(0);
2922 Short_t colmax = fCalibraMode->GetColMax(0);
2923 Float_t gf = fCurrentCoef[0];
2924 Float_t gfs = fCurrentCoef[1];
2925 Float_t gfE = fCurrentCoefE;
2926
2927 (*fDebugStreamer) << "FillFillCH" <<
2928 "detector=" << detector <<
2929 "caligroup=" << caligroup <<
2930 "rowmin=" << rowmin <<
2931 "rowmax=" << rowmax <<
2932 "colmin=" << colmin <<
2933 "colmax=" << colmax <<
2934 "gf=" << gf <<
2935 "gfs=" << gfs <<
2936 "gfE=" << gfE <<
2937 "\n";
2938
2939 }
2940 // Reset
2941 for (Int_t k = 0; k < 2304; k++) {
2942 fCurrentCoefDetector[k] = 0.0;
2943 }
2944
2945 }// loop detector
2946 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2947 }
2948 else {
2949
2950 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2951 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2952
2953 // Calcul the coef from the database choosen
2954 CalculChargeCoefMean(kFALSE);
2955
2956 //stack 2, not stack 2
2957 Int_t factor = 0;
2958 if(GetStack(fCountDet) == 2) factor = 12;
2959 else factor = 16;
2960
2961 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2962 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2963 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2964 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2965 }
2966 }
2967
2968 //Put default value negative
2969 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2970 fCurrentCoefE = 0.0;
2971
2972 FillFillCH(idect);
2973 }
2974
2975 return kTRUE;
2976}
2977
2978
2979//____________Functions for initialising the AliTRDCalibraFit in the code_________
2980Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2981{
2982 //
2983 // For the case where there are not enough entries in the histograms
2984 // of the calibration group, the value present in the choosen database
2985 // will be put. A negativ sign enables to know that a fit was not possible.
2986 //
2987 if (fDebugLevel == 1) {
2988 AliInfo("The element has not enough statistic to be fitted");
2989 }
2990 else if (fNbDet > 0) {
2991
2992 Int_t firstdetector = fCountDet;
2993 Int_t lastdetector = fCountDet+fNbDet;
2994 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2995 ,idect,firstdetector,lastdetector));
2996 // loop over detectors
2997 for(Int_t det = firstdetector; det < lastdetector; det++){
2998
2999 //Set the calibration object again
3000 fCountDet = det;
3001 SetCalROC(1);
3002
3003 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3004 // Put them at 1
3005 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3006 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3007 ,(Int_t) GetStack(fCountDet)
3008 ,(Int_t) GetSector(fCountDet),1);
3009 // Reconstruct row min row max
3010 ReconstructFitRowMinRowMax(idect,1);
3011
3012 // Calcul the coef from the database choosen for the detector
3013 CalculVdriftCoefMean();
3014 CalculT0CoefMean();
3015
3016 //stack 2, not stack 2
3017 Int_t factor = 0;
3018 if(GetStack(fCountDet) == 2) factor = 12;
3019 else factor = 16;
3020
3021 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3022 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3023 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3024 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3025 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3026 }
3027 }
3028
3029 //Put default value negative
3030 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3031 fCurrentCoefE = 0.0;
3032 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3033 fCurrentCoefE2 = 0.0;
3034
3035 // Fill the stuff
3036 FillVectorFit();
3037 FillVectorFit2();
3038 // Debug
3039 if(fDebugLevel > 1){
3040
3041 if ( !fDebugStreamer ) {
3042 //debug stream
3043 TDirectory *backup = gDirectory;
3044 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3045 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3046 }
3047
3048
3049 Int_t detector = fCountDet;
3050 Int_t caligroup = idect;
3051 Short_t rowmin = fCalibraMode->GetRowMin(1);
3052 Short_t rowmax = fCalibraMode->GetRowMax(1);
3053 Short_t colmin = fCalibraMode->GetColMin(1);
3054 Short_t colmax = fCalibraMode->GetColMax(1);
3055 Float_t vf = fCurrentCoef[0];
3056 Float_t vs = fCurrentCoef[1];
3057 Float_t vfE = fCurrentCoefE;
3058 Float_t t0f = fCurrentCoef2[0];
3059 Float_t t0s = fCurrentCoef2[1];
3060 Float_t t0E = fCurrentCoefE2;
3061
3062
3063
3064 (* fDebugStreamer) << "FillFillPH"<<
3065 "detector="<<detector<<
3066 "nentries="<<nentries<<
3067 "caligroup="<<caligroup<<
3068 "rowmin="<<rowmin<<
3069 "rowmax="<<rowmax<<
3070 "colmin="<<colmin<<
3071 "colmax="<<colmax<<
3072 "vf="<<vf<<
3073 "vs="<<vs<<
3074 "vfE="<<vfE<<
3075 "t0f="<<t0f<<
3076 "t0s="<<t0s<<
3077 "t0E="<<t0E<<
3078 "\n";
3079 }
3080 // Reset
3081 for (Int_t k = 0; k < 2304; k++) {
3082 fCurrentCoefDetector[k] = 0.0;
3083 fCurrentCoefDetector2[k] = 0.0;
3084 }
3085
3086 }// loop detector
3087 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3088 }
3089 else {
3090
3091 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3092 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3093
3094 CalculVdriftCoefMean();
3095 CalculT0CoefMean();
3096
3097 //stack 2 and not stack 2
3098 Int_t factor = 0;
3099 if(GetStack(fCountDet) == 2) factor = 12;
3100 else factor = 16;
3101
3102
3103 // Fill the fCurrentCoefDetector 2
3104 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3105 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3106 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3107 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3108 }
3109 }
3110
3111 // Put the default value
3112 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3113 fCurrentCoefE = 0.0;
3114 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3115 fCurrentCoefE2 = 0.0;
3116
3117 FillFillPH(idect,nentries);
3118
3119 }
3120
3121 return kTRUE;
3122
3123}
3124
3125
3126//____________Functions for initialising the AliTRDCalibraFit in the code_________
3127Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3128{
3129 //
3130 // For the case where there are not enough entries in the histograms
3131 // of the calibration group, the value present in the choosen database
3132 // will be put. A negativ sign enables to know that a fit was not possible.
3133 //
3134
3135 if (fDebugLevel == 1) {
3136 AliInfo("The element has not enough statistic to be fitted");
3137 }
3138 else if (fNbDet > 0){
3139
3140 Int_t firstdetector = fCountDet;
3141 Int_t lastdetector = fCountDet+fNbDet;
3142 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3143 ,idect,firstdetector,lastdetector));
3144
3145 // loop over detectors
3146 for(Int_t det = firstdetector; det < lastdetector; det++){
3147
3148 //Set the calibration object again
3149 fCountDet = det;
3150 SetCalROC(2);
3151
3152 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3153 // Put them at 1
3154 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3155 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3156 ,(Int_t) GetStack(fCountDet)
3157 ,(Int_t) GetSector(fCountDet),2);
3158 // Reconstruct row min row max
3159 ReconstructFitRowMinRowMax(idect,2);
3160
3161 // Calcul the coef from the database choosen for the detector
3162 CalculPRFCoefMean();
3163
3164 //stack 2, not stack 2
3165 Int_t factor = 0;
3166 if(GetStack(fCountDet) == 2) factor = 12;
3167 else factor = 16;
3168
3169 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3170 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3171 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3172 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3173 }
3174 }
3175
3176 //Put default value negative
3177 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3178 fCurrentCoefE = 0.0;
3179
3180 // Fill the stuff
3181 FillVectorFit();
3182 // Debug
3183 if(fDebugLevel > 1){
3184
3185 if ( !fDebugStreamer ) {
3186 //debug stream
3187 TDirectory *backup = gDirectory;
3188 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3189 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3190 }
3191
3192 Int_t detector = fCountDet;
3193 Int_t layer = GetLayer(fCountDet);
3194 Int_t caligroup = idect;
3195 Short_t rowmin = fCalibraMode->GetRowMin(2);
3196 Short_t rowmax = fCalibraMode->GetRowMax(2);
3197 Short_t colmin = fCalibraMode->GetColMin(2);
3198 Short_t colmax = fCalibraMode->GetColMax(2);
3199 Float_t widf = fCurrentCoef[0];
3200 Float_t wids = fCurrentCoef[1];
3201 Float_t widfE = fCurrentCoefE;
3202
3203 (* fDebugStreamer) << "FillFillPRF"<<
3204 "detector="<<detector<<
3205 "layer="<<layer<<
3206 "caligroup="<<caligroup<<
3207 "rowmin="<<rowmin<<
3208 "rowmax="<<rowmax<<
3209 "colmin="<<colmin<<
3210 "colmax="<<colmax<<
3211 "widf="<<widf<<
3212 "wids="<<wids<<
3213 "widfE="<<widfE<<
3214 "\n";
3215 }
3216 // Reset
3217 for (Int_t k = 0; k < 2304; k++) {
3218 fCurrentCoefDetector[k] = 0.0;
3219 }
3220
3221 }// loop detector
3222 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3223 }
3224 else {
3225
3226 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3227 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3228
3229 CalculPRFCoefMean();
3230
3231 // stack 2 and not stack 2
3232 Int_t factor = 0;
3233 if(GetStack(fCountDet) == 2) factor = 12;
3234 else factor = 16;
3235
3236
3237 // Fill the fCurrentCoefDetector
3238 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3239 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3240 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3241 }
3242 }
3243
3244 // Put the default value
3245 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3246 fCurrentCoefE = 0.0;
3247
3248 FillFillPRF(idect);
3249 }
3250
3251 return kTRUE;
3252
3253}
3254//____________Functions for initialising the AliTRDCalibraFit in the code_________
3255Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3256{
3257 //
3258 // For the case where there are not enough entries in the histograms
3259 // of the calibration group, the value present in the choosen database
3260 // will be put. A negativ sign enables to know that a fit was not possible.
3261 //
3262
3263 // Calcul the coef from the database choosen
3264 CalculVdriftLorentzCoef();
3265
3266 Int_t factor = 0;
3267 if(GetStack(fCountDet) == 2) factor = 1728;
3268 else factor = 2304;
3269
3270
3271 // Fill the fCurrentCoefDetector
3272 for (Int_t k = 0; k < factor; k++) {
3273 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3274 // should be negative
3275 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3276 }
3277
3278
3279 //Put default opposite sign
3280 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3281 fCurrentCoefE = 0.0;
3282 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3283 fCurrentCoefE2 = 0.0;
3284
3285 FillFillLinearFitter();
3286
3287 return kTRUE;
3288}
3289
3290//____________Functions for initialising the AliTRDCalibraFit in the code_________
3291Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3292{
3293 //
3294 // Fill the coefficients found with the fits or other
3295 // methods from the Fit functions
3296 //
3297
3298 if (fDebugLevel != 1) {
3299 if (fNbDet > 0){
3300 Int_t firstdetector = fCountDet;
3301 Int_t lastdetector = fCountDet+fNbDet;
3302 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3303 ,idect,firstdetector,lastdetector));
3304 // loop over detectors
3305 for(Int_t det = firstdetector; det < lastdetector; det++){
3306
3307 //Set the calibration object again
3308 fCountDet = det;
3309 SetCalROC(0);
3310
3311 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3312 // Put them at 1
3313 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3314 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3315 ,(Int_t) GetStack(fCountDet)
3316 ,(Int_t) GetSector(fCountDet),0);
3317 // Reconstruct row min row max
3318 ReconstructFitRowMinRowMax(idect,0);
3319
3320 // Calcul the coef from the database choosen for the detector
3321 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3322 else CalculChargeCoefMean(kTRUE);
3323
3324 //stack 2, not stack 2
3325 Int_t factor = 0;
3326 if(GetStack(fCountDet) == 2) factor = 12;
3327 else factor = 16;
3328
3329 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3330 Double_t coeftoput = 1.0;
3331 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3332 else coeftoput = fCurrentCoef[0];
3333 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3334 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3335 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3336 }
3337 }
3338
3339 // Fill the stuff
3340 FillVectorFit();
3341 // Debug
3342 if(fDebugLevel > 1){
3343
3344 if ( !fDebugStreamer ) {
3345 //debug stream
3346 TDirectory *backup = gDirectory;
3347 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3348 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3349 }
3350
3351 Int_t detector = fCountDet;
3352 Int_t caligroup = idect;
3353 Short_t rowmin = fCalibraMode->GetRowMin(0);
3354 Short_t rowmax = fCalibraMode->GetRowMax(0);
3355 Short_t colmin = fCalibraMode->GetColMin(0);
3356 Short_t colmax = fCalibraMode->GetColMax(0);
3357 Float_t gf = fCurrentCoef[0];
3358 Float_t gfs = fCurrentCoef[1];
3359 Float_t gfE = fCurrentCoefE;
3360
3361 (*fDebugStreamer) << "FillFillCH" <<
3362 "detector=" << detector <<
3363 "caligroup=" << caligroup <<
3364 "rowmin=" << rowmin <<
3365 "rowmax=" << rowmax <<
3366 "colmin=" << colmin <<
3367 "colmax=" << colmax <<
3368 "gf=" << gf <<
3369 "gfs=" << gfs <<
3370 "gfE=" << gfE <<
3371 "\n";
3372
3373 }
3374 // Reset
3375 for (Int_t k = 0; k < 2304; k++) {
3376 fCurrentCoefDetector[k] = 0.0;
3377 }
3378
3379 }// loop detector
3380 //printf("Check the count now: fCountDet %d\n",fCountDet);
3381 }
3382 else{
3383
3384 Int_t factor = 0;
3385 if(GetStack(fCountDet) == 2) factor = 12;
3386 else factor = 16;
3387
3388 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3389 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3390 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3391 }
3392 }
3393
3394 FillFillCH(idect);
3395 }
3396 }
3397
3398 return kTRUE;
3399
3400}
3401//____________Functions for initialising the AliTRDCalibraFit in the code_________
3402Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3403{
3404 //
3405 // Fill the coefficients found with the fits or other
3406 // methods from the Fit functions
3407 //
3408
3409 if (fDebugLevel != 1) {
3410 if (fNbDet > 0){
3411
3412 Int_t firstdetector = fCountDet;
3413 Int_t lastdetector = fCountDet+fNbDet;
3414 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3415 ,idect,firstdetector,lastdetector));
3416
3417 // loop over detectors
3418 for(Int_t det = firstdetector; det < lastdetector; det++){
3419
3420 //Set the calibration object again
3421 fCountDet = det;
3422 SetCalROC(1);
3423
3424 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3425 // Put them at 1
3426 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3427 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3428 ,(Int_t) GetStack(fCountDet)
3429 ,(Int_t) GetSector(fCountDet),1);
3430 // Reconstruct row min row max
3431 ReconstructFitRowMinRowMax(idect,1);
3432
3433 // Calcul the coef from the database choosen for the detector
3434 CalculVdriftCoefMean();
3435 CalculT0CoefMean();
3436
3437 //stack 2, not stack 2
3438 Int_t factor = 0;
3439 if(GetStack(fCountDet) == 2) factor = 12;
3440 else factor = 16;
3441
3442 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3443 Double_t coeftoput = 1.5;
3444 Double_t coeftoput2 = 0.0;
3445
3446 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3447 else coeftoput = fCurrentCoef[0];
3448
3449 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3450 else coeftoput2 = fCurrentCoef2[0];
3451
3452 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3453 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3454 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3455 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3456 }
3457 }
3458
3459 // Fill the stuff
3460 FillVectorFit();
3461 FillVectorFit2();
3462 // Debug
3463 if(fDebugLevel > 1){
3464
3465 if ( !fDebugStreamer ) {
3466 //debug stream
3467 TDirectory *backup = gDirectory;
3468 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3469 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3470 }
3471
3472
3473 Int_t detector = fCountDet;
3474 Int_t caligroup = idect;
3475 Short_t rowmin = fCalibraMode->GetRowMin(1);
3476 Short_t rowmax = fCalibraMode->GetRowMax(1);
3477 Short_t colmin = fCalibraMode->GetColMin(1);
3478 Short_t colmax = fCalibraMode->GetColMax(1);
3479 Float_t vf = fCurrentCoef[0];
3480 Float_t vs = fCurrentCoef[1];
3481 Float_t vfE = fCurrentCoefE;
3482 Float_t t0f = fCurrentCoef2[0];
3483 Float_t t0s = fCurrentCoef2[1];
3484 Float_t t0E = fCurrentCoefE2;
3485
3486
3487
3488 (* fDebugStreamer) << "FillFillPH"<<
3489 "detector="<<detector<<
3490 "nentries="<<nentries<<
3491 "caligroup="<<caligroup<<
3492 "rowmin="<<rowmin<<
3493 "rowmax="<<rowmax<<
3494 "colmin="<<colmin<<
3495 "colmax="<<colmax<<
3496 "vf="<<vf<<
3497 "vs="<<vs<<
3498 "vfE="<<vfE<<
3499 "t0f="<<t0f<<
3500 "t0s="<<t0s<<
3501 "t0E="<<t0E<<
3502 "\n";
3503 }
3504 // Reset
3505 for (Int_t k = 0; k < 2304; k++) {
3506 fCurrentCoefDetector[k] = 0.0;
3507 fCurrentCoefDetector2[k] = 0.0;
3508 }
3509
3510 }// loop detector
3511 //printf("Check the count now: fCountDet %d\n",fCountDet);
3512 }
3513 else {
3514
3515 Int_t factor = 0;
3516 if(GetStack(fCountDet) == 2) factor = 12;
3517 else factor = 16;
3518
3519 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3520 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3521 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3522 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3523 }
3524 }
3525
3526 FillFillPH(idect,nentries);
3527 }
3528 }
3529 return kTRUE;
3530}
3531//____________Functions for initialising the AliTRDCalibraFit in the code_________
3532Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3533{
3534 //
3535 // Fill the coefficients found with the fits or other
3536 // methods from the Fit functions
3537 //
3538
3539 if (fDebugLevel != 1) {
3540 if (fNbDet > 0){
3541
3542 Int_t firstdetector = fCountDet;
3543 Int_t lastdetector = fCountDet+fNbDet;
3544 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3545 ,idect,firstdetector,lastdetector));
3546
3547 // loop over detectors
3548 for(Int_t det = firstdetector; det < lastdetector; det++){
3549
3550 //Set the calibration object again
3551 fCountDet = det;
3552 SetCalROC(2);
3553
3554 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3555 // Put them at 1
3556 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3557 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3558 ,(Int_t) GetStack(fCountDet)
3559 ,(Int_t) GetSector(fCountDet),2);
3560 // Reconstruct row min row max
3561 ReconstructFitRowMinRowMax(idect,2);
3562
3563 // Calcul the coef from the database choosen for the detector
3564 CalculPRFCoefMean();
3565
3566 //stack 2, not stack 2
3567 Int_t factor = 0;
3568 if(GetStack(fCountDet) == 2) factor = 12;
3569 else factor = 16;
3570
3571 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3572 Double_t coeftoput = 1.0;
3573 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3574 else coeftoput = fCurrentCoef[0];
3575 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3576 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3577 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3578 }
3579 }
3580
3581 // Fill the stuff
3582 FillVectorFit();
3583 // Debug
3584 if(fDebugLevel > 1){
3585
3586 if ( !fDebugStreamer ) {
3587 //debug stream
3588 TDirectory *backup = gDirectory;
3589 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3590 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3591 }
3592
3593 Int_t detector = fCountDet;
3594 Int_t layer = GetLayer(fCountDet);
3595 Int_t caligroup = idect;
3596 Short_t rowmin = fCalibraMode->GetRowMin(2);
3597 Short_t rowmax = fCalibraMode->GetRowMax(2);
3598 Short_t colmin = fCalibraMode->GetColMin(2);
3599 Short_t colmax = fCalibraMode->GetColMax(2);
3600 Float_t widf = fCurrentCoef[0];
3601 Float_t wids = fCurrentCoef[1];
3602 Float_t widfE = fCurrentCoefE;
3603
3604 (* fDebugStreamer) << "FillFillPRF"<<
3605 "detector="<<detector<<
3606 "layer="<<layer<<
3607 "caligroup="<<caligroup<<
3608 "rowmin="<<rowmin<<
3609 "rowmax="<<rowmax<<
3610 "colmin="<<colmin<<
3611 "colmax="<<colmax<<
3612 "widf="<<widf<<
3613 "wids="<<wids<<
3614 "widfE="<<widfE<<
3615 "\n";
3616 }
3617 // Reset
3618 for (Int_t k = 0; k < 2304; k++) {
3619 fCurrentCoefDetector[k] = 0.0;
3620 }
3621
3622 }// loop detector
3623 //printf("Check the count now: fCountDet %d\n",fCountDet);
3624 }
3625 else {
3626
3627 Int_t factor = 0;
3628 if(GetStack(fCountDet) == 2) factor = 12;
3629 else factor = 16;
3630
3631 // Pointer to the branch
3632 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3633 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3634 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3635 }
3636 }
3637 FillFillPRF(idect);
3638 }
3639 }
3640
3641 return kTRUE;
3642
3643}
3644//____________Functions for initialising the AliTRDCalibraFit in the code_________
3645Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3646{
3647 //
3648 // Fill the coefficients found with the fits or other
3649 // methods from the Fit functions
3650 //
3651
3652 Int_t factor = 0;
3653 if(GetStack(fCountDet) == 2) factor = 1728;
3654 else factor = 2304;
3655
3656 // Pointer to the branch
3657 for (Int_t k = 0; k < factor; k++) {
3658 fCurrentCoefDetector[k] = fCurrentCoef[0];
3659 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3660 }
3661
3662 FillFillLinearFitter();
3663
3664 return kTRUE;
3665
3666}
3667//________________________________________________________________________________
3668void AliTRDCalibraFit::FillFillCH(Int_t idect)
3669{
3670 //
3671 // DebugStream and fVectorFit
3672 //
3673
3674 // End of one detector
3675 if ((idect == (fCount-1))) {
3676 FillVectorFit();
3677 // Reset
3678 for (Int_t k = 0; k < 2304; k++) {
3679 fCurrentCoefDetector[k] = 0.0;
3680 }
3681 }
3682
3683 if(fDebugLevel > 1){
3684
3685 if ( !fDebugStreamer ) {
3686 //debug stream
3687 TDirectory *backup = gDirectory;
3688 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3689 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3690 }
3691
3692 Int_t detector = fCountDet;
3693 Int_t caligroup = idect;
3694 Short_t rowmin = fCalibraMode->GetRowMin(0);
3695 Short_t rowmax = fCalibraMode->GetRowMax(0);
3696 Short_t colmin = fCalibraMode->GetColMin(0);
3697 Short_t colmax = fCalibraMode->GetColMax(0);
3698 Float_t gf = fCurrentCoef[0];
3699 Float_t gfs = fCurrentCoef[1];
3700 Float_t gfE = fCurrentCoefE;
3701
3702 (*fDebugStreamer) << "FillFillCH" <<
3703 "detector=" << detector <<
3704 "caligroup=" << caligroup <<
3705 "rowmin=" << rowmin <<
3706 "rowmax=" << rowmax <<
3707 "colmin=" << colmin <<
3708 "colmax=" << colmax <<
3709 "gf=" << gf <<
3710 "gfs=" << gfs <<
3711 "gfE=" << gfE <<
3712 "\n";
3713
3714 }
3715}
3716//________________________________________________________________________________
3717void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3718{
3719 //
3720 // DebugStream and fVectorFit and fVectorFit2
3721 //
3722
3723 // End of one detector
3724 if ((idect == (fCount-1))) {
3725 FillVectorFit();
3726 FillVectorFit2();
3727 // Reset
3728 for (Int_t k = 0; k < 2304; k++) {
3729 fCurrentCoefDetector[k] = 0.0;
3730 fCurrentCoefDetector2[k] = 0.0;
3731 }
3732 }
3733
3734 if(fDebugLevel > 1){
3735
3736 if ( !fDebugStreamer ) {
3737 //debug stream
3738 TDirectory *backup = gDirectory;
3739 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3740 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3741 }
3742
3743
3744 Int_t detector = fCountDet;
3745 Int_t caligroup = idect;
3746 Short_t rowmin = fCalibraMode->GetRowMin(1);
3747 Short_t rowmax = fCalibraMode->GetRowMax(1);
3748 Short_t colmin = fCalibraMode->GetColMin(1);
3749 Short_t colmax = fCalibraMode->GetColMax(1);
3750 Float_t vf = fCurrentCoef[0];
3751 Float_t vs = fCurrentCoef[1];
3752 Float_t vfE = fCurrentCoefE;
3753 Float_t t0f = fCurrentCoef2[0];
3754 Float_t t0s = fCurrentCoef2[1];
3755 Float_t t0E = fCurrentCoefE2;
3756
3757
3758
3759 (* fDebugStreamer) << "FillFillPH"<<
3760 "detector="<<detector<<
3761 "nentries="<<nentries<<
3762 "caligroup="<<caligroup<<
3763 "rowmin="<<rowmin<<
3764 "rowmax="<<rowmax<<
3765 "colmin="<<colmin<<
3766 "colmax="<<colmax<<
3767 "vf="<<vf<<
3768 "vs="<<vs<<
3769 "vfE="<<vfE<<
3770 "t0f="<<t0f<<
3771 "t0s="<<t0s<<
3772 "t0E="<<t0E<<
3773 "\n";
3774 }
3775
3776}
3777//________________________________________________________________________________
3778void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3779{
3780 //
3781 // DebugStream and fVectorFit
3782 //
3783
3784 // End of one detector
3785 if ((idect == (fCount-1))) {
3786 FillVectorFit();
3787 // Reset
3788 for (Int_t k = 0; k < 2304; k++) {
3789 fCurrentCoefDetector[k] = 0.0;
3790 }
3791 }
3792
3793
3794 if(fDebugLevel > 1){
3795
3796 if ( !fDebugStreamer ) {
3797 //debug stream
3798 TDirectory *backup = gDirectory;
3799 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3800 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3801 }
3802
3803 Int_t detector = fCountDet;
3804 Int_t layer = GetLayer(fCountDet);
3805 Int_t caligroup = idect;
3806 Short_t rowmin = fCalibraMode->GetRowMin(2);
3807 Short_t rowmax = fCalibraMode->GetRowMax(2);
3808 Short_t colmin = fCalibraMode->GetColMin(2);
3809 Short_t colmax = fCalibraMode->GetColMax(2);
3810 Float_t widf = fCurrentCoef[0];
3811 Float_t wids = fCurrentCoef[1];
3812 Float_t widfE = fCurrentCoefE;
3813
3814 (* fDebugStreamer) << "FillFillPRF"<<
3815 "detector="<<detector<<
3816 "layer="<<layer<<
3817 "caligroup="<<caligroup<<
3818 "rowmin="<<rowmin<<
3819 "rowmax="<<rowmax<<
3820 "colmin="<<colmin<<
3821 "colmax="<<colmax<<
3822 "widf="<<widf<<
3823 "wids="<<wids<<
3824 "widfE="<<widfE<<
3825 "\n";
3826 }
3827
3828}
3829//________________________________________________________________________________
3830void AliTRDCalibraFit::FillFillLinearFitter()
3831{
3832 //
3833 // DebugStream and fVectorFit
3834 //
3835
3836 // End of one detector
3837 FillVectorFit();
3838 FillVectorFit2();
3839
3840
3841 // Reset
3842 for (Int_t k = 0; k < 2304; k++) {
3843 fCurrentCoefDetector[k] = 0.0;
3844 fCurrentCoefDetector2[k] = 0.0;
3845 }
3846
3847
3848 if(fDebugLevel > 1){
3849
3850 if ( !fDebugStreamer ) {
3851 //debug stream
3852 TDirectory *backup = gDirectory;
3853 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3854 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3855 }
3856
3857 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3858 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3859 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3860 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3861 Float_t tiltangle = padplane->GetTiltingAngle();
3862 Int_t detector = fCountDet;
3863 Int_t stack = GetStack(fCountDet);
3864 Int_t layer = GetLayer(fCountDet);
3865 Float_t vf = fCurrentCoef[0];
3866 Float_t vs = fCurrentCoef[1];
3867 Float_t vfE = fCurrentCoefE;
3868 Float_t lorentzangler = fCurrentCoef2[0];
3869 Float_t elorentzangler = fCurrentCoefE2;
3870 Float_t lorentzangles = fCurrentCoef2[1];
3871
3872 (* fDebugStreamer) << "FillFillLinearFitter"<<
3873 "detector="<<detector<<
3874 "stack="<<stack<<
3875 "layer="<<layer<<
3876 "rowmd="<<rowmd<<
3877 "r="<<r<<
3878 "tiltangle="<<tiltangle<<
3879 "vf="<<vf<<
3880 "vs="<<vs<<
3881 "vfE="<<vfE<<
3882 "lorentzangler="<<lorentzangler<<
3883 "Elorentzangler="<<elorentzangler<<
3884 "lorentzangles="<<lorentzangles<<
3885 "\n";
3886 }
3887
3888}
3889//
3890//____________Calcul Coef Mean_________________________________________________
3891//
3892//_____________________________________________________________________________
3893Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3894{
3895 //
3896 // For the detector Dect calcul the mean time 0
3897 // for the calibration group idect from the choosen database
3898 //
3899
3900 fCurrentCoef2[1] = 0.0;
3901 if(fDebugLevel != 1){
3902 if(((fCalibraMode->GetNz(1) > 0) ||
3903 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3904
3905 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3906 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3907 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3908 }
3909 }
3910
3911 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3912
3913 }
3914 else {
3915
3916 if(!fAccCDB){
3917 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3918 }
3919 else{
3920
3921 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3922 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3923 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3924 }
3925 }
3926 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3927
3928 }
3929 }
3930 }
3931 return kTRUE;
3932}
3933
3934//_____________________________________________________________________________
3935Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3936{
3937 //
3938 // For the detector Dect calcul the mean gain factor
3939 // for the calibration group idect from the choosen database
3940 //
3941
3942 fCurrentCoef[1] = 0.0;
3943 if(fDebugLevel != 1){
3944 if (((fCalibraMode->GetNz(0) > 0) ||
3945 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3946 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3947 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3948 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3949 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3950 }
3951 }
3952 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3953 }
3954 else {
3955 //Per detectors
3956 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3957 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3958 }
3959 }
3960 return kTRUE;
3961}
3962//_____________________________________________________________________________
3963Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3964{
3965 //
3966 // For the detector Dect calcul the mean sigma of pad response
3967 // function for the calibration group idect from the choosen database
3968 //
3969
3970 fCurrentCoef[1] = 0.0;
3971 if(fDebugLevel != 1){
3972 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3973 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3974 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3975 }
3976 }
3977 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3978 }
3979 return kTRUE;
3980}
3981//_____________________________________________________________________________
3982Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3983{
3984 //
3985 // For the detector dect calcul the mean drift velocity for the
3986 // calibration group idect from the choosen database
3987 //
3988
3989 fCurrentCoef[1] = 0.0;
3990 if(fDebugLevel != 1){
3991 if (((fCalibraMode->GetNz(1) > 0) ||
3992 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3993
3994 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3995 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3996 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3997 }
3998 }
3999
4000 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4001
4002 }
4003 else {
4004 //per detectors
4005 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4006 }
4007 }
4008 return kTRUE;
4009}
4010//_____________________________________________________________________________
4011Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4012{
4013 //
4014 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4015 //
4016
4017 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
4018 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4019
4020 return kTRUE;
4021}
4022//_____________________________________________________________________________
4023Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4024{
4025 //
4026 // Default width of the PRF if there is no database as reference
4027 //
4028 switch(layer)
4029 {
4030 // default database
4031 //case 0: return 0.515;
4032 //case 1: return 0.502;
4033 //case 2: return 0.491;
4034 //case 3: return 0.481;
4035 //case 4: return 0.471;
4036 //case 5: return 0.463;
4037 //default: return 0.0;
4038
4039 // fit database
4040 case 0: return 0.538429;
4041 case 1: return 0.524302;
4042 case 2: return 0.511591;
4043 case 3: return 0.500140;
4044 case 4: return 0.489821;
4045 case 5: return 0.480524;
4046 default: return 0.0;
4047 }
4048}
4049//________________________________________________________________________________
4050void AliTRDCalibraFit::SetCalROC(Int_t i)
4051{
4052 //
4053 // Set the calib object for fCountDet
4054 //
4055
4056 Float_t value = 0.0;
4057
4058 //Get the CalDet object
4059 if(fAccCDB){
4060 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4061 if (!cal) {
4062 AliInfo("Could not get calibDB");
4063 return;
4064 }
4065 switch (i)
4066 {
4067 case 0:
4068 if( fCalROC ){
4069 fCalROC->~AliTRDCalROC();
4070 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4071 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4072 break;
4073 case 1:
4074 if( fCalROC ){
4075 fCalROC->~AliTRDCalROC();
4076 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4077 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4078 if( fCalROC2 ){
4079 fCalROC2->~AliTRDCalROC();
4080 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4081 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4082 break;
4083 case 2:
4084 if( fCalROC ){
4085 fCalROC->~AliTRDCalROC();
4086 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4087 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4088 break;
4089 default: return;
4090 }
4091 }
4092 else{
4093 switch (i)
4094 {
4095 case 0:
4096 if(fCalROC) delete fCalROC;
4097 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4098 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4099 fCalROC->SetValue(k,1.0);
4100 }
4101 break;
4102 case 1:
4103 if(fCalROC) delete fCalROC;
4104 if(fCalROC2) delete fCalROC2;
4105 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4106 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4107 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4108 fCalROC->SetValue(k,1.0);
4109 fCalROC2->SetValue(k,0.0);
4110 }
4111 break;
4112 case 2:
4113 if(fCalROC) delete fCalROC;
4114 value = GetPRFDefault(GetLayer(fCountDet));
4115 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4116 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4117 fCalROC->SetValue(k,value);
4118 }
4119 break;
4120 default: return;
4121 }
4122 }
4123
4124}
4125//____________Fit Methods______________________________________________________
4126
4127//_____________________________________________________________________________
4128void AliTRDCalibraFit::FitPente(TH1* projPH)
4129{
4130 //
4131 // Slope methode for the drift velocity
4132 //
4133
4134 // Constants
4135 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4136 Int_t binmax = 0;
4137 Int_t binmin = 0;
4138 fPhd[0] = 0.0;
4139 fPhd[1] = 0.0;
4140 fPhd[2] = 0.0;
4141 Int_t ju = 0;
4142 fCurrentCoefE = 0.0;
4143 fCurrentCoefE2 = 0.0;
4144 fCurrentCoef[0] = 0.0;
4145 fCurrentCoef2[0] = 0.0;
4146 TLine *line = new TLine();
4147
4148 // Some variables
4149 TAxis *xpph = projPH->GetXaxis();
4150 Int_t nbins = xpph->GetNbins();
4151 Double_t lowedge = xpph->GetBinLowEdge(1);
4152 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4153 Double_t widbins = (upedge - lowedge) / nbins;
4154 Double_t limit = upedge + 0.5 * widbins;
4155 Bool_t put = kTRUE;
4156
4157 // Beginning of the signal
4158 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4159 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4160 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4161 }
4162 binmax = (Int_t) pentea->GetMaximumBin();
4163 if (binmax <= 1) {
4164 binmax = 2;
4165 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4166 }
4167 if (binmax >= nbins) {
4168 binmax = nbins-1;
4169 put = kFALSE;
4170 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4171 }
4172 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4173 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4174 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4175 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4176 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4177 if (TMath::Abs(l3P2am) > 0.00000001) {
4178 fPhd[0] = -(l3P1am / (2 * l3P2am));
4179 }
4180 if(!fTakeTheMaxPH){
4181 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4182 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4183 }
4184 }
4185 // Amplification region
4186 binmax = 0;
4187 ju = 0;
4188 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4189 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4190 binmax = kbin;
4191 ju = 1;
4192 }
4193 }
4194 if (binmax <= 1) {
4195 binmax = 2;
4196 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4197 }
4198 if (binmax >= nbins) {
4199 binmax = nbins-1;
4200 put = kFALSE;
4201 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4202 }
4203 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4204 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4205 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4206 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4207 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4208 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4209 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4210 }
4211 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4212 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4213 }
4214 if(fTakeTheMaxPH){
4215 fCurrentCoefE2 = fCurrentCoefE;
4216 }
4217 // Drift region
4218 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4219 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4220 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4221 }
4222 binmin = 0;
4223 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4224 if (binmin <= 1) {
4225 binmin = 2;
4226 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4227 }
4228 if (binmin >= nbins) {
4229 binmin = nbins-1;
4230 put = kFALSE;
4231 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4232 }
4233 pente->Fit("pol2"
4234 ,"0MR"
4235 ,""
4236 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4237 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4238 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4239 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4240 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4241 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4242 if (TMath::Abs(l3P2dr) > 0.00000001) {
4243 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4244 }
4245 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4246 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4247 }
4248 Float_t fPhdt0 = 0.0;
4249 Float_t t0Shift = 0.0;
4250 if(fTakeTheMaxPH) {
4251 fPhdt0 = fPhd[1];
4252 t0Shift = fT0Shift1;
4253 }
4254 else {
4255 fPhdt0 = fPhd[0];
4256 t0Shift = fT0Shift0;
4257 }
4258
4259 if ((fPhd[2] > fPhd[0]) &&
4260 (fPhd[2] > fPhd[1]) &&
4261 (fPhd[1] > fPhd[0]) &&
4262 (put)) {
4263 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4264 fNumberFitSuccess++;
4265
4266 if (fPhdt0 >= 0.0) {
4267 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4268 if (fCurrentCoef2[0] < -1.0) {
4269 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4270 }
4271 }
4272 else {
4273 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4274 }
4275
4276 }
4277 else {
4278 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4279 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4280 }
4281
4282 if (fDebugLevel == 1) {
4283 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4284 cpentei->cd();
4285 projPH->Draw();
4286 line->SetLineColor(2);
4287 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4288 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4289 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4290 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4291 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4292 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4293 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4294 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4295 cpentei2->cd();
4296 pentea->Draw();
4297 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4298 cpentei3->cd();
4299 pente->Draw();
4300 }
4301 else {
4302 delete pentea;
4303 delete pente;
4304 }
4305}
4306//_____________________________________________________________________________
4307void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4308{
4309 //
4310 // Slope methode but with polynomes de Lagrange
4311 //
4312
4313 // Constants
4314 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4315 Int_t binmax = 0;
4316 Int_t binmin = 0;
4317 //Double_t *x = new Double_t[5];
4318 //Double_t *y = new Double_t[5];
4319 Double_t x[5];
4320 Double_t y[5];
4321 x[0] = 0.0;
4322 x[1] = 0.0;
4323 x[2] = 0.0;
4324 x[3] = 0.0;
4325 x[4] = 0.0;
4326 y[0] = 0.0;
4327 y[1] = 0.0;
4328 y[2] = 0.0;
4329 y[3] = 0.0;
4330 y[4] = 0.0;
4331 fPhd[0] = 0.0;
4332 fPhd[1] = 0.0;
4333 fPhd[2] = 0.0;
4334 Int_t ju = 0;
4335 fCurrentCoefE = 0.0;
4336 fCurrentCoefE2 = 1.0;
4337 fCurrentCoef[0] = 0.0;
4338 fCurrentCoef2[0] = 0.0;
4339 TLine *line = new TLine();
4340 TF1 * polynome = 0x0;
4341 TF1 * polynomea = 0x0;
4342 TF1 * polynomeb = 0x0;
4343 Double_t c0 = 0.0;
4344 Double_t c1 = 0.0;
4345 Double_t c2 = 0.0;
4346 Double_t c3 = 0.0;
4347 Double_t c4 = 0.0;
4348
4349 // Some variables
4350 TAxis *xpph = projPH->GetXaxis();
4351 Int_t nbins = xpph->GetNbins();
4352 Double_t lowedge = xpph->GetBinLowEdge(1);
4353 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4354 Double_t widbins = (upedge - lowedge) / nbins;
4355 Double_t limit = upedge + 0.5 * widbins;
4356
4357
4358 Bool_t put = kTRUE;
4359
4360 // Beginning of the signal
4361 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4362 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4363 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4364 }
4365
4366 binmax = (Int_t) pentea->GetMaximumBin();
4367
4368 Double_t minnn = 0.0;
4369 Double_t maxxx = 0.0;
4370
4371 Int_t kase = nbins-binmax;
4372
4373 switch(kase)
4374 {
4375 case 0:
4376 put = kFALSE;
4377 break;
4378 case 1:
4379 minnn = pentea->GetBinCenter(binmax-2);
4380 maxxx = pentea->GetBinCenter(binmax);
4381 x[0] = pentea->GetBinCenter(binmax-2);
4382 x[1] = pentea->GetBinCenter(binmax-1);
4383 x[2] = pentea->GetBinCenter(binmax);
4384 y[0] = pentea->GetBinContent(binmax-2);
4385 y[1] = pentea->GetBinContent(binmax-1);
4386 y[2] = pentea->GetBinContent(binmax);
4387 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4388 AliInfo("At the limit for beginning!");
4389 break;
4390 case 2:
4391 minnn = pentea->GetBinCenter(binmax-2);
4392 maxxx = pentea->GetBinCenter(binmax+1);
4393 x[0] = pentea->GetBinCenter(binmax-2);
4394 x[1] = pentea->GetBinCenter(binmax-1);
4395 x[2] = pentea->GetBinCenter(binmax);
4396 x[3] = pentea->GetBinCenter(binmax+1);
4397 y[0] = pentea->GetBinContent(binmax-2);
4398 y[1] = pentea->GetBinContent(binmax-1);
4399 y[2] = pentea->GetBinContent(binmax);
4400 y[3] = pentea->GetBinContent(binmax+1);
4401 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4402 break;
4403 default:
4404 switch(binmax){
4405 case 0:
4406 put = kFALSE;
4407 break;
4408 case 1:
4409 minnn = pentea->GetBinCenter(binmax);
4410 maxxx = pentea->GetBinCenter(binmax+2);
4411 x[0] = pentea->GetBinCenter(binmax);
4412 x[1] = pentea->GetBinCenter(binmax+1);
4413 x[2] = pentea->GetBinCenter(binmax+2);
4414 y[0] = pentea->GetBinContent(binmax);
4415 y[1] = pentea->GetBinContent(binmax+1);
4416 y[2] = pentea->GetBinContent(binmax+2);
4417 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4418 break;
4419 case 2:
4420 minnn = pentea->GetBinCenter(binmax-1);
4421 maxxx = pentea->GetBinCenter(binmax+2);
4422 x[0] = pentea->GetBinCenter(binmax-1);
4423 x[1] = pentea->GetBinCenter(binmax);
4424 x[2] = pentea->GetBinCenter(binmax+1);
4425 x[3] = pentea->GetBinCenter(binmax+2);
4426 y[0] = pentea->GetBinContent(binmax-1);
4427 y[1] = pentea->GetBinContent(binmax);
4428 y[2] = pentea->GetBinContent(binmax+1);
4429 y[3] = pentea->GetBinContent(binmax+2);
4430 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4431 break;
4432 default:
4433 minnn = pentea->GetBinCenter(binmax-2);
4434 maxxx = pentea->GetBinCenter(binmax+2);
4435 x[0] = pentea->GetBinCenter(binmax-2);
4436 x[1] = pentea->GetBinCenter(binmax-1);
4437 x[2] = pentea->GetBinCenter(binmax);
4438 x[3] = pentea->GetBinCenter(binmax+1);
4439 x[4] = pentea->GetBinCenter(binmax+2);
4440 y[0] = pentea->GetBinContent(binmax-2);
4441 y[1] = pentea->GetBinContent(binmax-1);
4442 y[2] = pentea->GetBinContent(binmax);
4443 y[3] = pentea->GetBinContent(binmax+1);
4444 y[4] = pentea->GetBinContent(binmax+2);
4445 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4446 break;
4447 }
4448 break;
4449 }
4450
4451
4452 if(put) {
4453 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4454 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4455
4456 Double_t step = (maxxx-minnn)/10000;
4457 Double_t l = minnn;
4458 Double_t maxvalue = 0.0;
4459 Double_t placemaximum = minnn;
4460 for(Int_t o = 0; o < 10000; o++){
4461 if(o == 0) maxvalue = polynomeb->Eval(l);
4462 if(maxvalue < (polynomeb->Eval(l))){
4463 maxvalue = polynomeb->Eval(l);
4464 placemaximum = l;
4465 }
4466 l += step;
4467 }
4468 fPhd[0] = placemaximum;
4469 }
4470
4471 // Amplification region
4472 binmax = 0;
4473 ju = 0;
4474 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4475 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4476 binmax = kbin;
4477 ju = 1;
4478 }
4479 }
4480
4481 Double_t minn = 0.0;
4482 Double_t maxx = 0.0;
4483 x[0] = 0.0;
4484 x[1] = 0.0;
4485 x[2] = 0.0;
4486 x[3] = 0.0;
4487 x[4] = 0.0;
4488 y[0] = 0.0;
4489 y[1] = 0.0;
4490 y[2] = 0.0;
4491 y[3] = 0.0;
4492 y[4] = 0.0;
4493
4494 Int_t kase1 = nbins - binmax;
4495
4496 //Determination of minn and maxx
4497 //case binmax = nbins
4498 //pol2
4499 switch(kase1)
4500 {
4501 case 0:
4502 minn = projPH->GetBinCenter(binmax-2);
4503 maxx = projPH->GetBinCenter(binmax);
4504 x[0] = projPH->GetBinCenter(binmax-2);
4505 x[1] = projPH->GetBinCenter(binmax-1);
4506 x[2] = projPH->GetBinCenter(binmax);
4507 y[0] = projPH->GetBinContent(binmax-2);
4508 y[1] = projPH->GetBinContent(binmax-1);
4509 y[2] = projPH->GetBinContent(binmax);
4510 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4511 //AliInfo("At the limit for the drift!");
4512 break;
4513 case 1:
4514 minn = projPH->GetBinCenter(binmax-2);
4515 maxx = projPH->GetBinCenter(binmax+1);
4516 x[0] = projPH->GetBinCenter(binmax-2);
4517 x[1] = projPH->GetBinCenter(binmax-1);
4518 x[2] = projPH->GetBinCenter(binmax);
4519 x[3] = projPH->GetBinCenter(binmax+1);
4520 y[0] = projPH->GetBinContent(binmax-2);
4521 y[1] = projPH->GetBinContent(binmax-1);
4522 y[2] = projPH->GetBinContent(binmax);
4523 y[3] = projPH->GetBinContent(binmax+1);
4524 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4525 break;
4526 default:
4527 switch(binmax)
4528 {
4529 case 0:
4530 put = kFALSE;
4531 break;
4532 case 1:
4533 minn = projPH->GetBinCenter(binmax);
4534 maxx = projPH->GetBinCenter(binmax+2);
4535 x[0] = projPH->GetBinCenter(binmax);
4536 x[1] = projPH->GetBinCenter(binmax+1);
4537 x[2] = projPH->GetBinCenter(binmax+2);
4538 y[0] = projPH->GetBinContent(binmax);
4539 y[1] = projPH->GetBinContent(binmax+1);
4540 y[2] = projPH->GetBinContent(binmax+2);
4541 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4542 break;
4543 case 2:
4544 minn = projPH->GetBinCenter(binmax-1);
4545 maxx = projPH->GetBinCenter(binmax+2);
4546 x[0] = projPH->GetBinCenter(binmax-1);
4547 x[1] = projPH->GetBinCenter(binmax);
4548 x[2] = projPH->GetBinCenter(binmax+1);
4549 x[3] = projPH->GetBinCenter(binmax+2);
4550 y[0] = projPH->GetBinContent(binmax-1);
4551 y[1] = projPH->GetBinContent(binmax);
4552 y[2] = projPH->GetBinContent(binmax+1);
4553 y[3] = projPH->GetBinContent(binmax+2);
4554 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4555 break;
4556 default:
4557 minn = projPH->GetBinCenter(binmax-2);
4558 maxx = projPH->GetBinCenter(binmax+2);
4559 x[0] = projPH->GetBinCenter(binmax-2);
4560 x[1] = projPH->GetBinCenter(binmax-1);
4561 x[2] = projPH->GetBinCenter(binmax);
4562 x[3] = projPH->GetBinCenter(binmax+1);
4563 x[4] = projPH->GetBinCenter(binmax+2);
4564 y[0] = projPH->GetBinContent(binmax-2);
4565 y[1] = projPH->GetBinContent(binmax-1);
4566 y[2] = projPH->GetBinContent(binmax);
4567 y[3] = projPH->GetBinContent(binmax+1);
4568 y[4] = projPH->GetBinContent(binmax+2);
4569 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4570 break;
4571 }
4572 break;
4573 }
4574
4575 if(put) {
4576 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4577 polynomea->SetParameters(c0,c1,c2,c3,c4);
4578
4579 Double_t step = (maxx-minn)/1000;
4580 Double_t l = minn;
4581 Double_t maxvalue = 0.0;
4582 Double_t placemaximum = minn;
4583 for(Int_t o = 0; o < 1000; o++){
4584 if(o == 0) maxvalue = polynomea->Eval(l);
4585 if(maxvalue < (polynomea->Eval(l))){
4586 maxvalue = polynomea->Eval(l);
4587 placemaximum = l;
4588 }
4589 l += step;
4590 }
4591 fPhd[1] = placemaximum;
4592 }
4593
4594 // Drift region
4595 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4596 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4597 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4598 }
4599 binmin = 0;
4600 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4601
4602 //should not happen
4603 if (binmin <= 1) {
4604 binmin = 2;
4605 put = 1;
4606 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4607 }
4608
4609 //check
4610 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4611 AliInfo("Too many fluctuations at the end!");
4612 put = kFALSE;
4613 }
4614 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4615 AliInfo("Too many fluctuations at the end!");
4616 put = kFALSE;
4617 }
4618 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4619 AliInfo("No entries for the next bin!");
4620 pente->SetBinContent(binmin,0);
4621 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4622 }
4623
4624
4625 x[0] = 0.0;
4626 x[1] = 0.0;
4627 x[2] = 0.0;
4628 x[3] = 0.0;
4629 x[4] = 0.0;
4630 y[0] = 0.0;
4631 y[1] = 0.0;
4632 y[2] = 0.0;
4633 y[3] = 0.0;
4634 y[4] = 0.0;
4635 Double_t min = 0.0;
4636 Double_t max = 0.0;
4637 Bool_t case1 = kFALSE;
4638 Bool_t case2 = kFALSE;
4639 Bool_t case4 = kFALSE;
4640
4641 //Determination of min and max
4642 //case binmin <= nbins-3
4643 //pol4 case 3
4644 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4645 min = pente->GetBinCenter(binmin-2);
4646 max = pente->GetBinCenter(binmin+2);
4647 x[0] = pente->GetBinCenter(binmin-2);
4648 x[1] = pente->GetBinCenter(binmin-1);
4649 x[2] = pente->GetBinCenter(binmin);
4650 x[3] = pente->GetBinCenter(binmin+1);
4651 x[4] = pente->GetBinCenter(binmin+2);
4652 y[0] = pente->GetBinContent(binmin-2);
4653 y[1] = pente->GetBinContent(binmin-1);
4654 y[2] = pente->GetBinContent(binmin);
4655 y[3] = pente->GetBinContent(binmin+1);
4656 y[4] = pente->GetBinContent(binmin+2);
4657 //Calcul the polynome de Lagrange
4658 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4659 //richtung +/-
4660 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4661 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4662 //AliInfo("polynome 4 false 1");
4663 put = kFALSE;
4664 }
4665 if(((binmin+3) <= (nbins-1)) &&
4666 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4667 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4668 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4669 AliInfo("polynome 4 false 2");
4670 put = kFALSE;
4671 }
4672 // poly 3
4673 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4674 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4675 //AliInfo("polynome 4 case 1");
4676 case1 = kTRUE;
4677 }
4678 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4679 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4680 //AliInfo("polynome 4 case 4");
4681 case4 = kTRUE;
4682 }
4683
4684 }
4685 //case binmin = nbins-2
4686 //pol3 case 1
4687 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4688 (case1)){
4689 min = pente->GetBinCenter(binmin-2);
4690 max = pente->GetBinCenter(binmin+1);
4691 x[0] = pente->GetBinCenter(binmin-2);
4692 x[1] = pente->GetBinCenter(binmin-1);
4693 x[2] = pente->GetBinCenter(binmin);
4694 x[3] = pente->GetBinCenter(binmin+1);
4695 y[0] = pente->GetBinContent(binmin-2);
4696 y[1] = pente->GetBinContent(binmin-1);
4697 y[2] = pente->GetBinContent(binmin);
4698 y[3] = pente->GetBinContent(binmin+1);
4699 //Calcul the polynome de Lagrange
4700 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4701 //richtung +: nothing
4702 //richtung -
4703 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4704 //AliInfo("polynome 3- case 2");
4705 case2 = kTRUE;
4706 }
4707 }
4708 //pol3 case 4
4709 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4710 (case4)){
4711 min = pente->GetBinCenter(binmin-1);
4712 max = pente->GetBinCenter(binmin+2);
4713 x[0] = pente->GetBinCenter(binmin-1);
4714 x[1] = pente->GetBinCenter(binmin);
4715 x[2] = pente->GetBinCenter(binmin+1);
4716 x[3] = pente->GetBinCenter(binmin+2);
4717 y[0] = pente->GetBinContent(binmin-1);
4718 y[1] = pente->GetBinContent(binmin);
4719 y[2] = pente->GetBinContent(binmin+1);
4720 y[3] = pente->GetBinContent(binmin+2);
4721 //Calcul the polynome de Lagrange
4722 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4723 //richtung +
4724 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4725 //AliInfo("polynome 3+ case 2");
4726 case2 = kTRUE;
4727 }
4728 }
4729 //pol2 case 5
4730 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4731 min = pente->GetBinCenter(binmin);
4732 max = pente->GetBinCenter(binmin+2);
4733 x[0] = pente->GetBinCenter(binmin);
4734 x[1] = pente->GetBinCenter(binmin+1);
4735 x[2] = pente->GetBinCenter(binmin+2);
4736 y[0] = pente->GetBinContent(binmin);
4737 y[1] = pente->GetBinContent(binmin+1);
4738 y[2] = pente->GetBinContent(binmin+2);
4739 //Calcul the polynome de Lagrange
4740 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4741 //richtung +
4742 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4743 //AliInfo("polynome 2+ false");
4744 put = kFALSE;
4745 }
4746 }
4747 //pol2 case 2
4748 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4749 (case2)){
4750 min = pente->GetBinCenter(binmin-1);
4751 max = pente->GetBinCenter(binmin+1);
4752 x[0] = pente->GetBinCenter(binmin-1);
4753 x[1] = pente->GetBinCenter(binmin);
4754 x[2] = pente->GetBinCenter(binmin+1);
4755 y[0] = pente->GetBinContent(binmin-1);
4756 y[1] = pente->GetBinContent(binmin);
4757 y[2] = pente->GetBinContent(binmin+1);
4758 //Calcul the polynome de Lagrange
4759 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4760 //richtung +: nothing
4761 //richtung -: nothing
4762 }
4763 //case binmin = nbins-1
4764 //pol2 case 0
4765 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4766 min = pente->GetBinCenter(binmin-2);
4767 max = pente->GetBinCenter(binmin);
4768 x[0] = pente->GetBinCenter(binmin-2);
4769 x[1] = pente->GetBinCenter(binmin-1);
4770 x[2] = pente->GetBinCenter(binmin);
4771 y[0] = pente->GetBinContent(binmin-2);
4772 y[1] = pente->GetBinContent(binmin-1);
4773 y[2] = pente->GetBinContent(binmin);
4774 //Calcul the polynome de Lagrange
4775 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4776 //AliInfo("At the limit for the drift!");
4777 //fluctuation too big!
4778 //richtung +: nothing
4779 //richtung -
4780 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4781 //AliInfo("polynome 2- false ");
4782 put = kFALSE;
4783 }
4784 }
4785 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4786 put = kFALSE;
4787 AliInfo("At the limit for the drift and not usable!");
4788 }
4789
4790 //pass
4791 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4792 put = kFALSE;
4793 AliInfo("For the drift...problem!");
4794 }
4795 //pass but should not happen
4796 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4797 put = kFALSE;
4798 AliInfo("For the drift...problem!");
4799 }
4800
4801 if(put) {
4802 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4803 polynome->SetParameters(c0,c1,c2,c3,c4);
4804 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4805 Double_t step = (max-min)/1000;
4806 Double_t l = min;
4807 Double_t minvalue = 0.0;
4808 Double_t placeminimum = min;
4809 for(Int_t o = 0; o < 1000; o++){
4810 if(o == 0) minvalue = polynome->Eval(l);
4811 if(minvalue > (polynome->Eval(l))){
4812 minvalue = polynome->Eval(l);
4813 placeminimum = l;
4814 }
4815 l += step;
4816 }
4817 fPhd[2] = placeminimum;
4818 }
4819 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4820 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4821 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4822
4823 Float_t fPhdt0 = 0.0;
4824 Float_t t0Shift = 0.0;
4825 if(fTakeTheMaxPH) {
4826 fPhdt0 = fPhd[1];
4827 t0Shift = fT0Shift1;
4828 }
4829 else {
4830 fPhdt0 = fPhd[0];
4831 t0Shift = fT0Shift0;
4832 }
4833
4834 if ((fPhd[2] > fPhd[0]) &&
4835 (fPhd[2] > fPhd[1]) &&
4836 (fPhd[1] > fPhd[0]) &&
4837 (put)) {
4838 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4839 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4840 else fNumberFitSuccess++;
4841 if (fPhdt0 >= 0.0) {
4842 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4843 if (fCurrentCoef2[0] < -1.0) {
4844 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4845 }
4846 }
4847 else {
4848 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4849 }
4850 }
4851 else {
4852 //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4853 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4854
4855 if((fPhd[1] > fPhd[0]) &&
4856 (put)) {
4857 if (fPhdt0 >= 0.0) {
4858 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4859 if (fCurrentCoef2[0] < -1.0) {
4860 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4861 }
4862 }
4863 else {
4864 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4865 }
4866 }
4867 else{
4868 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4869 //printf("Fit failed!\n");
4870 }
4871 }
4872
4873 if (fDebugLevel == 1) {
4874 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4875 cpentei->cd();
4876 projPH->Draw();
4877 line->SetLineColor(2);
4878 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4879 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4880 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4881 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4882 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4883 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4884 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4885 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4886 cpentei2->cd();
4887 pentea->Draw();
4888 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4889 cpentei3->cd();
4890 pente->Draw();
4891 }
4892 else {
4893 delete pentea;
4894 delete pente;
4895 if(polynome) delete polynome;
4896 if(polynomea) delete polynomea;
4897 if(polynomeb) delete polynomeb;
4898 //if(x) delete [] x;
4899 //if(y) delete [] y;
4900 if(line) delete line;
4901
4902 }
4903
4904 //Provisoire
4905 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4906 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4907
4908 projPH->SetDirectory(0);
4909
4910}
4911
4912//_____________________________________________________________________________
4913void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4914{
4915 //
4916 // Fit methode for the drift velocity
4917 //
4918
4919 // Constants
4920 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4921
4922 // Some variables
4923 TAxis *xpph = projPH->GetXaxis();
4924 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4925
4926 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4927 fPH->SetParameter(0,0.469); // Scaling
4928 fPH->SetParameter(1,0.18); // Start
4929 fPH->SetParameter(2,0.0857325); // AR
4930 fPH->SetParameter(3,1.89); // DR
4931 fPH->SetParameter(4,0.08); // QA/QD
4932 fPH->SetParameter(5,0.0); // Baseline
4933
4934 TLine *line = new TLine();
4935
4936 fCurrentCoef[0] = 0.0;
4937 fCurrentCoef2[0] = 0.0;
4938 fCurrentCoefE = 0.0;
4939 fCurrentCoefE2 = 0.0;
4940
4941 if (idect%fFitPHPeriode == 0) {
4942
4943 AliInfo(Form("The detector %d will be fitted",idect));
4944 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4945 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4946 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4947 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4948 fPH->SetParameter(4,0.225); // QA/QD
4949 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4950
4951 if (fDebugLevel != 1) {
4952 projPH->Fit(fPH,"0M","",0.0,upedge);
4953 }
4954 else {
4955 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4956 cpente->cd();
4957 projPH->Fit(fPH,"M+","",0.0,upedge);
4958 projPH->Draw("E0");
4959 line->SetLineColor(4);
4960 line->DrawLine(fPH->GetParameter(1)
4961 ,0
4962 ,fPH->GetParameter(1)
4963 ,projPH->GetMaximum());
4964 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4965 ,0
4966 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4967 ,projPH->GetMaximum());
4968 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4969 ,0
4970 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4971 ,projPH->GetMaximum());
4972 }
4973
4974 if (fPH->GetParameter(3) != 0) {
4975 fNumberFitSuccess++;
4976 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4977 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4978 fCurrentCoef2[0] = fPH->GetParameter(1);
4979 fCurrentCoefE2 = fPH->GetParError(1);
4980 }
4981 else {
4982 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4983 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4984 }
4985
4986 }
4987 else {
4988
4989 // Put the default value
4990 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4991 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4992 }
4993
4994 if (fDebugLevel != 1) {
4995 delete fPH;
4996 }
4997
4998}
4999//_____________________________________________________________________________
5000Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5001{
5002 //
5003 // Fit methode for the sigma of the pad response function
5004 //
5005
5006 TVectorD param(3);
5007
5008 fCurrentCoef[0] = 0.0;
5009 fCurrentCoefE = 0.0;
5010
5011 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
5012
5013 if(TMath::Abs(ret+4) <= 0.000000001){
5014 fCurrentCoef[0] = -fCurrentCoef[1];
5015 return kFALSE;
5016 }
5017 else {
5018 fNumberFitSuccess++;
5019 fCurrentCoef[0] = param[2];
5020 fCurrentCoefE = ret;
5021 return kTRUE;
5022 }
5023}
5024//_____________________________________________________________________________
5025Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, Bool_t bError)
5026{
5027 //
5028 // Fit methode for the sigma of the pad response function
5029 //
5030
5031 //We should have at least 3 points
5032 if(nBins <=3) return -4.0;
5033
5034 TLinearFitter fitter(3,"pol2");
5035 fitter.StoreData(kFALSE);
5036 fitter.ClearPoints();
5037 TVectorD par(3);
5038 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5039 Float_t entries = 0;
5040 Int_t nbbinwithentries = 0;
5041 for (Int_t i=0; i<nBins; i++){
5042 entries+=arraye[i];
5043 if(arraye[i] > 15) nbbinwithentries++;
5044 //printf("entries for i %d: %f\n",i,arraye[i]);
5045 }
5046 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5047 //printf("entries %f\n",entries);
5048 //printf("nbbinwithentries %d\n",nbbinwithentries);
5049
5050 Int_t npoints=0;
5051 Float_t errorm = 0.0;
5052 Float_t errorn = 0.0;
5053 Float_t error = 0.0;
5054
5055 //
5056 for (Int_t ibin=0;ibin<nBins; ibin++){
5057 Float_t entriesI = arraye[ibin];
5058 Float_t valueI = arraym[ibin];
5059 Double_t xcenter = 0.0;
5060 Float_t val = 0.0;
5061 if ((entriesI>15) && (valueI>0.0)){
5062 xcenter = xMin+(ibin+0.5)*binWidth;
5063 errorm = 0.0;
5064 errorn = 0.0;
5065 error = 0.0;
5066 if(!bError){
5067 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5068 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5069 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5070 }
5071 else{
5072 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5073 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5074 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5075 }
5076 if(TMath::Abs(error) < 0.000000001) continue;
5077 val = TMath::Log(Float_t(valueI));
5078 fitter.AddPoint(&xcenter,val,error);
5079 npoints++;
5080 }
5081
5082 if(fDebugLevel > 1){
5083
5084 if ( !fDebugStreamer ) {
5085 //debug stream
5086 TDirectory *backup = gDirectory;
5087 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5088 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5089 }
5090
5091 Int_t detector = fCountDet;
5092 Int_t layer = GetLayer(fCountDet);
5093 Int_t group = ibin;
5094
5095 (* fDebugStreamer) << "FitGausMIFill"<<
5096 "detector="<<detector<<
5097 "layer="<<layer<<
5098 "nbins="<<nBins<<
5099 "group="<<group<<
5100 "entriesI="<<entriesI<<
5101 "valueI="<<valueI<<
5102 "val="<<val<<
5103 "xcenter="<<xcenter<<
5104 "errorm="<<errorm<<
5105 "errorn="<<errorn<<
5106 "error="<<error<<
5107 "bError="<<bError<<
5108 "\n";
5109 }
5110
5111 }
5112
5113 if(npoints <=3) return -4.0;
5114
5115 Double_t chi2 = 0;
5116 if (npoints>3){
5117 fitter.Eval();
5118 fitter.GetParameters(par);
5119 chi2 = fitter.GetChisquare()/Float_t(npoints);
5120
5121
5122 if (!param) param = new TVectorD(3);
5123 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5124 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5125 Double_t deltax = (fitter.GetParError(2))/x;
5126 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5127 chi2 = errorparam2;
5128
5129 (*param)[1] = par[1]/(-2.*par[2]);
5130 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5131 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5132 if ( lnparam0>307 ) return -4;
5133 (*param)[0] = TMath::Exp(lnparam0);
5134
5135 if(fDebugLevel > 1){
5136
5137 if ( !fDebugStreamer ) {
5138 //debug stream
5139 TDirectory *backup = gDirectory;
5140 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5141 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5142 }
5143
5144 Int_t detector = fCountDet;
5145 Int_t layer = GetLayer(fCountDet);
5146
5147
5148 (* fDebugStreamer) << "FitGausMIFit"<<
5149 "detector="<<detector<<
5150 "layer="<<layer<<
5151 "nbins="<<nBins<<
5152 "errorsigma="<<chi2<<
5153 "mean="<<(*param)[1]<<
5154 "sigma="<<(*param)[2]<<
5155 "constant="<<(*param)[0]<<
5156 "\n";
5157 }
5158 }
5159
5160 if((chi2/(*param)[2]) > 0.1){
5161 if(bError){
5162 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5163 }
5164 else return -4.0;
5165 }
5166
5167 if(fDebugLevel == 1){
5168 TString name("PRF");
5169 name += (Int_t)xMin;
5170 name += (Int_t)xMax;
5171 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5172 c1->cd();
5173 name += "histo";
5174 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5175 for(Int_t k = 0; k < nBins; k++){
5176 histo->SetBinContent(k+1,arraym[k]);
5177 histo->SetBinError(k+1,arrayme[k]);
5178 }
5179 histo->Draw();
5180 name += "functionf";
5181 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5182 f1->SetParameter(0, (*param)[0]);
5183 f1->SetParameter(1, (*param)[1]);
5184 f1->SetParameter(2, (*param)[2]);
5185 f1->Draw("same");
5186 }
5187
5188
5189 return chi2;
5190
5191}
5192//_____________________________________________________________________________
5193void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5194{
5195 //
5196 // Fit methode for the sigma of the pad response function
5197 //
5198
5199 fCurrentCoef[0] = 0.0;
5200 fCurrentCoefE = 0.0;
5201
5202 if (fDebugLevel != 1) {
5203 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5204 }
5205 else {
5206 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5207 cfit->cd();
5208 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5209 projPRF->Draw();
5210 }
5211 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5212 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5213 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5214 else {
5215 fNumberFitSuccess++;
5216 }
5217}
5218//_____________________________________________________________________________
5219void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5220{
5221 //
5222 // Fit methode for the sigma of the pad response function
5223 //
5224 fCurrentCoef[0] = 0.0;
5225 fCurrentCoefE = 0.0;
5226 if (fDebugLevel == 1) {
5227 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5228 cfit->cd();
5229 projPRF->Draw();
5230 }
5231 fCurrentCoef[0] = projPRF->GetRMS();
5232 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5233 else {
5234 fNumberFitSuccess++;
5235 }
5236}
5237//_____________________________________________________________________________
5238void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5239{
5240 //
5241 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5242 //
5243
5244 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5245
5246
5247 Int_t nbins = (Int_t)(nybins/(2*nbg));
5248 Float_t lowedge = -3.0*nbg;
5249 Float_t upedge = lowedge + 3.0;
5250 Int_t offset = 0;
5251 Int_t npoints = 0;
5252 Double_t xvalues = -0.2*nbg+0.1;
5253 Double_t y = 0.0;
5254 Int_t total = 2*nbg;
5255
5256
5257 for(Int_t k = 0; k < total; k++){
5258 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5259 npoints++;
5260 y = fCurrentCoef[0]*fCurrentCoef[0];
5261 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5262 }
5263
5264 if(fDebugLevel > 1){
5265
5266 if ( !fDebugStreamer ) {
5267 //debug stream
5268 TDirectory *backup = gDirectory;
5269 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5270 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5271 }
5272
5273 Int_t detector = fCountDet;
5274 Int_t layer = GetLayer(fCountDet);
5275 Int_t nbtotal = total;
5276 Int_t group = k;
5277 Float_t low = lowedge;
5278 Float_t up = upedge;
5279 Float_t tnp = xvalues;
5280 Float_t wid = fCurrentCoef[0];
5281 Float_t widfE = fCurrentCoefE;
5282
5283 (* fDebugStreamer) << "FitTnpRange0"<<
5284 "detector="<<detector<<
5285 "layer="<<layer<<
5286 "nbtotal="<<nbtotal<<
5287 "group="<<group<<
5288 "low="<<low<<
5289 "up="<<up<<
5290 "offset="<<offset<<
5291 "tnp="<<tnp<<
5292 "wid="<<wid<<
5293 "widfE="<<widfE<<
5294 "\n";
5295 }
5296
5297 offset += nbins;
5298 lowedge += 3.0;
5299 upedge += 3.0;
5300 xvalues += 0.2;
5301
5302 }
5303
5304 fCurrentCoefE = 0.0;
5305 fCurrentCoef[0] = 0.0;
5306
5307 //printf("npoints\n",npoints);
5308
5309 if(npoints < 3){
5310 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5311 }
5312 else{
5313
5314 TVectorD pars0;
5315 linearfitter.Eval();
5316 linearfitter.GetParameters(pars0);
5317 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5318 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5319 Double_t min0 = 0.0;
5320 Double_t ermin0 = 0.0;
5321 //Double_t prfe0 = 0.0;
5322 Double_t prf0 = 0.0;
5323 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5324 min0 = -pars0[1]/(2*pars0[2]);
5325 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5326 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5327 if(prf0 > 0.0) {
5328 /*
5329 prfe0 = linearfitter->GetParError(0)*pointError0
5330 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5331 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5332 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5333 fCurrentCoefE = (Float_t) prfe0;
5334 */
5335 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5336 }
5337 else{
5338 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5339 }
5340 }
5341 else {
5342 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5343 }
5344
5345 if(fDebugLevel > 1){
5346
5347 if ( !fDebugStreamer ) {
5348 //debug stream
5349 TDirectory *backup = gDirectory;
5350 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5351 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5352 }
5353
5354 Int_t detector = fCountDet;
5355 Int_t layer = GetLayer(fCountDet);
5356 Int_t nbtotal = total;
5357 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5358 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5359
5360 (* fDebugStreamer) << "FitTnpRange1"<<
5361 "detector="<<detector<<
5362 "layer="<<layer<<
5363 "nbtotal="<<nbtotal<<
5364 "par0="<<pars0[0]<<
5365 "par1="<<pars0[1]<<
5366 "par2="<<pars0[2]<<
5367 "npoints="<<npoints<<
5368 "sigmax="<<sigmax<<
5369 "tan="<<min0<<
5370 "sigmaprf="<<fCurrentCoef[0]<<
5371 "sigprf="<<fCurrentCoef[1]<<
5372 "\n";
5373 }
5374
5375 }
5376
5377}
5378//_____________________________________________________________________________
5379void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5380{
5381 //
5382 // Only mean methode for the gain factor
5383 //
5384
5385 fCurrentCoef[0] = mean;
5386 fCurrentCoefE = 0.0;
5387 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5388 if (fDebugLevel == 1) {
5389 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5390 cpmean->cd();
5391 projch->Draw();
5392 }
5393 CalculChargeCoefMean(kTRUE);
5394 fNumberFitSuccess++;
5395}
5396//_____________________________________________________________________________
5397void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5398{
5399 //
5400 // mean w methode for the gain factor
5401 //
5402
5403 //Number of bins
5404 Int_t nybins = projch->GetNbinsX();
5405
5406 //The weight function
5407 Double_t a = 0.00228515;
5408 Double_t b = -0.00231487;
5409 Double_t c = 0.00044298;
5410 Double_t d = -0.00379239;
5411 Double_t e = 0.00338349;
5412
5413// 0 |0.00228515
5414// 1 |-0.00231487
5415// 2 |0.00044298
5416// 3 |-0.00379239
5417// 4 |0.00338349
5418
5419
5420
5421 //A arbitrary error for the moment
5422 fCurrentCoefE = 0.0;
5423 fCurrentCoef[0] = 0.0;
5424
5425 //Calcul
5426 Double_t sumw = 0.0;
5427 Double_t sum = 0.0;
5428 Float_t sumAll = (Float_t) nentries;
5429 Int_t sumCurrent = 0;
5430 for(Int_t k = 0; k <nybins; k++){
5431 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5432 if (fraction>0.95) break;
5433 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5434 e*fraction*fraction*fraction*fraction;
5435 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5436 sum += weight*projch->GetBinContent(k+1);
5437 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5438 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5439 }
5440 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5441
5442 if (fDebugLevel == 1) {
5443 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5444 cpmeanw->cd();
5445 projch->Draw();
5446 }
5447 fNumberFitSuccess++;
5448 CalculChargeCoefMean(kTRUE);
5449}
5450//_____________________________________________________________________________
5451void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5452{
5453 //
5454 // mean w methode for the gain factor
5455 //
5456
5457 //Number of bins
5458 Int_t nybins = projch->GetNbinsX();
5459
5460 //The weight function
5461 Double_t a = 0.00228515;
5462 Double_t b = -0.00231487;
5463 Double_t c = 0.00044298;
5464 Double_t d = -0.00379239;
5465 Double_t e = 0.00338349;
5466
5467// 0 |0.00228515
5468// 1 |-0.00231487
5469// 2 |0.00044298
5470// 3 |-0.00379239
5471// 4 |0.00338349
5472
5473
5474
5475 //A arbitrary error for the moment
5476 fCurrentCoefE = 0.0;
5477 fCurrentCoef[0] = 0.0;
5478
5479 //Calcul
5480 Double_t sumw = 0.0;
5481 Double_t sum = 0.0;
5482 Int_t sumCurrent = 0;
5483 for(Int_t k = 0; k <nybins; k++){
5484 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5485 if (fraction>0.95) break;
5486 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5487 e*fraction*fraction*fraction*fraction;
5488 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5489 sum += weight*projch->GetBinContent(k+1);
5490 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5491 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5492 }
5493 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5494
5495 if (fDebugLevel == 1) {
5496 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5497 cpmeanw->cd();
5498 projch->Draw();
5499 }
5500 fNumberFitSuccess++;
5501}
5502//_____________________________________________________________________________
5503void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5504{
5505 //
5506 // Fit methode for the gain factor
5507 //
5508
5509 fCurrentCoef[0] = 0.0;
5510 fCurrentCoefE = 0.0;
5511 Double_t chisqrl = 0.0;
5512 Double_t chisqrg = 0.0;
5513 Double_t chisqr = 0.0;
5514 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5515
5516 projch->Fit("landau","0",""
5517 ,(Double_t) mean/fBeginFitCharge
5518 ,projch->GetBinCenter(projch->GetNbinsX()));
5519 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5520 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5521 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5522 chisqrl = projch->GetFunction("landau")->GetChisquare();
5523
5524 projch->Fit("gaus","0",""
5525 ,(Double_t) mean/fBeginFitCharge
5526 ,projch->GetBinCenter(projch->GetNbinsX()));
5527 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5528 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5529 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5530
5531 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5532 if (fDebugLevel != 1) {
5533 projch->Fit("fLandauGaus","0",""
5534 ,(Double_t) mean/fBeginFitCharge
5535 ,projch->GetBinCenter(projch->GetNbinsX()));
5536 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5537 }
5538 else {
5539 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5540 cp->cd();
5541 projch->Fit("fLandauGaus","+",""
5542 ,(Double_t) mean/fBeginFitCharge
5543 ,projch->GetBinCenter(projch->GetNbinsX()));
5544 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5545 projch->Draw();
5546 fLandauGaus->Draw("same");
5547 }
5548
5549 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5550 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5551 fNumberFitSuccess++;
5552 CalculChargeCoefMean(kTRUE);
5553 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5554 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5555 }
5556 else {
5557 CalculChargeCoefMean(kFALSE);
5558 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5559 }
5560
5561 if (fDebugLevel != 1) {
5562 delete fLandauGaus;
5563 }
5564
5565}
5566//_____________________________________________________________________________
5567void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5568{
5569 //
5570 // Fit methode for the gain factor more time consuming
5571 //
5572
5573
5574 //Some parameters to initialise
5575 Double_t widthLandau, widthGaus, mPV, integral;
5576 Double_t chisquarel = 0.0;
5577 Double_t chisquareg = 0.0;
5578 projch->Fit("landau","0M+",""
5579 ,(Double_t) mean/6
5580 ,projch->GetBinCenter(projch->GetNbinsX()));
5581 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5582 chisquarel = projch->GetFunction("landau")->GetChisquare();
5583 projch->Fit("gaus","0M+",""
5584 ,(Double_t) mean/6
5585 ,projch->GetBinCenter(projch->GetNbinsX()));
5586 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5587 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5588
5589 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5590 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5591
5592 // Setting fit range and start values
5593 Double_t fr[2];
5594 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5595 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5596 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5597 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5598 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5599 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5600 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5601 fr[0] = 0.3 * mean;
5602 fr[1] = 3.0 * mean;
5603 fCurrentCoef[0] = 0.0;
5604 fCurrentCoefE = 0.0;
5605
5606 Double_t chisqr;
5607 Int_t ndf;
5608 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5609 ,&pllo[0],&plhi[0]
5610 ,&fp[0],&fpe[0]
5611 ,&chisqr,&ndf);
5612
5613 Double_t projchPeak;
5614 Double_t projchFWHM;
5615 LanGauPro(fp,projchPeak,projchFWHM);
5616
5617 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5618 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5619 fNumberFitSuccess++;
5620 CalculChargeCoefMean(kTRUE);
5621 fCurrentCoef[0] = fp[1];
5622 fCurrentCoefE = fpe[1];
5623 //chargeCoefE2 = chisqr;
5624 }
5625 else {
5626 CalculChargeCoefMean(kFALSE);
5627 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5628 }
5629 if (fDebugLevel == 1) {
5630 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5631 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5632 cpy->cd();
5633 projch->Draw();
5634 fitsnr->Draw("same");
5635 }
5636 else {
5637 delete fitsnr;
5638 }
5639}
5640//_____________________________________________________________________________
5641void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5642{
5643 //
5644 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5645 //
5646 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5647 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5648 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5649
5650 c4 = 0.0;
5651 c3 = 0.0;
5652 c2 = x0+x1+x2;
5653 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5654 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5655
5656}
5657
5658//_____________________________________________________________________________
5659void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5660{
5661 //
5662 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5663 //
5664 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5665 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5666 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5667 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5668
5669 c4 = 0.0;
5670 c3 = x0+x1+x2+x3;
5671 c2 = -(x0*(x[1]+x[2]+x[3])
5672 +x1*(x[0]+x[2]+x[3])
5673 +x2*(x[0]+x[1]+x[3])
5674 +x3*(x[0]+x[1]+x[2]));
5675 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5676 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5677 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5678 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5679
5680 c0 = -(x0*x[1]*x[2]*x[3]
5681 +x1*x[0]*x[2]*x[3]
5682 +x2*x[0]*x[1]*x[3]
5683 +x3*x[0]*x[1]*x[2]);
5684
5685
5686}
5687
5688//_____________________________________________________________________________
5689void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5690{
5691 //
5692 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5693 //
5694 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5695 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5696 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5697 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5698 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5699
5700
5701 c4 = x0+x1+x2+x3+x4;
5702 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5703 +x1*(x[0]+x[2]+x[3]+x[4])
5704 +x2*(x[0]+x[1]+x[3]+x[4])
5705 +x3*(x[0]+x[1]+x[2]+x[4])
5706 +x4*(x[0]+x[1]+x[2]+x[3]));
5707 c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5708 +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5709 +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
5710 +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
5711 +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
5712
5713 c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5714 +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
5715 +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
5716 +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
5717 +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
5718
5719 c0 = (x0*x[1]*x[2]*x[3]*x[4]
5720 +x1*x[0]*x[2]*x[3]*x[4]
5721 +x2*x[0]*x[1]*x[3]*x[4]
5722 +x3*x[0]*x[1]*x[2]*x[4]
5723 +x4*x[0]*x[1]*x[2]*x[3]);
5724
5725}
5726//_____________________________________________________________________________
5727void AliTRDCalibraFit::NormierungCharge()
5728{
5729 //
5730 // Normalisation of the gain factor resulting for the fits
5731 //
5732
5733 // Calcul of the mean of choosen method by fFitChargeNDB
5734 Double_t sum = 0.0;
5735 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5736 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5737 Int_t total = 0;
5738 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5739 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5740 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5741 if (GetStack(detector) == 2) {
5742 total = 1728;
5743 }
5744 if (GetStack(detector) != 2) {
5745 total = 2304;
5746 }
5747 for (Int_t j = 0; j < total; j++) {
5748 if (coef[j] >= 0) {
5749 sum += coef[j];
5750 }
5751 }
5752 }
5753
5754 if (sum > 0) {
5755 fScaleFitFactor = fScaleFitFactor / sum;
5756 }
5757 else {
5758 fScaleFitFactor = 1.0;
5759 }
5760
5761 //methode de boeuf mais bon...
5762 Double_t scalefactor = fScaleFitFactor;
5763
5764 if(fDebugLevel > 1){
5765
5766 if ( !fDebugStreamer ) {
5767 //debug stream
5768 TDirectory *backup = gDirectory;
5769 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5770 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5771 }
5772 (* fDebugStreamer) << "NormierungCharge"<<
5773 "scalefactor="<<scalefactor<<
5774 "\n";
5775 }
5776}
5777//_____________________________________________________________________________
5778TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5779{
5780 //
5781 // Rebin of the 1D histo for the gain calibration if needed.
5782 // you have to choose fRebin, divider of fNumberBinCharge
5783 //
5784
5785 TAxis *xhist = hist->GetXaxis();
5786 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5787 ,xhist->GetBinLowEdge(1)
5788 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5789
5790 AliInfo(Form("fRebin: %d",fRebin));
5791 Int_t i = 1;
5792 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5793 Double_t sum = 0.0;
5794 for (Int_t ji = i; ji < i+fRebin; ji++) {
5795 sum += hist->GetBinContent(ji);
5796 }
5797 sum = sum / fRebin;
5798 rehist->SetBinContent(k,sum);
5799 i += fRebin;
5800 }
5801
5802 return rehist;
5803
5804}
5805
5806//_____________________________________________________________________________
5807TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5808{
5809 //
5810 // Rebin of the 1D histo for the gain calibration if needed
5811 // you have to choose fRebin divider of fNumberBinCharge
5812 //
5813
5814 TAxis *xhist = hist->GetXaxis();
5815 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5816 ,xhist->GetBinLowEdge(1)
5817 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5818
5819 AliInfo(Form("fRebin: %d",fRebin));
5820 Int_t i = 1;
5821 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5822 Double_t sum = 0.0;
5823 for (Int_t ji = i; ji < i+fRebin; ji++) {
5824 sum += hist->GetBinContent(ji);
5825 }
5826 sum = sum/fRebin;
5827 rehist->SetBinContent(k,sum);
5828 i += fRebin;
5829 }
5830
5831 return rehist;
5832
5833}
5834//
5835//____________Some basic geometry function_____________________________________
5836//
5837
5838//_____________________________________________________________________________
5839Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5840{
5841 //
5842 // Reconstruct the plane number from the detector number
5843 //
5844
5845 return ((Int_t) (d % 6));
5846
5847}
5848
5849//_____________________________________________________________________________
5850Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5851{
5852 //
5853 // Reconstruct the stack number from the detector number
5854 //
5855 const Int_t kNlayer = 6;
5856
5857 return ((Int_t) (d % 30) / kNlayer);
5858
5859}
5860
5861//_____________________________________________________________________________
5862Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5863{
5864 //
5865 // Reconstruct the sector number from the detector number
5866 //
5867 Int_t fg = 30;
5868
5869 return ((Int_t) (d / fg));
5870
5871}
5872
5873//
5874//____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5875//
5876//_______________________________________________________________________________
5877void AliTRDCalibraFit::ResetVectorFit()
5878{
5879 //
5880 // Reset the VectorFits
5881 //
5882
5883 fVectorFit.SetOwner();
5884 fVectorFit.Clear();
5885 fVectorFit2.SetOwner();
5886 fVectorFit2.Clear();
5887
5888}
5889//
5890//____________Private Functions________________________________________________
5891//
5892
5893//_____________________________________________________________________________
5894Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5895{
5896 //
5897 // Function for the fit
5898 //
5899
5900 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5901
5902 //PARAMETERS FOR FIT PH
5903 // PASAv.4
5904 //fAsymmGauss->SetParameter(0,0.113755);
5905 //fAsymmGauss->SetParameter(1,0.350706);
5906 //fAsymmGauss->SetParameter(2,0.0604244);
5907 //fAsymmGauss->SetParameter(3,7.65596);
5908 //fAsymmGauss->SetParameter(4,1.00124);
5909 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5910
5911 Double_t xx = x[0];
5912
5913 if (xx < par[1]) {
5914 return par[5];
5915 }
5916
5917 Double_t dx = 0.005;
5918 Double_t xs = par[1];
5919 Double_t ss = 0.0;
5920 Double_t paras[2] = { 0.0, 0.0 };
5921
5922 while (xs < xx) {
5923 if ((xs >= par[1]) &&
5924 (xs < (par[1]+par[2]))) {
5925 //fAsymmGauss->SetParameter(0,par[0]);
5926 //fAsymmGauss->SetParameter(1,xs);
5927 //ss += fAsymmGauss->Eval(xx);
5928 paras[0] = par[0];
5929 paras[1] = xs;
5930 ss += AsymmGauss(&xx,paras);
5931 }
5932 if ((xs >= (par[1]+par[2])) &&
5933 (xs < (par[1]+par[2]+par[3]))) {
5934 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5935 //fAsymmGauss->SetParameter(1,xs);
5936 //ss += fAsymmGauss->Eval(xx);
5937 paras[0] = par[0]*par[4];
5938 paras[1] = xs;
5939 ss += AsymmGauss(&xx,paras);
5940 }
5941 xs += dx;
5942 }
5943
5944 return ss + par[5];
5945
5946}
5947
5948//_____________________________________________________________________________
5949Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5950{
5951 //
5952 // Function for the fit
5953 //
5954
5955 //par[0] = normalization
5956 //par[1] = mean
5957 //par[2] = sigma
5958 //norm0 = 1
5959 //par[3] = lambda0
5960 //par[4] = norm1
5961 //par[5] = lambda1
5962
5963 Double_t par1save = par[1];
5964 //Double_t par2save = par[2];
5965 Double_t par2save = 0.0604244;
5966 //Double_t par3save = par[3];
5967 Double_t par3save = 7.65596;
5968 //Double_t par5save = par[5];
5969 Double_t par5save = 0.870597;
5970 Double_t dx = x[0] - par1save;
5971
5972 Double_t sigma2 = par2save*par2save;
5973 Double_t sqrt2 = TMath::Sqrt(2.0);
5974 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5975 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5976 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5977 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5978
5979 //return par[0]*(exp1+par[4]*exp2);
5980 return par[0] * (exp1 + 1.00124 * exp2);
5981
5982}
5983
5984//_____________________________________________________________________________
5985Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5986{
5987 //
5988 // Sum Landau + Gaus with identical mean
5989 //
5990
5991 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5992 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5993 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5994 Double_t val = valLandau + valGaus;
5995
5996 return val;
5997
5998}
5999
6000//_____________________________________________________________________________
6001Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6002{
6003 //
6004 // Function for the fit
6005 //
6006 // Fit parameters:
6007 // par[0]=Width (scale) parameter of Landau density
6008 // par[1]=Most Probable (MP, location) parameter of Landau density
6009 // par[2]=Total area (integral -inf to inf, normalization constant)
6010 // par[3]=Width (sigma) of convoluted Gaussian function
6011 //
6012 // In the Landau distribution (represented by the CERNLIB approximation),
6013 // the maximum is located at x=-0.22278298 with the location parameter=0.
6014 // This shift is corrected within this function, so that the actual
6015 // maximum is identical to the MP parameter.
6016 //
6017
6018 // Numeric constants
6019 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6020 Double_t mpshift = -0.22278298; // Landau maximum location
6021
6022 // Control constants
6023 Double_t np = 100.0; // Number of convolution steps
6024 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6025
6026 // Variables
6027 Double_t xx;
6028 Double_t mpc;
6029 Double_t fland;
6030 Double_t sum = 0.0;
6031 Double_t xlow;
6032 Double_t xupp;
6033 Double_t step;
6034 Double_t i;
6035
6036 // MP shift correction
6037 mpc = par[1] - mpshift * par[0];
6038
6039 // Range of convolution integral
6040 xlow = x[0] - sc * par[3];
6041 xupp = x[0] + sc * par[3];
6042
6043 step = (xupp - xlow) / np;
6044
6045 // Convolution integral of Landau and Gaussian by sum
6046 for (i = 1.0; i <= np/2; i++) {
6047
6048 xx = xlow + (i-.5) * step;
6049 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6050 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6051
6052 xx = xupp - (i-.5) * step;
6053 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6054 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6055
6056 }
6057
6058 return (par[2] * step * sum * invsq2pi / par[3]);
6059
6060}
6061//_____________________________________________________________________________
6062TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6063 , const Double_t *parlimitslo, const Double_t *parlimitshi
6064 , Double_t *fitparams, Double_t *fiterrors
6065 , Double_t *chiSqr, Int_t *ndf) const
6066{
6067 //
6068 // Function for the fit
6069 //
6070
6071 Int_t i;
6072 Char_t funname[100];
6073
6074 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6075 if (ffitold) {
6076 delete ffitold;
6077 }
6078
6079 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6080 ffit->SetParameters(startvalues);
6081 ffit->SetParNames("Width","MP","Area","GSigma");
6082
6083 for (i = 0; i < 4; i++) {
6084 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6085 }
6086
6087 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6088
6089 ffit->GetParameters(fitparams); // Obtain fit parameters
6090 for (i = 0; i < 4; i++) {
6091 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6092 }
6093 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6094 ndf[0] = ffit->GetNDF(); // Obtain ndf
6095
6096 return (ffit); // Return fit function
6097
6098}
6099
6100//_____________________________________________________________________________
6101Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6102{
6103 //
6104 // Function for the fit
6105 //
6106
6107 Double_t p;
6108 Double_t x;
6109 Double_t fy;
6110 Double_t fxr;
6111 Double_t fxl;
6112 Double_t step;
6113 Double_t l;
6114 Double_t lold;
6115
6116 Int_t i = 0;
6117 Int_t maxcalls = 10000;
6118
6119 // Search for maximum
6120 p = params[1] - 0.1 * params[0];
6121 step = 0.05 * params[0];
6122 lold = -2.0;
6123 l = -1.0;
6124
6125 while ((l != lold) && (i < maxcalls)) {
6126 i++;
6127 lold = l;
6128 x = p + step;
6129 l = LanGauFun(&x,params);
6130 if (l < lold) {
6131 step = -step / 10.0;
6132 }
6133 p += step;
6134 }
6135
6136 if (i == maxcalls) {
6137 return (-1);
6138 }
6139 maxx = x;
6140 fy = l / 2.0;
6141
6142 // Search for right x location of fy
6143 p = maxx + params[0];
6144 step = params[0];
6145 lold = -2.0;
6146 l = -1e300;
6147 i = 0;
6148
6149 while ( (l != lold) && (i < maxcalls) ) {
6150 i++;
6151
6152 lold = l;
6153 x = p + step;
6154 l = TMath::Abs(LanGauFun(&x,params) - fy);
6155
6156 if (l > lold)
6157 step = -step/10;
6158
6159 p += step;
6160 }
6161
6162 if (i == maxcalls)
6163 return (-2);
6164
6165 fxr = x;
6166
6167
6168 // Search for left x location of fy
6169
6170 p = maxx - 0.5 * params[0];
6171 step = -params[0];
6172 lold = -2.0;
6173 l = -1.0e300;
6174 i = 0;
6175
6176 while ((l != lold) && (i < maxcalls)) {
6177 i++;
6178 lold = l;
6179 x = p + step;
6180 l = TMath::Abs(LanGauFun(&x,params) - fy);
6181 if (l > lold) {
6182 step = -step / 10.0;
6183 }
6184 p += step;
6185 }
6186
6187 if (i == maxcalls) {
6188 return (-3);
6189 }
6190
6191 fxl = x;
6192 fwhm = fxr - fxl;
6193
6194 return (0);
6195}
6196//_____________________________________________________________________________
6197Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6198{
6199 //
6200 // Gaus with identical mean
6201 //
6202
6203 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
6204
6205 return gauss;
6206
6207}
6208
6209
6210