1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
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
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
62 #include <TTreeStream.h>
63 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
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"
85 ClassImp(AliTRDCalibraFit)
87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
90 //_____________singleton implementation_________________________________________________
91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
94 // Singleton implementation
97 if (fgTerminated != kFALSE) {
101 if (fgInstance == 0) {
102 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
124 //______________________________________________________________________________________
125 AliTRDCalibraFit::AliTRDCalibraFit()
128 ,fNumberOfBinsExpected(0)
130 ,fBeginFitCharge(3.5)
132 ,fTakeTheMaxPH(kTRUE)
141 ,fNumberFitSuccess(0)
148 ,fCalibraMode(new AliTRDCalibraMode())
153 ,fScaleFitFactor(0.0)
162 ,fCalDetVdriftUsed(0x0)
164 ,fCurrentCoefDetector(0x0)
165 ,fCurrentCoefDetector2(0x0)
170 // Default constructor
173 fGeo = new AliTRDgeometry();
175 // Current variables initialised
176 for (Int_t k = 0; k < 2; k++) {
177 fCurrentCoef[k] = 0.0;
178 fCurrentCoef2[k] = 0.0;
180 for (Int_t i = 0; i < 3; i++) {
186 //______________________________________________________________________________________
187 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
190 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
192 ,fBeginFitCharge(c.fBeginFitCharge)
193 ,fFitPHPeriode(c.fFitPHPeriode)
194 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
195 ,fT0Shift0(c.fT0Shift0)
196 ,fT0Shift1(c.fT0Shift1)
197 ,fRangeFitPRF(c.fRangeFitPRF)
199 ,fMinEntries(c.fMinEntries)
201 ,fScaleGain(c.fScaleGain)
202 ,fNumberFit(c.fNumberFit)
203 ,fNumberFitSuccess(c.fNumberFitSuccess)
204 ,fNumberEnt(c.fNumberEnt)
205 ,fStatisticMean(c.fStatisticMean)
207 ,fDebugLevel(c.fDebugLevel)
208 ,fFitVoir(c.fFitVoir)
209 ,fMagneticField(c.fMagneticField)
211 ,fCurrentCoefE(c.fCurrentCoefE)
212 ,fCurrentCoefE2(c.fCurrentCoefE2)
215 ,fScaleFitFactor(c.fScaleFitFactor)
216 ,fEntriesCurrent(c.fEntriesCurrent)
217 ,fCountDet(c.fCountDet)
224 ,fCalDetVdriftUsed(0x0)
226 ,fCurrentCoefDetector(0x0)
227 ,fCurrentCoefDetector2(0x0)
235 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
237 //Current variables initialised
238 for (Int_t k = 0; k < 2; k++) {
239 fCurrentCoef[k] = 0.0;
240 fCurrentCoef2[k] = 0.0;
242 for (Int_t i = 0; i < 3; i++) {
246 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
247 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
249 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
250 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
252 if(c.fCalDetVdriftUsed) fCalDetVdriftUsed = new AliTRDCalDet(*c.fCalDetVdriftUsed);
253 if(c.fCalDetExBUsed) fCalDetExBUsed = new AliTRDCalDet(*c.fCalDetExBUsed);
255 fVectorFit.SetName(c.fVectorFit.GetName());
256 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
257 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
258 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
260 if (GetStack(detector) == 2) {
266 Float_t *coef = new Float_t[ntotal];
267 for (Int_t i = 0; i < ntotal; i++) {
268 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
270 fitInfo->SetCoef(coef);
271 fitInfo->SetDetector(detector);
272 fVectorFit.Add((TObject *) fitInfo);
274 fVectorFit.SetName(c.fVectorFit.GetName());
275 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
276 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
277 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
279 if (GetStack(detector) == 2) {
285 Float_t *coef = new Float_t[ntotal];
286 for (Int_t i = 0; i < ntotal; i++) {
287 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
289 fitInfo->SetCoef(coef);
290 fitInfo->SetDetector(detector);
291 fVectorFit2.Add((TObject *) fitInfo);
296 fGeo = new AliTRDgeometry();
299 //____________________________________________________________________________________
300 AliTRDCalibraFit::~AliTRDCalibraFit()
303 // AliTRDCalibraFit destructor
305 if ( fDebugStreamer ) delete fDebugStreamer;
306 if ( fCalDet ) delete fCalDet;
307 if ( fCalDet2 ) delete fCalDet2;
308 if ( fCalROC ) delete fCalROC;
309 if ( fCalROC2 ) delete fCalROC2;
310 if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed;
311 if ( fCalDetExBUsed) delete fCalDetExBUsed;
312 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
313 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
315 fVectorFit2.Delete();
321 //_____________________________________________________________________________
322 void AliTRDCalibraFit::Destroy()
334 //_____________________________________________________________________________
335 void AliTRDCalibraFit::DestroyDebugStreamer()
338 // Delete DebugStreamer
341 if ( fDebugStreamer ) delete fDebugStreamer;
342 fDebugStreamer = 0x0;
345 //__________________________________________________________________________________
346 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
349 // From the drift velocity and t0
350 // return the position of the peak and maximum negative slope
353 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
354 Double_t widbins = 0.1; // 0.1 mus
356 //peak and maxnegslope in mus
357 Double_t begind = t0*widbins + fT0Shift0;
358 Double_t peakd = t0*widbins + fT0Shift1;
359 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
361 // peak and maxnegslope in timebin
362 begin = TMath::Nint(begind*widbins);
363 peak = TMath::Nint(peakd*widbins);
364 end = TMath::Nint(maxnegslope*widbins);
367 //____________Functions fit Online CH2d________________________________________
368 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
371 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
372 // calibration group normalized the resulted coefficients (to 1 normally)
375 // Set the calibration mode
376 //const char *name = ch->GetTitle();
377 TString name = ch->GetTitle();
378 if(!SetModeCalibration(name,0)) return kFALSE;
380 // Number of Ybins (detectors or groups of pads)
381 Int_t nbins = ch->GetNbinsX();// charge
382 Int_t nybins = ch->GetNbinsY();// groups number
383 if (!InitFit(nybins,0)) {
389 fStatisticMean = 0.0;
391 fNumberFitSuccess = 0;
393 // Init fCountDet and fCount
394 InitfCountDetAndfCount(0);
395 // Beginning of the loop betwwen dect1 and dect2
396 for (Int_t idect = fDect1; idect < fDect2; idect++) {
397 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
398 UpdatefCountDetAndfCount(idect,0);
399 ReconstructFitRowMinRowMax(idect, 0);
401 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
402 projch->SetDirectory(0);
403 // Number of entries for this calibration group
404 Double_t nentries = 0.0;
406 for (Int_t k = 0; k < nbins; k++) {
407 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
408 nentries += ch->GetBinContent(binnb);
409 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
410 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
412 projch->SetEntries(nentries);
413 //printf("The number of entries for the group %d is %f\n",idect,nentries);
418 // Rebin and statistic stuff
420 projch = ReBin((TH1I *) projch);
422 // This detector has not enough statistics or was off
423 if (nentries <= fMinEntries) {
424 NotEnoughStatisticCH(idect);
425 if (fDebugLevel != 1) {
430 // Statistics of the group fitted
431 fStatisticMean += nentries;
436 case 0: FitMeanW((TH1 *) projch, nentries); break;
437 case 1: FitMean((TH1 *) projch, nentries, mean); break;
438 case 2: FitCH((TH1 *) projch, mean); break;
439 case 3: FitBisCH((TH1 *) projch, mean); break;
440 default: return kFALSE;
443 FillInfosFitCH(idect);
445 if (fDebugLevel != 1) {
450 if (fDebugLevel != 1) {
454 if (fNumberFit > 0) {
455 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));
456 fStatisticMean = fStatisticMean / fNumberFit;
459 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
461 delete fDebugStreamer;
462 fDebugStreamer = 0x0;
466 //____________Functions fit Online CH2d________________________________________
467 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
470 // Reconstruct a 1D histo from the vectorCH for each calibration group,
471 // fit the histo, normalized the resulted coefficients (to 1 normally)
474 // Set the calibraMode
475 //const char *name = calvect->GetNameCH();
476 TString name = calvect->GetNameCH();
477 if(!SetModeCalibration(name,0)) return kFALSE;
479 // Number of Xbins (detectors or groups of pads)
480 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
486 fStatisticMean = 0.0;
488 fNumberFitSuccess = 0;
490 // Init fCountDet and fCount
491 InitfCountDetAndfCount(0);
492 // Beginning of the loop between dect1 and dect2
493 for (Int_t idect = fDect1; idect < fDect2; idect++) {
494 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
495 UpdatefCountDetAndfCount(idect,0);
496 ReconstructFitRowMinRowMax(idect,0);
498 Double_t nentries = 0.0;
500 if(!calvect->GetCHEntries(fCountDet)) {
501 NotEnoughStatisticCH(idect);
507 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
508 projch->SetDirectory(0);
509 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
510 nentries += projch->GetBinContent(k+1);
511 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
517 //printf("The number of entries for the group %d is %f\n",idect,nentries);
520 projch = ReBin((TH1F *) projch);
522 // This detector has not enough statistics or was not found in VectorCH
523 if (nentries <= fMinEntries) {
524 NotEnoughStatisticCH(idect);
527 // Statistic of the histos fitted
528 fStatisticMean += nentries;
533 case 0: FitMeanW((TH1 *) projch, nentries); break;
534 case 1: FitMean((TH1 *) projch, nentries, mean); break;
535 case 2: FitCH((TH1 *) projch, mean); break;
536 case 3: FitBisCH((TH1 *) projch, mean); break;
537 default: return kFALSE;
540 FillInfosFitCH(idect);
543 if (fDebugLevel != 1) {
547 if (fNumberFit > 0) {
548 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));
549 fStatisticMean = fStatisticMean / fNumberFit;
552 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
554 delete fDebugStreamer;
555 fDebugStreamer = 0x0;
558 //____________Functions fit Online CH2d________________________________________
559 Double_t AliTRDCalibraFit::AnalyseCHAllTogether(const TH2I *ch)
562 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
563 // calibration group normalized the resulted coefficients (to 1 normally)
565 Int_t nbins = ch->GetNbinsX();// charge
566 Int_t nybins = ch->GetNbinsY();// groups number
568 TH1I *projch = (TH1I *) ch->ProjectionX("projch",1,nybins+1,(Option_t *)"e");
569 projch->SetDirectory(0);
570 // Number of entries for this calibration group
571 Double_t nentries = 0.0;
573 for (Int_t k = 0; k < nbins; k++) {
574 nentries += projch->GetBinContent(k+1);
575 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
576 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
578 projch->SetEntries(nentries);
579 //printf("The number of entries for the group %d is %f\n",idect,nentries);
583 // This detector has not enough statistics or was off
584 if (nentries <= fMinEntries) {
586 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
592 case 0: FitMeanW((TH1 *) projch, nentries); break;
593 case 1: FitMean((TH1 *) projch, nentries, mean); break;
594 case 2: FitCH((TH1 *) projch, mean); break;
595 case 3: FitBisCH((TH1 *) projch, mean); break;
596 default: return -100.0;
598 delete fDebugStreamer;
599 fDebugStreamer = 0x0;
601 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
605 //________________functions fit Online PH2d____________________________________
606 Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph)
609 // Take the 1D profiles (average pulse height), projections of the 2D PH
610 // on the Xaxis, for each calibration group
611 // Reconstruct a drift velocity
612 // A first calibration of T0 is also made using the same method
615 // Number of Xbins (detectors or groups of pads)
616 Int_t nbins = ph->GetNbinsX();// time
617 Int_t nybins = ph->GetNbinsY();// calibration group
620 TH1D *projph = (TH1D *) ph->ProjectionX("projph",1,nybins+1,(Option_t *) "e");
621 projph->SetDirectory(0);
622 // Number of entries for this calibration group
623 Double_t nentries = 0;
624 for(Int_t idect = 0; idect < nybins; idect++){
625 for (Int_t k = 0; k < nbins; k++) {
626 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
627 nentries += ph->GetBinEntries(binnb);
630 //printf("AnalysePHAllTogether:: the number of entries is %f\n",nentries);
631 // This detector has not enough statistics or was off
632 if (nentries <= fMinEntries) {
633 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
634 if (fDebugLevel != 1) {
640 //printf("Method\n");
643 case 0: FitLagrangePoly((TH1 *) projph); break;
644 case 1: FitPente((TH1 *) projph); break;
645 case 2: FitPH((TH1 *) projph,0); break;
646 default: return -100.0;
649 if (fDebugLevel != 1) {
652 delete fDebugStreamer;
653 fDebugStreamer = 0x0;
655 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
659 //____________Functions fit Online PH2d________________________________________
660 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
663 // Reconstruct the average pulse height from the vectorPH for each
665 // Reconstruct a drift velocity
666 // A first calibration of T0 is also made using the same method (slope method)
669 // Set the calibration mode
670 //const char *name = calvect->GetNamePH();
671 TString name = calvect->GetNamePH();
672 if(!SetModeCalibration(name,1)) return kFALSE;
674 // Number of Xbins (detectors or groups of pads)
675 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
681 fStatisticMean = 0.0;
683 fNumberFitSuccess = 0;
685 // Init fCountDet and fCount
686 InitfCountDetAndfCount(1);
687 // Beginning of the loop
688 for (Int_t idect = fDect1; idect < fDect2; idect++) {
689 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
690 UpdatefCountDetAndfCount(idect,1);
691 ReconstructFitRowMinRowMax(idect,1);
694 if(!calvect->GetPHEntries(fCountDet)) {
695 NotEnoughStatisticPH(idect,fEntriesCurrent);
700 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
701 projph->SetDirectory(0);
702 if(fEntriesCurrent > 0) fNumberEnt++;
703 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
704 // This detector has not enough statistics or was off
705 if (fEntriesCurrent <= fMinEntries) {
706 //printf("Not enough stat!\n");
707 NotEnoughStatisticPH(idect,fEntriesCurrent);
710 // Statistic of the histos fitted
712 fStatisticMean += fEntriesCurrent;
713 // Calcul of "real" coef
714 CalculVdriftCoefMean();
719 case 0: FitLagrangePoly((TH1 *) projph); break;
720 case 1: FitPente((TH1 *) projph); break;
721 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
722 default: return kFALSE;
724 // Fill the tree if end of a detector or only the pointer to the branch!!!
725 FillInfosFitPH(idect,fEntriesCurrent);
729 if (fNumberFit > 0) {
730 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));
731 fStatisticMean = fStatisticMean / fNumberFit;
734 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
736 delete fDebugStreamer;
737 fDebugStreamer = 0x0;
740 //________________functions fit Online PH2d____________________________________
741 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
744 // Take the 1D profiles (average pulse height), projections of the 2D PH
745 // on the Xaxis, for each calibration group
746 // Reconstruct a drift velocity
747 // A first calibration of T0 is also made using the same method
750 // Set the calibration mode
751 //const char *name = ph->GetTitle();
752 TString name = ph->GetTitle();
753 if(!SetModeCalibration(name,1)) return kFALSE;
755 //printf("Mode calibration set\n");
757 // Number of Xbins (detectors or groups of pads)
758 Int_t nbins = ph->GetNbinsX();// time
759 Int_t nybins = ph->GetNbinsY();// calibration group
760 if (!InitFit(nybins,1)) {
764 //printf("Init fit\n");
770 //printf("Init fit PH\n");
772 fStatisticMean = 0.0;
774 fNumberFitSuccess = 0;
776 // Init fCountDet and fCount
777 InitfCountDetAndfCount(1);
778 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
780 // Beginning of the loop
781 for (Int_t idect = fDect1; idect < fDect2; idect++) {
782 //printf("idect = %d\n",idect);
783 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
784 UpdatefCountDetAndfCount(idect,1);
785 ReconstructFitRowMinRowMax(idect,1);
787 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
788 projph->SetDirectory(0);
789 // Number of entries for this calibration group
790 Double_t nentries = 0;
791 for (Int_t k = 0; k < nbins; k++) {
792 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
793 nentries += ph->GetBinEntries(binnb);
798 //printf("The number of entries for the group %d is %f\n",idect,nentries);
799 // This detector has not enough statistics or was off
800 if (nentries <= fMinEntries) {
801 //printf("Not enough statistic!\n");
802 NotEnoughStatisticPH(idect,nentries);
803 if (fDebugLevel != 1) {
808 // Statistics of the histos fitted
810 fStatisticMean += nentries;
811 // Calcul of "real" coef
812 CalculVdriftCoefMean();
815 //printf("Method\n");
818 case 0: FitLagrangePoly((TH1 *) projph); break;
819 case 1: FitPente((TH1 *) projph); break;
820 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
821 default: return kFALSE;
823 // Fill the tree if end of a detector or only the pointer to the branch!!!
824 FillInfosFitPH(idect,nentries);
826 if (fDebugLevel != 1) {
831 if (fNumberFit > 0) {
832 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));
833 fStatisticMean = fStatisticMean / fNumberFit;
836 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
838 delete fDebugStreamer;
839 fDebugStreamer = 0x0;
842 //____________Functions fit Online PRF2d_______________________________________
843 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
846 // Take the 1D profiles (pad response function), projections of the 2D PRF
847 // on the Xaxis, for each calibration group
848 // Fit with a gaussian to reconstruct the sigma of the pad response function
851 // Set the calibration mode
852 //const char *name = prf->GetTitle();
853 TString name = prf->GetTitle();
854 if(!SetModeCalibration(name,2)) return kFALSE;
856 // Number of Ybins (detectors or groups of pads)
857 Int_t nybins = prf->GetNbinsY();// calibration groups
858 Int_t nbins = prf->GetNbinsX();// bins
859 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
860 if((nbg > 0) || (nbg == -1)) return kFALSE;
861 if (!InitFit(nybins,2)) {
867 fStatisticMean = 0.0;
869 fNumberFitSuccess = 0;
871 // Init fCountDet and fCount
872 InitfCountDetAndfCount(2);
873 // Beginning of the loop
874 for (Int_t idect = fDect1; idect < fDect2; idect++) {
875 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
876 UpdatefCountDetAndfCount(idect,2);
877 ReconstructFitRowMinRowMax(idect,2);
879 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
880 projprf->SetDirectory(0);
881 // Number of entries for this calibration group
882 Double_t nentries = 0;
883 for (Int_t k = 0; k < nbins; k++) {
884 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
885 nentries += prf->GetBinEntries(binnb);
887 if(nentries > 0) fNumberEnt++;
888 // This detector has not enough statistics or was off
889 if (nentries <= fMinEntries) {
890 NotEnoughStatisticPRF(idect);
891 if (fDebugLevel != 1) {
896 // Statistics of the histos fitted
898 fStatisticMean += nentries;
899 // Calcul of "real" coef
904 case 0: FitPRF((TH1 *) projprf); break;
905 case 1: RmsPRF((TH1 *) projprf); break;
906 default: return kFALSE;
908 // Fill the tree if end of a detector or only the pointer to the branch!!!
909 FillInfosFitPRF(idect);
911 if (fDebugLevel != 1) {
916 if (fNumberFit > 0) {
917 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
918 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
919 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
920 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
921 fStatisticMean = fStatisticMean / fNumberFit;
924 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
926 delete fDebugStreamer;
927 fDebugStreamer = 0x0;
930 //____________Functions fit Online PRF2d_______________________________________
931 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
934 // Take the 1D profiles (pad response function), projections of the 2D PRF
935 // on the Xaxis, for each calibration group
936 // Fit with a gaussian to reconstruct the sigma of the pad response function
939 // Set the calibration mode
940 //const char *name = prf->GetTitle();
941 TString name = prf->GetTitle();
942 if(!SetModeCalibration(name,2)) return kFALSE;
944 // Number of Ybins (detectors or groups of pads)
945 TAxis *xprf = prf->GetXaxis();
946 TAxis *yprf = prf->GetYaxis();
947 Int_t nybins = yprf->GetNbins();// calibration groups
948 Int_t nbins = xprf->GetNbins();// bins
949 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
950 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
951 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
952 if(nbg == -1) return kFALSE;
953 if(nbg > 0) fMethod = 1;
955 if (!InitFit(nybins,2)) {
961 fStatisticMean = 0.0;
963 fNumberFitSuccess = 0;
965 // Init fCountDet and fCount
966 InitfCountDetAndfCount(2);
967 // Beginning of the loop
968 for (Int_t idect = fDect1; idect < fDect2; idect++) {
969 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
970 UpdatefCountDetAndfCount(idect,2);
971 ReconstructFitRowMinRowMax(idect,2);
972 // Build the array of entries and sum
973 TArrayD arraye = TArrayD(nbins);
974 TArrayD arraym = TArrayD(nbins);
975 TArrayD arrayme = TArrayD(nbins);
976 Double_t nentries = 0;
977 //printf("nbins %d\n",nbins);
978 for (Int_t k = 0; k < nbins; k++) {
979 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
980 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
981 Double_t mean = (Double_t)prf->GetBinContent(binnb);
982 Double_t error = (Double_t)prf->GetBinError(binnb);
983 //printf("for %d we have %f\n",k,entries);
985 arraye.AddAt(entries,k);
986 arraym.AddAt(mean,k);
987 arrayme.AddAt(error,k);
989 if(nentries > 0) fNumberEnt++;
990 //printf("The number of entries for the group %d is %f\n",idect,nentries);
991 // This detector has not enough statistics or was off
992 if (nentries <= fMinEntries) {
993 NotEnoughStatisticPRF(idect);
996 // Statistics of the histos fitted
998 fStatisticMean += nentries;
999 // Calcul of "real" coef
1000 CalculPRFCoefMean();
1004 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
1005 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
1006 default: return kFALSE;
1008 // Fill the tree if end of a detector or only the pointer to the branch!!!
1009 FillInfosFitPRF(idect);
1012 if (fNumberFit > 0) {
1013 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1014 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1015 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1016 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1017 fStatisticMean = fStatisticMean / fNumberFit;
1020 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1022 delete fDebugStreamer;
1023 fDebugStreamer = 0x0;
1026 //____________Functions fit Online PRF2d_______________________________________
1027 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
1030 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1031 // each calibration group
1032 // Fit with a gaussian to reconstruct the sigma of the pad response function
1035 // Set the calibra mode
1036 //const char *name = calvect->GetNamePRF();
1037 TString name = calvect->GetNamePRF();
1038 if(!SetModeCalibration(name,2)) return kFALSE;
1039 //printf("test0 %s\n",name);
1041 // Number of Xbins (detectors or groups of pads)
1042 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1043 //printf("test1\n");
1046 if (!InitFitPRF()) {
1047 ///printf("test2\n");
1050 fStatisticMean = 0.0;
1052 fNumberFitSuccess = 0;
1054 // Init fCountDet and fCount
1055 InitfCountDetAndfCount(2);
1056 // Beginning of the loop
1057 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1058 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
1059 UpdatefCountDetAndfCount(idect,2);
1060 ReconstructFitRowMinRowMax(idect,2);
1062 fEntriesCurrent = 0;
1063 if(!calvect->GetPRFEntries(fCountDet)) {
1064 NotEnoughStatisticPRF(idect);
1067 TString tname("PRF");
1069 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
1070 projprf->SetDirectory(0);
1071 if(fEntriesCurrent > 0) fNumberEnt++;
1072 // This detector has not enough statistics or was off
1073 if (fEntriesCurrent <= fMinEntries) {
1074 NotEnoughStatisticPRF(idect);
1077 // Statistic of the histos fitted
1079 fStatisticMean += fEntriesCurrent;
1080 // Calcul of "real" coef
1081 CalculPRFCoefMean();
1085 case 1: FitPRF((TH1 *) projprf); break;
1086 case 2: RmsPRF((TH1 *) projprf); break;
1087 default: return kFALSE;
1089 // Fill the tree if end of a detector or only the pointer to the branch!!!
1090 FillInfosFitPRF(idect);
1093 if (fNumberFit > 0) {
1094 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));
1097 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1099 delete fDebugStreamer;
1100 fDebugStreamer = 0x0;
1103 //____________Functions fit Online PRF2d_______________________________________
1104 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
1107 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1108 // each calibration group
1109 // Fit with a gaussian to reconstruct the sigma of the pad response function
1112 // Set the calibra mode
1113 //const char *name = calvect->GetNamePRF();
1114 TString name = calvect->GetNamePRF();
1115 if(!SetModeCalibration(name,2)) return kFALSE;
1116 //printf("test0 %s\n",name);
1117 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1118 //printf("test1 %d\n",nbg);
1119 if(nbg == -1) return kFALSE;
1120 if(nbg > 0) fMethod = 1;
1122 // Number of Xbins (detectors or groups of pads)
1123 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1124 //printf("test2\n");
1127 if (!InitFitPRF()) {
1128 //printf("test3\n");
1131 fStatisticMean = 0.0;
1133 fNumberFitSuccess = 0;
1137 Double_t *arrayx = 0;
1138 Double_t *arraye = 0;
1139 Double_t *arraym = 0;
1140 Double_t *arrayme = 0;
1141 Float_t lowedge = 0.0;
1142 Float_t upedge = 0.0;
1143 // Init fCountDet and fCount
1144 InitfCountDetAndfCount(2);
1145 // Beginning of the loop
1146 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1147 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1148 UpdatefCountDetAndfCount(idect,2);
1149 ReconstructFitRowMinRowMax(idect,2);
1151 fEntriesCurrent = 0;
1152 if(!calvect->GetPRFEntries(fCountDet)) {
1153 NotEnoughStatisticPRF(idect);
1156 TString tname("PRF");
1158 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1159 nbins = projprftree->GetN();
1160 arrayx = (Double_t *)projprftree->GetX();
1161 arraye = (Double_t *)projprftree->GetEX();
1162 arraym = (Double_t *)projprftree->GetY();
1163 arrayme = (Double_t *)projprftree->GetEY();
1164 Float_t step = arrayx[1]-arrayx[0];
1165 lowedge = arrayx[0] - step/2.0;
1166 upedge = arrayx[(nbins-1)] + step/2.0;
1167 //printf("nbins est %d\n",nbins);
1168 for(Int_t k = 0; k < nbins; k++){
1169 fEntriesCurrent += (Int_t)arraye[k];
1170 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1171 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1173 if(fEntriesCurrent > 0) fNumberEnt++;
1174 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1175 // This detector has not enough statistics or was off
1176 if (fEntriesCurrent <= fMinEntries) {
1177 NotEnoughStatisticPRF(idect);
1180 // Statistic of the histos fitted
1182 fStatisticMean += fEntriesCurrent;
1183 // Calcul of "real" coef
1184 CalculPRFCoefMean();
1188 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1189 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1190 default: return kFALSE;
1192 // Fill the tree if end of a detector or only the pointer to the branch!!!
1193 FillInfosFitPRF(idect);
1196 if (fNumberFit > 0) {
1197 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));
1200 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1202 delete fDebugStreamer;
1203 fDebugStreamer = 0x0;
1206 //____________Functions fit Online CH2d________________________________________
1207 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1210 // The linear method
1213 fStatisticMean = 0.0;
1215 fNumberFitSuccess = 0;
1217 if(!InitFitLinearFitter()) return kFALSE;
1220 for(Int_t idet = 0; idet < 540; idet++){
1223 //printf("detector number %d\n",idet);
1228 fEntriesCurrent = 0;
1230 Bool_t here = calivdli->GetParam(idet,¶m);
1231 Bool_t heree = calivdli->GetError(idet,&error);
1232 //printf("here %d and heree %d\n",here, heree);
1234 fEntriesCurrent = (Int_t) error[2];
1237 //printf("Number of entries %d\n",fEntriesCurrent);
1238 // Nothing found or not enough statistic
1239 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1240 NotEnoughStatisticLinearFitter();
1247 fStatisticMean += fEntriesCurrent;
1250 if((-(param[1])) <= 0.0) {
1251 NotEnoughStatisticLinearFitter();
1255 // CalculDatabaseVdriftandTan
1256 CalculVdriftLorentzCoef();
1257 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFit detector %d, vdrift %f and %f and exB %f and %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCurrentCoef[1],fCalDetExBUsed->GetValue(idet),fCurrentCoef2[1]);
1260 fNumberFitSuccess ++;
1262 // Put the fCurrentCoef
1263 fCurrentCoef[0] = -param[1];
1264 // here the database must be the one of the reconstruction for the lorentz angle....
1265 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1266 fCurrentCoefE = error[1];
1267 fCurrentCoefE2 = error[0];
1268 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1269 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1273 FillInfosFitLinearFitter();
1278 if (fNumberFit > 0) {
1279 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));
1282 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1284 delete fDebugStreamer;
1285 fDebugStreamer = 0x0;
1289 //____________Functions fit Online CH2d________________________________________
1290 void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall)
1293 // The linear method
1296 // Get the mean vdrift and exb used
1297 Double_t meanvdriftused = 0.0;
1298 Double_t meanexbused = 0.0;
1299 Double_t counterdet = 0.0;
1300 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) {
1301 vdriftoverall = -100.0;
1308 TH2S *linearfitterhisto = 0x0;
1310 for(Int_t idet = 0; idet < 540; idet++){
1312 TH2S * u = calivdli->GetLinearFitterHistoForce(idet);
1313 Double_t detectorentries = u->Integral();
1314 meanvdriftused += fCalDetVdriftUsed->GetValue(idet)*detectorentries;
1315 meanexbused += fCalDetExBUsed->GetValue(idet)*detectorentries;
1316 counterdet += detectorentries;
1318 //printf("detectorentries %f\n",detectorentries);
1320 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether detector %d, vdrift %f and exB %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCalDetExBUsed->GetValue(idet));
1322 if(idet == 0) linearfitterhisto = u;
1323 else linearfitterhisto->Add(u);
1326 if(counterdet > 0.0){
1327 meanvdriftused = meanvdriftused/counterdet;
1328 meanexbused = meanexbused/counterdet;
1331 vdriftoverall = -100.0;
1337 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether MEAN vdrift %f and exB %f\n",meanvdriftused,meanexbused);
1342 TAxis *xaxis = linearfitterhisto->GetXaxis();
1343 TAxis *yaxis = linearfitterhisto->GetYaxis();
1344 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
1346 Double_t integral = linearfitterhisto->Integral();
1347 //printf("Integral is %f\n",integral);
1348 Bool_t securitybreaking = kFALSE;
1349 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
1350 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
1351 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
1352 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
1353 Double_t x = xaxis->GetBinCenter(ibinx+1);
1354 Double_t y = yaxis->GetBinCenter(ibiny+1);
1356 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
1357 if(!securitybreaking){
1358 linearfitter.AddPoint(&x,y);
1363 linearfitter.AddPoint(&x,y);
1373 //printf("AnalyseLinearFittersAllTogether::Find %d entries\n",entries);
1374 //printf("Minstats %d\n",fMinEntries);
1378 // Eval the linear fitter
1379 if(entries > fMinEntries){
1380 TVectorD par = TVectorD(2);
1382 if((linearfitter.EvalRobust(0.8)==0)) {
1383 //printf("Take the param\n");
1384 linearfitter.GetParameters(par);
1387 //printf("Finish\n");
1388 // Put the fCurrentCoef
1389 fCurrentCoef[0] = -par[1];
1390 // here the database must be the one of the reconstruction for the lorentz angle....
1391 if(fCurrentCoef[0] > 0.0) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
1392 else fCurrentCoef2[0] = 100.0;
1397 fCurrentCoef[0] = -100.0;
1398 fCurrentCoef2[0] = 100.0;
1406 fCurrentCoef[0] = -100.0;
1407 fCurrentCoef2[0] = 100.0;
1411 vdriftoverall = fCurrentCoef[0];
1412 exboverall = fCurrentCoef2[0];
1415 delete linearfitterhisto;
1416 delete fDebugStreamer;
1417 fDebugStreamer = 0x0;
1420 //____________Functions for seeing if the pad is really okey___________________
1421 //_____________________________________________________________________________
1422 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1425 // Get numberofgroupsprf
1429 const Char_t *pattern0 = "Ngp0";
1430 const Char_t *pattern1 = "Ngp1";
1431 const Char_t *pattern2 = "Ngp2";
1432 const Char_t *pattern3 = "Ngp3";
1433 const Char_t *pattern4 = "Ngp4";
1434 const Char_t *pattern5 = "Ngp5";
1435 const Char_t *pattern6 = "Ngp6";
1438 if (strstr(nametitle.Data(),pattern0)) {
1441 if (strstr(nametitle.Data(),pattern1)) {
1444 if (strstr(nametitle.Data(),pattern2)) {
1447 if (strstr(nametitle.Data(),pattern3)) {
1450 if (strstr(nametitle.Data(),pattern4)) {
1453 if (strstr(nametitle.Data(),pattern5)) {
1456 if (strstr(nametitle.Data(),pattern6)){
1463 //_____________________________________________________________________________
1464 Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1467 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1468 // corresponding to the given name
1471 if(!SetNzFromTObject(name,i)) return kFALSE;
1472 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1477 //_____________________________________________________________________________
1478 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1481 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1482 // corresponding to the given TObject
1486 const Char_t *patternrphi0 = "Nrphi0";
1487 const Char_t *patternrphi1 = "Nrphi1";
1488 const Char_t *patternrphi2 = "Nrphi2";
1489 const Char_t *patternrphi3 = "Nrphi3";
1490 const Char_t *patternrphi4 = "Nrphi4";
1491 const Char_t *patternrphi5 = "Nrphi5";
1492 const Char_t *patternrphi6 = "Nrphi6";
1495 const Char_t *patternrphi10 = "Nrphi10";
1496 const Char_t *patternrphi100 = "Nrphi100";
1497 const Char_t *patternz10 = "Nz10";
1498 const Char_t *patternz100 = "Nz100";
1501 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1502 fCalibraMode->SetAllTogether(i);
1504 if (fDebugLevel > 1) {
1505 AliInfo(Form("fNbDet %d and 100",fNbDet));
1509 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1510 fCalibraMode->SetPerSuperModule(i);
1512 if (fDebugLevel > 1) {
1513 AliInfo(Form("fNDet %d and 100",fNbDet));
1518 if (strstr(name.Data(),patternrphi0)) {
1519 fCalibraMode->SetNrphi(i ,0);
1520 if (fDebugLevel > 1) {
1521 AliInfo(Form("fNbDet %d and 0",fNbDet));
1525 if (strstr(name.Data(),patternrphi1)) {
1526 fCalibraMode->SetNrphi(i, 1);
1527 if (fDebugLevel > 1) {
1528 AliInfo(Form("fNbDet %d and 1",fNbDet));
1532 if (strstr(name.Data(),patternrphi2)) {
1533 fCalibraMode->SetNrphi(i, 2);
1534 if (fDebugLevel > 1) {
1535 AliInfo(Form("fNbDet %d and 2",fNbDet));
1539 if (strstr(name.Data(),patternrphi3)) {
1540 fCalibraMode->SetNrphi(i, 3);
1541 if (fDebugLevel > 1) {
1542 AliInfo(Form("fNbDet %d and 3",fNbDet));
1546 if (strstr(name.Data(),patternrphi4)) {
1547 fCalibraMode->SetNrphi(i, 4);
1548 if (fDebugLevel > 1) {
1549 AliInfo(Form("fNbDet %d and 4",fNbDet));
1553 if (strstr(name.Data(),patternrphi5)) {
1554 fCalibraMode->SetNrphi(i, 5);
1555 if (fDebugLevel > 1) {
1556 AliInfo(Form("fNbDet %d and 5",fNbDet));
1560 if (strstr(name.Data(),patternrphi6)) {
1561 fCalibraMode->SetNrphi(i, 6);
1562 if (fDebugLevel > 1) {
1563 AliInfo(Form("fNbDet %d and 6",fNbDet));
1568 if (fDebugLevel > 1) {
1569 AliInfo(Form("fNbDet %d and rest",fNbDet));
1571 fCalibraMode->SetNrphi(i ,0);
1575 //_____________________________________________________________________________
1576 Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1579 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1580 // corresponding to the given TObject
1584 const Char_t *patternz0 = "Nz0";
1585 const Char_t *patternz1 = "Nz1";
1586 const Char_t *patternz2 = "Nz2";
1587 const Char_t *patternz3 = "Nz3";
1588 const Char_t *patternz4 = "Nz4";
1590 const Char_t *patternrphi10 = "Nrphi10";
1591 const Char_t *patternrphi100 = "Nrphi100";
1592 const Char_t *patternz10 = "Nz10";
1593 const Char_t *patternz100 = "Nz100";
1595 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1596 fCalibraMode->SetAllTogether(i);
1598 if (fDebugLevel > 1) {
1599 AliInfo(Form("fNbDet %d and 100",fNbDet));
1603 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1604 fCalibraMode->SetPerSuperModule(i);
1606 if (fDebugLevel > 1) {
1607 AliInfo(Form("fNbDet %d and 10",fNbDet));
1611 if (strstr(name.Data(),patternz0)) {
1612 fCalibraMode->SetNz(i, 0);
1613 if (fDebugLevel > 1) {
1614 AliInfo(Form("fNbDet %d and 0",fNbDet));
1618 if (strstr(name.Data(),patternz1)) {
1619 fCalibraMode->SetNz(i ,1);
1620 if (fDebugLevel > 1) {
1621 AliInfo(Form("fNbDet %d and 1",fNbDet));
1625 if (strstr(name.Data(),patternz2)) {
1626 fCalibraMode->SetNz(i ,2);
1627 if (fDebugLevel > 1) {
1628 AliInfo(Form("fNbDet %d and 2",fNbDet));
1632 if (strstr(name.Data(),patternz3)) {
1633 fCalibraMode->SetNz(i ,3);
1634 if (fDebugLevel > 1) {
1635 AliInfo(Form("fNbDet %d and 3",fNbDet));
1639 if (strstr(name.Data(),patternz4)) {
1640 fCalibraMode->SetNz(i ,4);
1641 if (fDebugLevel > 1) {
1642 AliInfo(Form("fNbDet %d and 4",fNbDet));
1647 if (fDebugLevel > 1) {
1648 AliInfo(Form("fNbDet %d and rest",fNbDet));
1650 fCalibraMode->SetNz(i ,0);
1653 //______________________________________________________________________
1654 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1656 // Remove the results too far from the mean value and rms
1657 // type: 0 gain, 1 vdrift
1661 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1663 AliInfo("The Vector Fit is not complete!");
1666 Int_t detector = -1;
1668 Float_t value = 0.0;
1670 /////////////////////////////////
1671 // Calculate the mean values
1672 ////////////////////////////////
1674 ////////////////////////
1675 Double_t meanAll = 0.0;
1676 Double_t rmsAll = 0.0;
1681 for (Int_t k = 0; k < loop; k++) {
1682 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1683 sector = GetSector(detector);
1685 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1687 rmsAll += value*value;
1693 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1694 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1695 for (Int_t row = 0; row < rowMax; row++) {
1696 for (Int_t col = 0; col < colMax; col++) {
1697 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1699 rmsAll += value*value;
1709 meanAll = meanAll/countAll;
1710 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1712 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1713 /////////////////////////////////////////////////
1715 ////////////////////////////////////////////////
1716 Double_t defaultvalue = -1.0;
1717 if(type==1) defaultvalue = -1.5;
1718 for (Int_t k = 0; k < loop; k++) {
1719 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1720 sector = GetSector(detector);
1721 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1722 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1723 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1725 // remove the results too far away
1726 for (Int_t row = 0; row < rowMax; row++) {
1727 for (Int_t col = 0; col < colMax; col++) {
1728 value = coef[(Int_t)(col*rowMax+row)];
1729 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1730 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1736 //______________________________________________________________________
1737 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1739 // Remove the results too far from the mean and rms
1743 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1745 AliInfo("The Vector Fit is not complete!");
1748 Int_t detector = -1;
1750 Float_t value = 0.0;
1752 /////////////////////////////////
1753 // Calculate the mean values
1754 ////////////////////////////////
1756 ////////////////////////
1757 Double_t meanAll = 0.0;
1758 Double_t rmsAll = 0.0;
1763 for (Int_t k = 0; k < loop; k++) {
1764 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1765 sector = GetSector(detector);
1767 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1770 rmsAll += value*value;
1775 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1776 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1777 for (Int_t row = 0; row < rowMax; row++) {
1778 for (Int_t col = 0; col < colMax; col++) {
1779 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1781 rmsAll += value*value;
1790 meanAll = meanAll/countAll;
1791 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1793 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1794 /////////////////////////////////////////////////
1796 ////////////////////////////////////////////////
1797 for (Int_t k = 0; k < loop; k++) {
1798 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1799 sector = GetSector(detector);
1800 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1801 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1802 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1804 // remove the results too far away
1805 for (Int_t row = 0; row < rowMax; row++) {
1806 for (Int_t col = 0; col < colMax; col++) {
1807 value = coef[(Int_t)(col*rowMax+row)];
1808 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) {
1809 //printf("value outlier %f\n",value);
1810 coef[(Int_t)(col*rowMax+row)] = 100.0;
1816 //______________________________________________________________________
1817 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1819 // ofwhat is equaled to 0: mean value of all passing detectors
1820 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1823 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1825 AliInfo("The Vector Fit is not complete!");
1828 Int_t detector = -1;
1830 Float_t value = 0.0;
1832 /////////////////////////////////
1833 // Calculate the mean values
1834 ////////////////////////////////
1836 ////////////////////////
1837 Double_t meanAll = 0.0;
1838 Double_t meanSupermodule[18];
1839 Double_t meanDetector[540];
1840 Double_t rmsAll = 0.0;
1841 Double_t rmsSupermodule[18];
1842 Double_t rmsDetector[540];
1844 Int_t countSupermodule[18];
1845 Int_t countDetector[540];
1846 for(Int_t sm = 0; sm < 18; sm++){
1847 rmsSupermodule[sm] = 0.0;
1848 meanSupermodule[sm] = 0.0;
1849 countSupermodule[sm] = 0;
1851 for(Int_t det = 0; det < 540; det++){
1852 rmsDetector[det] = 0.0;
1853 meanDetector[det] = 0.0;
1854 countDetector[det] = 0;
1859 for (Int_t k = 0; k < loop; k++) {
1860 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1861 sector = GetSector(detector);
1863 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1865 rmsDetector[detector] += value*value;
1866 meanDetector[detector] += value;
1867 countDetector[detector]++;
1868 rmsSupermodule[sector] += value*value;
1869 meanSupermodule[sector] += value;
1870 countSupermodule[sector]++;
1871 rmsAll += value*value;
1877 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1878 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1879 for (Int_t row = 0; row < rowMax; row++) {
1880 for (Int_t col = 0; col < colMax; col++) {
1881 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1883 rmsDetector[detector] += value*value;
1884 meanDetector[detector] += value;
1885 countDetector[detector]++;
1886 rmsSupermodule[sector] += value*value;
1887 meanSupermodule[sector] += value;
1888 countSupermodule[sector]++;
1889 rmsAll += value*value;
1899 meanAll = meanAll/countAll;
1900 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1902 for(Int_t sm = 0; sm < 18; sm++){
1903 if(countSupermodule[sm] > 0) {
1904 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1905 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1908 for(Int_t det = 0; det < 540; det++){
1909 if(countDetector[det] > 0) {
1910 meanDetector[det] = meanDetector[det]/countDetector[det];
1911 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1914 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1915 ///////////////////////////////////////////////
1916 // Put the mean value for the no-fitted
1917 /////////////////////////////////////////////
1918 for (Int_t k = 0; k < loop; k++) {
1919 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1920 sector = GetSector(detector);
1921 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1922 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1923 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1925 for (Int_t row = 0; row < rowMax; row++) {
1926 for (Int_t col = 0; col < colMax; col++) {
1927 value = coef[(Int_t)(col*rowMax+row)];
1929 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1931 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1932 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1933 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1937 if(fDebugLevel > 1){
1939 if ( !fDebugStreamer ) {
1941 TDirectory *backup = gDirectory;
1942 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1943 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1946 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1948 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1949 "detector="<<detector<<
1961 //______________________________________________________________________
1962 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1964 // ofwhat is equaled to 0: mean value of all passing detectors
1965 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1968 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1970 AliInfo("The Vector Fit is not complete!");
1973 Int_t detector = -1;
1975 Float_t value = 0.0;
1977 /////////////////////////////////
1978 // Calculate the mean values
1979 ////////////////////////////////
1981 ////////////////////////
1982 Double_t meanAll = 0.0;
1983 Double_t rmsAll = 0.0;
1984 Double_t meanSupermodule[18];
1985 Double_t rmsSupermodule[18];
1986 Double_t meanDetector[540];
1987 Double_t rmsDetector[540];
1989 Int_t countSupermodule[18];
1990 Int_t countDetector[540];
1991 for(Int_t sm = 0; sm < 18; sm++){
1992 rmsSupermodule[sm] = 0.0;
1993 meanSupermodule[sm] = 0.0;
1994 countSupermodule[sm] = 0;
1996 for(Int_t det = 0; det < 540; det++){
1997 rmsDetector[det] = 0.0;
1998 meanDetector[det] = 0.0;
1999 countDetector[det] = 0;
2003 for (Int_t k = 0; k < loop; k++) {
2004 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2005 sector = GetSector(detector);
2007 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
2009 rmsDetector[detector] += value*value;
2010 meanDetector[detector] += value;
2011 countDetector[detector]++;
2012 rmsSupermodule[sector] += value*value;
2013 meanSupermodule[sector] += value;
2014 countSupermodule[sector]++;
2016 rmsAll += value*value;
2021 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2022 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2023 for (Int_t row = 0; row < rowMax; row++) {
2024 for (Int_t col = 0; col < colMax; col++) {
2025 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2027 rmsDetector[detector] += value*value;
2028 meanDetector[detector] += value;
2029 countDetector[detector]++;
2030 rmsSupermodule[sector] += value*value;
2031 meanSupermodule[sector] += value;
2032 countSupermodule[sector]++;
2033 rmsAll += value*value;
2043 meanAll = meanAll/countAll;
2044 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
2046 for(Int_t sm = 0; sm < 18; sm++){
2047 if(countSupermodule[sm] > 0) {
2048 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
2049 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
2052 for(Int_t det = 0; det < 540; det++){
2053 if(countDetector[det] > 0) {
2054 meanDetector[det] = meanDetector[det]/countDetector[det];
2055 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
2058 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
2059 ////////////////////////////////////////////
2060 // Put the mean value for the no-fitted
2061 /////////////////////////////////////////////
2062 for (Int_t k = 0; k < loop; k++) {
2063 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2064 sector = GetSector(detector);
2065 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2066 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2067 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
2069 for (Int_t row = 0; row < rowMax; row++) {
2070 for (Int_t col = 0; col < colMax; col++) {
2071 value = coef[(Int_t)(col*rowMax+row)];
2073 if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2075 if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
2076 else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
2077 else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2081 if(fDebugLevel > 1){
2083 if ( !fDebugStreamer ) {
2085 TDirectory *backup = gDirectory;
2086 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2087 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2090 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2092 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
2093 "detector="<<detector<<
2106 //_____________________________________________________________________________
2107 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
2110 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2111 // It takes the mean value of the coefficients per detector
2112 // This object has to be written in the database
2115 // Create the DetObject
2116 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2118 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2119 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2120 Int_t detector = -1;
2121 Float_t value = 0.0;
2124 for (Int_t k = 0; k < loop; k++) {
2125 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2128 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
2132 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2133 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2134 for (Int_t row = 0; row < rowMax; row++) {
2135 for (Int_t col = 0; col < colMax; col++) {
2136 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2137 mean += TMath::Abs(value);
2141 if(count > 0) mean = mean/count;
2143 object->SetValue(detector,mean);
2148 //_____________________________________________________________________________
2149 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
2152 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2153 // It takes the mean value of the coefficients per detector
2154 // This object has to be written in the database
2157 // Create the DetObject
2158 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2160 fScaleGain = scaleFitFactor;
2162 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2163 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2164 Int_t detector = -1;
2165 Float_t value = 0.0;
2167 for (Int_t k = 0; k < loop; k++) {
2168 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2171 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2172 if(!meanOtherBefore){
2173 if(value > 0) value = value*scaleFitFactor;
2175 else value = value*scaleFitFactor;
2176 mean = TMath::Abs(value);
2180 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2181 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2182 for (Int_t row = 0; row < rowMax; row++) {
2183 for (Int_t col = 0; col < colMax; col++) {
2184 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2185 if(!meanOtherBefore) {
2186 if(value > 0) value = value*scaleFitFactor;
2188 else value = value*scaleFitFactor;
2189 mean += TMath::Abs(value);
2193 if(count > 0) mean = mean/count;
2195 if(mean < 0.1) mean = 0.1;
2196 object->SetValue(detector,mean);
2201 //_____________________________________________________________________________
2202 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
2205 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2206 // It takes the min value of the coefficients per detector
2207 // This object has to be written in the database
2210 // Create the DetObject
2211 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2213 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2214 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2215 Int_t detector = -1;
2216 Float_t value = 0.0;
2218 for (Int_t k = 0; k < loop; k++) {
2219 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2220 Float_t min = 100.0;
2222 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2223 //printf("Create det object %f for %d\n",value,k);
2225 if(value > 70.0) value = value-100.0;
2230 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2231 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2232 for (Int_t row = 0; row < rowMax; row++) {
2233 for (Int_t col = 0; col < colMax; col++) {
2234 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2236 if(value > 70.0) value = value-100.0;
2238 if(min > value) min = value;
2242 object->SetValue(detector,min);
2248 //_____________________________________________________________________________
2249 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2252 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2253 // It takes the min value of the coefficients per detector
2254 // This object has to be written in the database
2257 // Create the DetObject
2258 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2261 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2262 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2263 Int_t detector = -1;
2264 Float_t value = 0.0;
2266 for (Int_t k = 0; k < loop; k++) {
2267 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2269 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2270 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2271 Float_t min = 100.0;
2272 for (Int_t row = 0; row < rowMax; row++) {
2273 for (Int_t col = 0; col < colMax; col++) {
2274 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2275 mean += -TMath::Abs(value);
2279 if(count > 0) mean = mean/count;
2281 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2282 if(value > 70.0) value = value-100.0;
2283 object->SetValue(detector,value);
2289 //_____________________________________________________________________________
2290 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2293 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2294 // You need first to create the object for the detectors,
2295 // where the mean value is put.
2296 // This object has to be written in the database
2299 // Create the DetObject
2300 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2303 for(Int_t k = 0; k < 540; k++){
2304 AliTRDCalROC *calROC = object->GetCalROC(k);
2305 Int_t nchannels = calROC->GetNchannels();
2306 for(Int_t ch = 0; ch < nchannels; ch++){
2307 calROC->SetValue(ch,1.0);
2313 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2314 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2315 Int_t detector = -1;
2316 Float_t value = 0.0;
2318 for (Int_t k = 0; k < loop; k++) {
2319 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2320 AliTRDCalROC *calROC = object->GetCalROC(detector);
2321 Float_t mean = detobject->GetValue(detector);
2322 if(TMath::Abs(mean) <= 0.0000000001) continue;
2323 Int_t rowMax = calROC->GetNrows();
2324 Int_t colMax = calROC->GetNcols();
2325 for (Int_t row = 0; row < rowMax; row++) {
2326 for (Int_t col = 0; col < colMax; col++) {
2327 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2328 if(value > 0) value = value*scaleFitFactor;
2329 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2337 //_____________________________________________________________________________
2338 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2341 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2342 // You need first to create the object for the detectors,
2343 // where the mean value is put.
2344 // This object has to be written in the database
2347 // Create the DetObject
2348 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2351 for(Int_t k = 0; k < 540; k++){
2352 AliTRDCalROC *calROC = object->GetCalROC(k);
2353 Int_t nchannels = calROC->GetNchannels();
2354 for(Int_t ch = 0; ch < nchannels; ch++){
2355 calROC->SetValue(ch,1.0);
2361 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2362 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2363 Int_t detector = -1;
2364 Float_t value = 0.0;
2366 for (Int_t k = 0; k < loop; k++) {
2367 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2368 AliTRDCalROC *calROC = object->GetCalROC(detector);
2369 Float_t mean = detobject->GetValue(detector);
2370 if(mean == 0) continue;
2371 Int_t rowMax = calROC->GetNrows();
2372 Int_t colMax = calROC->GetNcols();
2373 for (Int_t row = 0; row < rowMax; row++) {
2374 for (Int_t col = 0; col < colMax; col++) {
2375 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2376 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2384 //_____________________________________________________________________________
2385 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2388 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2389 // You need first to create the object for the detectors,
2390 // where the mean value is put.
2391 // This object has to be written in the database
2394 // Create the DetObject
2395 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2398 for(Int_t k = 0; k < 540; k++){
2399 AliTRDCalROC *calROC = object->GetCalROC(k);
2400 Int_t nchannels = calROC->GetNchannels();
2401 for(Int_t ch = 0; ch < nchannels; ch++){
2402 calROC->SetValue(ch,0.0);
2408 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2409 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2410 Int_t detector = -1;
2411 Float_t value = 0.0;
2413 for (Int_t k = 0; k < loop; k++) {
2414 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2415 AliTRDCalROC *calROC = object->GetCalROC(detector);
2416 Float_t min = detobject->GetValue(detector);
2417 Int_t rowMax = calROC->GetNrows();
2418 Int_t colMax = calROC->GetNcols();
2419 for (Int_t row = 0; row < rowMax; row++) {
2420 for (Int_t col = 0; col < colMax; col++) {
2421 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2423 if(value > 70.0) value = value - 100.0;
2425 calROC->SetValue(col,row,value-min);
2433 //_____________________________________________________________________________
2434 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2437 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2438 // This object has to be written in the database
2441 // Create the DetObject
2442 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2444 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2445 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2446 Int_t detector = -1;
2447 Float_t value = 0.0;
2449 for (Int_t k = 0; k < loop; k++) {
2450 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2451 AliTRDCalROC *calROC = object->GetCalROC(detector);
2452 Int_t rowMax = calROC->GetNrows();
2453 Int_t colMax = calROC->GetNcols();
2454 for (Int_t row = 0; row < rowMax; row++) {
2455 for (Int_t col = 0; col < colMax; col++) {
2456 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2457 calROC->SetValue(col,row,TMath::Abs(value));
2465 //_____________________________________________________________________________
2466 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2469 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2470 // 0 successful fit 1 not successful fit
2471 // mean is the mean value over the successful fit
2472 // do not use it for t0: no meaning
2475 // Create the CalObject
2476 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2480 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2482 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2483 for(Int_t k = 0; k < 540; k++){
2484 object->SetValue(k,1.0);
2487 Int_t detector = -1;
2488 Float_t value = 0.0;
2490 for (Int_t k = 0; k < loop; k++) {
2491 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2492 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2493 if(value <= 0) object->SetValue(detector,1.0);
2495 object->SetValue(detector,0.0);
2500 if(count > 0) mean /= count;
2503 //_____________________________________________________________________________
2504 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2507 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2508 // 0 not successful fit 1 successful fit
2509 // mean mean value over the successful fit
2512 // Create the CalObject
2513 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2517 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2519 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2520 for(Int_t k = 0; k < 540; k++){
2521 AliTRDCalROC *calROC = object->GetCalROC(k);
2522 Int_t nchannels = calROC->GetNchannels();
2523 for(Int_t ch = 0; ch < nchannels; ch++){
2524 calROC->SetValue(ch,1.0);
2528 Int_t detector = -1;
2529 Float_t value = 0.0;
2531 for (Int_t k = 0; k < loop; k++) {
2532 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2533 AliTRDCalROC *calROC = object->GetCalROC(detector);
2534 Int_t nchannels = calROC->GetNchannels();
2535 for (Int_t ch = 0; ch < nchannels; ch++) {
2536 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2537 if(value <= 0) calROC->SetValue(ch,1.0);
2539 calROC->SetValue(ch,0.0);
2545 if(count > 0) mean /= count;
2548 //_____________________________________________________________________________
2549 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2552 // Set FitPH if 1 then each detector will be fitted
2555 if (periodeFitPH > 0) {
2556 fFitPHPeriode = periodeFitPH;
2559 AliInfo("periodeFitPH must be higher than 0!");
2563 //_____________________________________________________________________________
2564 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2567 // The fit of the deposited charge distribution begins at
2568 // histo->Mean()/beginFitCharge
2569 // You can here set beginFitCharge
2572 if (beginFitCharge > 0) {
2573 fBeginFitCharge = beginFitCharge;
2576 AliInfo("beginFitCharge must be strict positif!");
2581 //_____________________________________________________________________________
2582 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2585 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2586 // You can here set t0Shift0
2590 fT0Shift0 = t0Shift;
2593 AliInfo("t0Shift0 must be strict positif!");
2598 //_____________________________________________________________________________
2599 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2602 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2603 // You can here set t0Shift1
2607 fT0Shift1 = t0Shift;
2610 AliInfo("t0Shift must be strict positif!");
2615 //_____________________________________________________________________________
2616 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2619 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2620 // You can here set rangeFitPRF
2623 if ((rangeFitPRF > 0) &&
2624 (rangeFitPRF <= 1.5)) {
2625 fRangeFitPRF = rangeFitPRF;
2628 AliInfo("rangeFitPRF must be between 0 and 1.0");
2633 //_____________________________________________________________________________
2634 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2637 // Minimum entries for fitting
2640 if (minEntries > 0) {
2641 fMinEntries = minEntries;
2644 AliInfo("fMinEntries must be >= 0.");
2649 //_____________________________________________________________________________
2650 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2653 // Rebin with rebin time less bins the Ch histo
2654 // You can set here rebin that should divide the number of bins of CH histo
2659 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2662 AliInfo("You have to choose a positiv value!");
2666 //_____________________________________________________________________________
2667 Bool_t AliTRDCalibraFit::FillVectorFit()
2670 // For the Fit functions fill the vector Fit
2673 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2676 if (GetStack(fCountDet) == 2) {
2683 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2684 Float_t *coef = new Float_t[ntotal];
2685 for (Int_t i = 0; i < ntotal; i++) {
2686 coef[i] = fCurrentCoefDetector[i];
2689 Int_t detector = fCountDet;
2691 fitInfo->SetCoef(coef);
2692 fitInfo->SetDetector(detector);
2693 fVectorFit.Add((TObject *) fitInfo);
2698 //_____________________________________________________________________________
2699 Bool_t AliTRDCalibraFit::FillVectorFit2()
2702 // For the Fit functions fill the vector Fit
2705 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2708 if (GetStack(fCountDet) == 2) {
2715 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2716 Float_t *coef = new Float_t[ntotal];
2717 for (Int_t i = 0; i < ntotal; i++) {
2718 coef[i] = fCurrentCoefDetector2[i];
2721 Int_t detector = fCountDet;
2723 fitInfo->SetCoef(coef);
2724 fitInfo->SetDetector(detector);
2725 fVectorFit2.Add((TObject *) fitInfo);
2730 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2731 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2734 // Init the number of expected bins and fDect1[i] fDect2[i]
2737 gStyle->SetPalette(1);
2738 gStyle->SetOptStat(1111);
2739 gStyle->SetPadBorderMode(0);
2740 gStyle->SetCanvasColor(10);
2741 gStyle->SetPadLeftMargin(0.13);
2742 gStyle->SetPadRightMargin(0.01);
2744 // Mode groups of pads: the total number of bins!
2745 CalculNumberOfBinsExpected(i);
2747 // Quick verification that we have the good pad calibration mode!
2748 if (fNumberOfBinsExpected != nbins) {
2749 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2753 // Security for fDebug 3 and 4
2754 if ((fDebugLevel >= 3) &&
2758 AliInfo("This detector doesn't exit!");
2762 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2763 CalculDect1Dect2(i);
2768 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2769 Bool_t AliTRDCalibraFit::InitFitCH()
2772 // Init the fVectorFitCH for normalisation
2773 // Init the histo for debugging
2778 fScaleFitFactor = 0.0;
2779 fCurrentCoefDetector = new Float_t[2304];
2780 for (Int_t k = 0; k < 2304; k++) {
2781 fCurrentCoefDetector[k] = 0.0;
2783 fVectorFit.SetName("gainfactorscoefficients");
2785 // fDebug == 0 nothing
2786 // fDebug == 1 and fFitVoir no histo
2787 if (fDebugLevel == 1) {
2788 if(!CheckFitVoir()) return kFALSE;
2790 //Get the CalDet object
2792 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2794 AliInfo("Could not get calibDB");
2797 if(fCalDet) delete fCalDet;
2798 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2801 Float_t devalue = 1.0;
2802 if(fCalDet) delete fCalDet;
2803 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2804 for(Int_t k = 0; k < 540; k++){
2805 fCalDet->SetValue(k,devalue);
2811 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2812 Bool_t AliTRDCalibraFit::InitFitPH()
2815 // Init the arrays of results
2816 // Init the histos for debugging
2820 fVectorFit.SetName("driftvelocitycoefficients");
2821 fVectorFit2.SetName("t0coefficients");
2823 fCurrentCoefDetector = new Float_t[2304];
2824 for (Int_t k = 0; k < 2304; k++) {
2825 fCurrentCoefDetector[k] = 0.0;
2828 fCurrentCoefDetector2 = new Float_t[2304];
2829 for (Int_t k = 0; k < 2304; k++) {
2830 fCurrentCoefDetector2[k] = 0.0;
2833 //fDebug == 0 nothing
2834 // fDebug == 1 and fFitVoir no histo
2835 if (fDebugLevel == 1) {
2836 if(!CheckFitVoir()) return kFALSE;
2838 //Get the CalDet object
2840 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2842 AliInfo("Could not get calibDB");
2845 if(fCalDet) delete fCalDet;
2846 if(fCalDet2) delete fCalDet2;
2847 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2848 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2851 Float_t devalue = 1.5;
2852 Float_t devalue2 = 0.0;
2853 if(fCalDet) delete fCalDet;
2854 if(fCalDet2) delete fCalDet2;
2855 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2856 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2857 for(Int_t k = 0; k < 540; k++){
2858 fCalDet->SetValue(k,devalue);
2859 fCalDet2->SetValue(k,devalue2);
2864 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2865 Bool_t AliTRDCalibraFit::InitFitPRF()
2868 // Init the calibration mode (Nz, Nrphi), the histograms for
2869 // debugging the fit methods if fDebug > 0,
2873 fVectorFit.SetName("prfwidthcoefficients");
2875 fCurrentCoefDetector = new Float_t[2304];
2876 for (Int_t k = 0; k < 2304; k++) {
2877 fCurrentCoefDetector[k] = 0.0;
2880 // fDebug == 0 nothing
2881 // fDebug == 1 and fFitVoir no histo
2882 if (fDebugLevel == 1) {
2883 if(!CheckFitVoir()) return kFALSE;
2887 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2888 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2891 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2896 fCurrentCoefDetector = new Float_t[2304];
2897 fCurrentCoefDetector2 = new Float_t[2304];
2898 for (Int_t k = 0; k < 2304; k++) {
2899 fCurrentCoefDetector[k] = 0.0;
2900 fCurrentCoefDetector2[k] = 0.0;
2903 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) return kFALSE;
2907 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2908 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2911 // Init the current detector where we are fCountDet and the
2912 // next fCount for the functions Fit...
2915 // Loop on the Xbins of ch!!
2916 fCountDet = -1; // Current detector
2917 fCount = 0; // To find the next detector
2920 if (fDebugLevel >= 3) {
2921 // Set countdet to the detector
2922 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2923 // Set counter to write at the end of the detector
2925 // Get the right calib objects
2928 if(fDebugLevel == 1) {
2930 fCalibraMode->CalculXBins(fCountDet,i);
2931 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2932 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2934 fCalibraMode->CalculXBins(fCountDet,i);
2935 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2941 fCount = fCalibraMode->GetXbins(i);
2943 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2944 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2945 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2946 ,(Int_t) GetStack(fCountDet)
2947 ,(Int_t) GetSector(fCountDet),i);
2950 //_______________________________________________________________________________
2951 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2954 // Calculate the number of bins expected (calibration groups)
2957 fNumberOfBinsExpected = 0;
2959 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2960 fNumberOfBinsExpected = 1;
2964 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2965 fNumberOfBinsExpected = 18;
2969 fCalibraMode->ModePadCalibration(2,i);
2970 fCalibraMode->ModePadFragmentation(0,2,0,i);
2971 fCalibraMode->SetDetChamb2(i);
2972 if (fDebugLevel > 1) {
2973 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2975 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2976 fCalibraMode->ModePadCalibration(0,i);
2977 fCalibraMode->ModePadFragmentation(0,0,0,i);
2978 fCalibraMode->SetDetChamb0(i);
2979 if (fDebugLevel > 1) {
2980 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2982 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2985 //_______________________________________________________________________________
2986 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2989 // Calculate the range of fits
2994 if (fDebugLevel == 1) {
2998 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
3000 fDect2 = fNumberOfBinsExpected;
3002 if (fDebugLevel >= 3) {
3003 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3004 fCalibraMode->CalculXBins(fCountDet,i);
3005 fDect1 = fCalibraMode->GetXbins(i);
3006 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3007 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3008 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3009 ,(Int_t) GetStack(fCountDet)
3010 ,(Int_t) GetSector(fCountDet),i);
3011 // Set for the next detector
3012 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3015 //_______________________________________________________________________________
3016 Bool_t AliTRDCalibraFit::CheckFitVoir()
3019 // Check if fFitVoir is in the range
3022 if (fFitVoir < fNumberOfBinsExpected) {
3023 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
3026 AliInfo("fFitVoir is out of range of the histo!");
3031 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3032 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
3035 // See if we are in a new detector and update the
3036 // variables fNfragZ and fNfragRphi if yes
3037 // Will never happen for only one detector (3 and 4)
3038 // Doesn't matter for 2
3040 if (fCount == idect) {
3041 // On en est au detector (or first detector in the group)
3043 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
3044 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3045 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3046 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3047 ,(Int_t) GetStack(fCountDet)
3048 ,(Int_t) GetSector(fCountDet),i);
3049 // Set for the next detector
3050 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3055 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3056 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
3059 // Reconstruct the min pad row, max pad row, min pad col and
3060 // max pad col of the calibration group for the Fit functions
3061 // idect is the calibration group inside the detector
3063 if (fDebugLevel != 1) {
3064 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
3066 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
3067 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
3069 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3070 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
3073 // For the case where there are not enough entries in the histograms
3074 // of the calibration group, the value present in the choosen database
3075 // will be put. A negativ sign enables to know that a fit was not possible.
3078 if (fDebugLevel == 1) {
3079 AliInfo("The element has not enough statistic to be fitted");
3081 else if (fNbDet > 0){
3082 Int_t firstdetector = fCountDet;
3083 Int_t lastdetector = fCountDet+fNbDet;
3084 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3085 // loop over detectors
3086 for(Int_t det = firstdetector; det < lastdetector; det++){
3088 //Set the calibration object again
3092 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3094 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3095 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3096 ,(Int_t) GetStack(fCountDet)
3097 ,(Int_t) GetSector(fCountDet),0);
3098 // Reconstruct row min row max
3099 ReconstructFitRowMinRowMax(idect,0);
3101 // Calcul the coef from the database choosen for the detector
3102 CalculChargeCoefMean(kFALSE);
3104 //stack 2, not stack 2
3106 if(GetStack(fCountDet) == 2) factor = 12;
3109 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3110 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3111 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3112 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3116 //Put default value negative
3117 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3118 fCurrentCoefE = 0.0;
3123 if(fDebugLevel > 1){
3125 if ( !fDebugStreamer ) {
3127 TDirectory *backup = gDirectory;
3128 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3129 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3132 Int_t detector = fCountDet;
3133 Int_t caligroup = idect;
3134 Short_t rowmin = fCalibraMode->GetRowMin(0);
3135 Short_t rowmax = fCalibraMode->GetRowMax(0);
3136 Short_t colmin = fCalibraMode->GetColMin(0);
3137 Short_t colmax = fCalibraMode->GetColMax(0);
3138 Float_t gf = fCurrentCoef[0];
3139 Float_t gfs = fCurrentCoef[1];
3140 Float_t gfE = fCurrentCoefE;
3142 (*fDebugStreamer) << "FillFillCH" <<
3143 "detector=" << detector <<
3144 "caligroup=" << caligroup <<
3145 "rowmin=" << rowmin <<
3146 "rowmax=" << rowmax <<
3147 "colmin=" << colmin <<
3148 "colmax=" << colmax <<
3156 for (Int_t k = 0; k < 2304; k++) {
3157 fCurrentCoefDetector[k] = 0.0;
3161 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3165 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
3167 // Calcul the coef from the database choosen
3168 CalculChargeCoefMean(kFALSE);
3170 //stack 2, not stack 2
3172 if(GetStack(fCountDet) == 2) factor = 12;
3175 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3176 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3177 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3178 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3182 //Put default value negative
3183 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3184 fCurrentCoefE = 0.0;
3193 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3194 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
3197 // For the case where there are not enough entries in the histograms
3198 // of the calibration group, the value present in the choosen database
3199 // will be put. A negativ sign enables to know that a fit was not possible.
3201 if (fDebugLevel == 1) {
3202 AliInfo("The element has not enough statistic to be fitted");
3204 else if (fNbDet > 0) {
3206 Int_t firstdetector = fCountDet;
3207 Int_t lastdetector = fCountDet+fNbDet;
3208 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3209 // loop over detectors
3210 for(Int_t det = firstdetector; det < lastdetector; det++){
3212 //Set the calibration object again
3216 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3218 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3219 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3220 ,(Int_t) GetStack(fCountDet)
3221 ,(Int_t) GetSector(fCountDet),1);
3222 // Reconstruct row min row max
3223 ReconstructFitRowMinRowMax(idect,1);
3225 // Calcul the coef from the database choosen for the detector
3226 CalculVdriftCoefMean();
3229 //stack 2, not stack 2
3231 if(GetStack(fCountDet) == 2) factor = 12;
3234 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3235 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3236 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3237 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3238 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3242 //Put default value negative
3243 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3244 fCurrentCoefE = 0.0;
3245 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3246 fCurrentCoefE2 = 0.0;
3252 if(fDebugLevel > 1){
3254 if ( !fDebugStreamer ) {
3256 TDirectory *backup = gDirectory;
3257 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3258 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3262 Int_t detector = fCountDet;
3263 Int_t caligroup = idect;
3264 Short_t rowmin = fCalibraMode->GetRowMin(1);
3265 Short_t rowmax = fCalibraMode->GetRowMax(1);
3266 Short_t colmin = fCalibraMode->GetColMin(1);
3267 Short_t colmax = fCalibraMode->GetColMax(1);
3268 Float_t vf = fCurrentCoef[0];
3269 Float_t vs = fCurrentCoef[1];
3270 Float_t vfE = fCurrentCoefE;
3271 Float_t t0f = fCurrentCoef2[0];
3272 Float_t t0s = fCurrentCoef2[1];
3273 Float_t t0E = fCurrentCoefE2;
3277 (* fDebugStreamer) << "FillFillPH"<<
3278 "detector="<<detector<<
3279 "nentries="<<nentries<<
3280 "caligroup="<<caligroup<<
3294 for (Int_t k = 0; k < 2304; k++) {
3295 fCurrentCoefDetector[k] = 0.0;
3296 fCurrentCoefDetector2[k] = 0.0;
3300 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3304 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3306 CalculVdriftCoefMean();
3309 //stack 2 and not stack 2
3311 if(GetStack(fCountDet) == 2) factor = 12;
3315 // Fill the fCurrentCoefDetector 2
3316 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3317 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3318 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3319 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3323 // Put the default value
3324 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3325 fCurrentCoefE = 0.0;
3326 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3327 fCurrentCoefE2 = 0.0;
3329 FillFillPH(idect,nentries);
3338 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3339 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3342 // For the case where there are not enough entries in the histograms
3343 // of the calibration group, the value present in the choosen database
3344 // will be put. A negativ sign enables to know that a fit was not possible.
3347 if (fDebugLevel == 1) {
3348 AliInfo("The element has not enough statistic to be fitted");
3350 else if (fNbDet > 0){
3352 Int_t firstdetector = fCountDet;
3353 Int_t lastdetector = fCountDet+fNbDet;
3354 // AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3356 // loop over detectors
3357 for(Int_t det = firstdetector; det < lastdetector; det++){
3359 //Set the calibration object again
3363 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3365 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3366 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3367 ,(Int_t) GetStack(fCountDet)
3368 ,(Int_t) GetSector(fCountDet),2);
3369 // Reconstruct row min row max
3370 ReconstructFitRowMinRowMax(idect,2);
3372 // Calcul the coef from the database choosen for the detector
3373 CalculPRFCoefMean();
3375 //stack 2, not stack 2
3377 if(GetStack(fCountDet) == 2) factor = 12;
3380 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3381 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3382 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3383 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3387 //Put default value negative
3388 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3389 fCurrentCoefE = 0.0;
3394 if(fDebugLevel > 1){
3396 if ( !fDebugStreamer ) {
3398 TDirectory *backup = gDirectory;
3399 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3400 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3403 Int_t detector = fCountDet;
3404 Int_t layer = GetLayer(fCountDet);
3405 Int_t caligroup = idect;
3406 Short_t rowmin = fCalibraMode->GetRowMin(2);
3407 Short_t rowmax = fCalibraMode->GetRowMax(2);
3408 Short_t colmin = fCalibraMode->GetColMin(2);
3409 Short_t colmax = fCalibraMode->GetColMax(2);
3410 Float_t widf = fCurrentCoef[0];
3411 Float_t wids = fCurrentCoef[1];
3412 Float_t widfE = fCurrentCoefE;
3414 (* fDebugStreamer) << "FillFillPRF"<<
3415 "detector="<<detector<<
3417 "caligroup="<<caligroup<<
3428 for (Int_t k = 0; k < 2304; k++) {
3429 fCurrentCoefDetector[k] = 0.0;
3433 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3437 // AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3439 CalculPRFCoefMean();
3441 // stack 2 and not stack 2
3443 if(GetStack(fCountDet) == 2) factor = 12;
3447 // Fill the fCurrentCoefDetector
3448 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3449 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3450 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3454 // Put the default value
3455 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3456 fCurrentCoefE = 0.0;
3464 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3465 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3468 // For the case where there are not enough entries in the histograms
3469 // of the calibration group, the value present in the choosen database
3470 // will be put. A negativ sign enables to know that a fit was not possible.
3473 // Calcul the coef from the database choosen
3474 CalculVdriftLorentzCoef();
3477 if(GetStack(fCountDet) == 2) factor = 1728;
3481 // Fill the fCurrentCoefDetector
3482 for (Int_t k = 0; k < factor; k++) {
3483 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3484 // should be negative
3485 fCurrentCoefDetector2[k] = fCurrentCoef2[1]+100.0;
3489 //Put default opposite sign only for vdrift
3490 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3491 fCurrentCoefE = 0.0;
3492 fCurrentCoef2[0] = fCurrentCoef2[1]+100.0;
3493 fCurrentCoefE2 = 0.0;
3495 FillFillLinearFitter();
3500 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3501 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3504 // Fill the coefficients found with the fits or other
3505 // methods from the Fit functions
3508 if (fDebugLevel != 1) {
3510 Int_t firstdetector = fCountDet;
3511 Int_t lastdetector = fCountDet+fNbDet;
3512 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3513 // loop over detectors
3514 for(Int_t det = firstdetector; det < lastdetector; det++){
3516 //Set the calibration object again
3520 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3522 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3523 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3524 ,(Int_t) GetStack(fCountDet)
3525 ,(Int_t) GetSector(fCountDet),0);
3526 // Reconstruct row min row max
3527 ReconstructFitRowMinRowMax(idect,0);
3529 // Calcul the coef from the database choosen for the detector
3530 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3531 else CalculChargeCoefMean(kTRUE);
3533 //stack 2, not stack 2
3535 if(GetStack(fCountDet) == 2) factor = 12;
3538 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3539 Double_t coeftoput = 1.0;
3540 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3541 else coeftoput = fCurrentCoef[0];
3542 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3543 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3544 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3551 if(fDebugLevel > 1){
3553 if ( !fDebugStreamer ) {
3555 TDirectory *backup = gDirectory;
3556 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3557 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3560 Int_t detector = fCountDet;
3561 Int_t caligroup = idect;
3562 Short_t rowmin = fCalibraMode->GetRowMin(0);
3563 Short_t rowmax = fCalibraMode->GetRowMax(0);
3564 Short_t colmin = fCalibraMode->GetColMin(0);
3565 Short_t colmax = fCalibraMode->GetColMax(0);
3566 Float_t gf = fCurrentCoef[0];
3567 Float_t gfs = fCurrentCoef[1];
3568 Float_t gfE = fCurrentCoefE;
3570 (*fDebugStreamer) << "FillFillCH" <<
3571 "detector=" << detector <<
3572 "caligroup=" << caligroup <<
3573 "rowmin=" << rowmin <<
3574 "rowmax=" << rowmax <<
3575 "colmin=" << colmin <<
3576 "colmax=" << colmax <<
3584 for (Int_t k = 0; k < 2304; k++) {
3585 fCurrentCoefDetector[k] = 0.0;
3589 //printf("Check the count now: fCountDet %d\n",fCountDet);
3594 if(GetStack(fCountDet) == 2) factor = 12;
3597 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3598 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3599 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3610 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3611 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3614 // Fill the coefficients found with the fits or other
3615 // methods from the Fit functions
3618 if (fDebugLevel != 1) {
3621 Int_t firstdetector = fCountDet;
3622 Int_t lastdetector = fCountDet+fNbDet;
3623 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3625 // loop over detectors
3626 for(Int_t det = firstdetector; det < lastdetector; det++){
3628 //Set the calibration object again
3632 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3634 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3635 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3636 ,(Int_t) GetStack(fCountDet)
3637 ,(Int_t) GetSector(fCountDet),1);
3638 // Reconstruct row min row max
3639 ReconstructFitRowMinRowMax(idect,1);
3641 // Calcul the coef from the database choosen for the detector
3642 CalculVdriftCoefMean();
3645 //stack 2, not stack 2
3647 if(GetStack(fCountDet) == 2) factor = 12;
3650 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3651 Double_t coeftoput = 1.5;
3652 Double_t coeftoput2 = 0.0;
3654 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3655 else coeftoput = fCurrentCoef[0];
3657 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3658 else coeftoput2 = fCurrentCoef2[0];
3660 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3661 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3662 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3663 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3671 if(fDebugLevel > 1){
3673 if ( !fDebugStreamer ) {
3675 TDirectory *backup = gDirectory;
3676 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3677 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3681 Int_t detector = fCountDet;
3682 Int_t caligroup = idect;
3683 Short_t rowmin = fCalibraMode->GetRowMin(1);
3684 Short_t rowmax = fCalibraMode->GetRowMax(1);
3685 Short_t colmin = fCalibraMode->GetColMin(1);
3686 Short_t colmax = fCalibraMode->GetColMax(1);
3687 Float_t vf = fCurrentCoef[0];
3688 Float_t vs = fCurrentCoef[1];
3689 Float_t vfE = fCurrentCoefE;
3690 Float_t t0f = fCurrentCoef2[0];
3691 Float_t t0s = fCurrentCoef2[1];
3692 Float_t t0E = fCurrentCoefE2;
3696 (* fDebugStreamer) << "FillFillPH"<<
3697 "detector="<<detector<<
3698 "nentries="<<nentries<<
3699 "caligroup="<<caligroup<<
3713 for (Int_t k = 0; k < 2304; k++) {
3714 fCurrentCoefDetector[k] = 0.0;
3715 fCurrentCoefDetector2[k] = 0.0;
3719 //printf("Check the count now: fCountDet %d\n",fCountDet);
3724 if(GetStack(fCountDet) == 2) factor = 12;
3727 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3728 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3729 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3730 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3734 FillFillPH(idect,nentries);
3739 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3740 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3743 // Fill the coefficients found with the fits or other
3744 // methods from the Fit functions
3747 if (fDebugLevel != 1) {
3750 Int_t firstdetector = fCountDet;
3751 Int_t lastdetector = fCountDet+fNbDet;
3752 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3754 // loop over detectors
3755 for(Int_t det = firstdetector; det < lastdetector; det++){
3757 //Set the calibration object again
3761 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3763 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3764 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3765 ,(Int_t) GetStack(fCountDet)
3766 ,(Int_t) GetSector(fCountDet),2);
3767 // Reconstruct row min row max
3768 ReconstructFitRowMinRowMax(idect,2);
3770 // Calcul the coef from the database choosen for the detector
3771 CalculPRFCoefMean();
3773 //stack 2, not stack 2
3775 if(GetStack(fCountDet) == 2) factor = 12;
3778 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3779 Double_t coeftoput = 1.0;
3780 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3781 else coeftoput = fCurrentCoef[0];
3782 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3783 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3784 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3791 if(fDebugLevel > 1){
3793 if ( !fDebugStreamer ) {
3795 TDirectory *backup = gDirectory;
3796 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3797 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3800 Int_t detector = fCountDet;
3801 Int_t layer = GetLayer(fCountDet);
3802 Int_t caligroup = idect;
3803 Short_t rowmin = fCalibraMode->GetRowMin(2);
3804 Short_t rowmax = fCalibraMode->GetRowMax(2);
3805 Short_t colmin = fCalibraMode->GetColMin(2);
3806 Short_t colmax = fCalibraMode->GetColMax(2);
3807 Float_t widf = fCurrentCoef[0];
3808 Float_t wids = fCurrentCoef[1];
3809 Float_t widfE = fCurrentCoefE;
3811 (* fDebugStreamer) << "FillFillPRF"<<
3812 "detector="<<detector<<
3814 "caligroup="<<caligroup<<
3825 for (Int_t k = 0; k < 2304; k++) {
3826 fCurrentCoefDetector[k] = 0.0;
3830 //printf("Check the count now: fCountDet %d\n",fCountDet);
3835 if(GetStack(fCountDet) == 2) factor = 12;
3838 // Pointer to the branch
3839 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3840 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3841 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3851 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3852 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3855 // Fill the coefficients found with the fits or other
3856 // methods from the Fit functions
3860 if(GetStack(fCountDet) == 2) factor = 1728;
3863 // Pointer to the branch
3864 for (Int_t k = 0; k < factor; k++) {
3865 fCurrentCoefDetector[k] = fCurrentCoef[0];
3866 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3869 FillFillLinearFitter();
3874 //________________________________________________________________________________
3875 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3878 // DebugStream and fVectorFit
3881 // End of one detector
3882 if ((idect == (fCount-1))) {
3885 for (Int_t k = 0; k < 2304; k++) {
3886 fCurrentCoefDetector[k] = 0.0;
3890 if(fDebugLevel > 1){
3892 if ( !fDebugStreamer ) {
3894 TDirectory *backup = gDirectory;
3895 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3896 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3899 Int_t detector = fCountDet;
3900 Int_t caligroup = idect;
3901 Short_t rowmin = fCalibraMode->GetRowMin(0);
3902 Short_t rowmax = fCalibraMode->GetRowMax(0);
3903 Short_t colmin = fCalibraMode->GetColMin(0);
3904 Short_t colmax = fCalibraMode->GetColMax(0);
3905 Float_t gf = fCurrentCoef[0];
3906 Float_t gfs = fCurrentCoef[1];
3907 Float_t gfE = fCurrentCoefE;
3909 (*fDebugStreamer) << "FillFillCH" <<
3910 "detector=" << detector <<
3911 "caligroup=" << caligroup <<
3912 "rowmin=" << rowmin <<
3913 "rowmax=" << rowmax <<
3914 "colmin=" << colmin <<
3915 "colmax=" << colmax <<
3923 //________________________________________________________________________________
3924 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3927 // DebugStream and fVectorFit and fVectorFit2
3930 // End of one detector
3931 if ((idect == (fCount-1))) {
3935 for (Int_t k = 0; k < 2304; k++) {
3936 fCurrentCoefDetector[k] = 0.0;
3937 fCurrentCoefDetector2[k] = 0.0;
3941 if(fDebugLevel > 1){
3943 if ( !fDebugStreamer ) {
3945 TDirectory *backup = gDirectory;
3946 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3947 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3951 Int_t detector = fCountDet;
3952 Int_t caligroup = idect;
3953 Short_t rowmin = fCalibraMode->GetRowMin(1);
3954 Short_t rowmax = fCalibraMode->GetRowMax(1);
3955 Short_t colmin = fCalibraMode->GetColMin(1);
3956 Short_t colmax = fCalibraMode->GetColMax(1);
3957 Float_t vf = fCurrentCoef[0];
3958 Float_t vs = fCurrentCoef[1];
3959 Float_t vfE = fCurrentCoefE;
3960 Float_t t0f = fCurrentCoef2[0];
3961 Float_t t0s = fCurrentCoef2[1];
3962 Float_t t0E = fCurrentCoefE2;
3966 (* fDebugStreamer) << "FillFillPH"<<
3967 "detector="<<detector<<
3968 "nentries="<<nentries<<
3969 "caligroup="<<caligroup<<
3984 //________________________________________________________________________________
3985 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3988 // DebugStream and fVectorFit
3991 // End of one detector
3992 if ((idect == (fCount-1))) {
3995 for (Int_t k = 0; k < 2304; k++) {
3996 fCurrentCoefDetector[k] = 0.0;
4001 if(fDebugLevel > 1){
4003 if ( !fDebugStreamer ) {
4005 TDirectory *backup = gDirectory;
4006 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4007 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4010 Int_t detector = fCountDet;
4011 Int_t layer = GetLayer(fCountDet);
4012 Int_t caligroup = idect;
4013 Short_t rowmin = fCalibraMode->GetRowMin(2);
4014 Short_t rowmax = fCalibraMode->GetRowMax(2);
4015 Short_t colmin = fCalibraMode->GetColMin(2);
4016 Short_t colmax = fCalibraMode->GetColMax(2);
4017 Float_t widf = fCurrentCoef[0];
4018 Float_t wids = fCurrentCoef[1];
4019 Float_t widfE = fCurrentCoefE;
4021 (* fDebugStreamer) << "FillFillPRF"<<
4022 "detector="<<detector<<
4024 "caligroup="<<caligroup<<
4036 //________________________________________________________________________________
4037 void AliTRDCalibraFit::FillFillLinearFitter()
4040 // DebugStream and fVectorFit
4043 // End of one detector
4049 for (Int_t k = 0; k < 2304; k++) {
4050 fCurrentCoefDetector[k] = 0.0;
4051 fCurrentCoefDetector2[k] = 0.0;
4055 if(fDebugLevel > 1){
4057 if ( !fDebugStreamer ) {
4059 TDirectory *backup = gDirectory;
4060 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
4061 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4064 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4065 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4066 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4067 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4068 Float_t tiltangle = padplane->GetTiltingAngle();
4069 Int_t detector = fCountDet;
4070 Int_t stack = GetStack(fCountDet);
4071 Int_t layer = GetLayer(fCountDet);
4072 Float_t vf = fCurrentCoef[0];
4073 Float_t vs = fCurrentCoef[1];
4074 Float_t vfE = fCurrentCoefE;
4075 Float_t lorentzangler = fCurrentCoef2[0];
4076 Float_t elorentzangler = fCurrentCoefE2;
4077 Float_t lorentzangles = fCurrentCoef2[1];
4079 (* fDebugStreamer) << "FillFillLinearFitter"<<
4080 "detector="<<detector<<
4085 "tiltangle="<<tiltangle<<
4089 "lorentzangler="<<lorentzangler<<
4090 "Elorentzangler="<<elorentzangler<<
4091 "lorentzangles="<<lorentzangles<<
4097 //____________Calcul Coef Mean_________________________________________________
4099 //_____________________________________________________________________________
4100 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
4103 // For the detector Dect calcul the mean time 0
4104 // for the calibration group idect from the choosen database
4107 fCurrentCoef2[1] = 0.0;
4108 if(fDebugLevel != 1){
4109 if(((fCalibraMode->GetNz(1) > 0) ||
4110 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4112 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4113 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4114 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4118 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4124 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4128 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
4129 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
4130 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4133 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
4141 //_____________________________________________________________________________
4142 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
4145 // For the detector Dect calcul the mean gain factor
4146 // for the calibration group idect from the choosen database
4149 fCurrentCoef[1] = 0.0;
4150 if(fDebugLevel != 1){
4151 if (((fCalibraMode->GetNz(0) > 0) ||
4152 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
4153 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4154 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4155 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4156 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4159 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4163 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4164 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
4169 //_____________________________________________________________________________
4170 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
4173 // For the detector Dect calcul the mean sigma of pad response
4174 // function for the calibration group idect from the choosen database
4177 fCurrentCoef[1] = 0.0;
4178 if(fDebugLevel != 1){
4179 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4180 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4181 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
4184 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
4188 //_____________________________________________________________________________
4189 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
4192 // For the detector dect calcul the mean drift velocity for the
4193 // calibration group idect from the choosen database
4196 fCurrentCoef[1] = 0.0;
4197 if(fDebugLevel != 1){
4198 if (((fCalibraMode->GetNz(1) > 0) ||
4199 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4201 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4202 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4203 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4207 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4212 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4217 //_____________________________________________________________________________
4218 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4221 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4224 fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet);
4225 fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet);
4229 //_____________________________________________________________________________
4230 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4233 // Default width of the PRF if there is no database as reference
4238 //case 0: return 0.515;
4239 //case 1: return 0.502;
4240 //case 2: return 0.491;
4241 //case 3: return 0.481;
4242 //case 4: return 0.471;
4243 //case 5: return 0.463;
4244 //default: return 0.0;
4247 case 0: return 0.538429;
4248 case 1: return 0.524302;
4249 case 2: return 0.511591;
4250 case 3: return 0.500140;
4251 case 4: return 0.489821;
4252 case 5: return 0.480524;
4253 default: return 0.0;
4256 //________________________________________________________________________________
4257 void AliTRDCalibraFit::SetCalROC(Int_t i)
4260 // Set the calib object for fCountDet
4263 Float_t value = 0.0;
4265 //Get the CalDet object
4267 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4269 AliInfo("Could not get calibDB");
4276 fCalROC->~AliTRDCalROC();
4277 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4278 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4282 fCalROC->~AliTRDCalROC();
4283 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4284 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4286 fCalROC2->~AliTRDCalROC();
4287 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4288 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4292 fCalROC->~AliTRDCalROC();
4293 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4294 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4303 if(fCalROC) delete fCalROC;
4304 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4305 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4306 fCalROC->SetValue(k,1.0);
4310 if(fCalROC) delete fCalROC;
4311 if(fCalROC2) delete fCalROC2;
4312 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4313 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4314 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4315 fCalROC->SetValue(k,1.0);
4316 fCalROC2->SetValue(k,0.0);
4320 if(fCalROC) delete fCalROC;
4321 value = GetPRFDefault(GetLayer(fCountDet));
4322 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4323 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4324 fCalROC->SetValue(k,value);
4332 //____________Fit Methods______________________________________________________
4334 //_____________________________________________________________________________
4335 void AliTRDCalibraFit::FitPente(TH1* projPH)
4338 // Slope methode for the drift velocity
4342 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4349 fCurrentCoefE = 0.0;
4350 fCurrentCoefE2 = 0.0;
4351 fCurrentCoef[0] = 0.0;
4352 fCurrentCoef2[0] = 0.0;
4353 TLine *line = new TLine();
4356 TAxis *xpph = projPH->GetXaxis();
4357 Int_t nbins = xpph->GetNbins();
4358 Double_t lowedge = xpph->GetBinLowEdge(1);
4359 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4360 Double_t widbins = (upedge - lowedge) / nbins;
4361 Double_t limit = upedge + 0.5 * widbins;
4364 // Beginning of the signal
4365 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4366 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4367 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4369 binmax = (Int_t) pentea->GetMaximumBin();
4372 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4374 if (binmax >= nbins) {
4377 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4379 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4380 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4381 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4382 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4383 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4384 if (TMath::Abs(l3P2am) > 0.00000001) {
4385 fPhd[0] = -(l3P1am / (2 * l3P2am));
4388 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4389 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4392 // Amplification region
4395 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4396 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4403 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4405 if (binmax >= nbins) {
4408 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4410 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4411 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4412 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4413 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4414 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4415 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4416 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4418 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4419 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4422 fCurrentCoefE2 = fCurrentCoefE;
4425 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4426 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4427 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4430 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4433 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4435 if (binmin >= nbins) {
4438 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4443 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4444 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4445 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4446 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4447 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4448 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4449 if (TMath::Abs(l3P2dr) > 0.00000001) {
4450 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4452 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4453 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4455 Float_t fPhdt0 = 0.0;
4456 Float_t t0Shift = 0.0;
4459 t0Shift = fT0Shift1;
4463 t0Shift = fT0Shift0;
4466 if ((fPhd[2] > fPhd[0]) &&
4467 (fPhd[2] > fPhd[1]) &&
4468 (fPhd[1] > fPhd[0]) &&
4470 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4471 fNumberFitSuccess++;
4473 if (fPhdt0 >= 0.0) {
4474 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4475 if (fCurrentCoef2[0] < -3.0) {
4476 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4480 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4485 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4486 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4489 if (fDebugLevel == 1) {
4490 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4493 line->SetLineColor(2);
4494 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4495 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4496 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4497 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4498 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4499 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4500 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4501 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4504 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4513 //_____________________________________________________________________________
4514 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4517 // Slope methode but with polynomes de Lagrange
4521 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4524 //Double_t *x = new Double_t[5];
4525 //Double_t *y = new Double_t[5];
4542 fCurrentCoefE = 0.0;
4543 fCurrentCoefE2 = 1.0;
4544 fCurrentCoef[0] = 0.0;
4545 fCurrentCoef2[0] = 0.0;
4546 TLine *line = new TLine();
4547 TF1 * polynome = 0x0;
4548 TF1 * polynomea = 0x0;
4549 TF1 * polynomeb = 0x0;
4557 TAxis *xpph = projPH->GetXaxis();
4558 Int_t nbins = xpph->GetNbins();
4559 Double_t lowedge = xpph->GetBinLowEdge(1);
4560 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4561 Double_t widbins = (upedge - lowedge) / nbins;
4562 Double_t limit = upedge + 0.5 * widbins;
4567 // Beginning of the signal
4568 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4569 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4570 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4573 binmax = (Int_t) pentea->GetMaximumBin();
4575 Double_t minnn = 0.0;
4576 Double_t maxxx = 0.0;
4578 Int_t kase = nbins-binmax;
4586 minnn = pentea->GetBinCenter(binmax-2);
4587 maxxx = pentea->GetBinCenter(binmax);
4588 x[0] = pentea->GetBinCenter(binmax-2);
4589 x[1] = pentea->GetBinCenter(binmax-1);
4590 x[2] = pentea->GetBinCenter(binmax);
4591 y[0] = pentea->GetBinContent(binmax-2);
4592 y[1] = pentea->GetBinContent(binmax-1);
4593 y[2] = pentea->GetBinContent(binmax);
4594 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4595 AliInfo("At the limit for beginning!");
4598 minnn = pentea->GetBinCenter(binmax-2);
4599 maxxx = pentea->GetBinCenter(binmax+1);
4600 x[0] = pentea->GetBinCenter(binmax-2);
4601 x[1] = pentea->GetBinCenter(binmax-1);
4602 x[2] = pentea->GetBinCenter(binmax);
4603 x[3] = pentea->GetBinCenter(binmax+1);
4604 y[0] = pentea->GetBinContent(binmax-2);
4605 y[1] = pentea->GetBinContent(binmax-1);
4606 y[2] = pentea->GetBinContent(binmax);
4607 y[3] = pentea->GetBinContent(binmax+1);
4608 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4616 minnn = pentea->GetBinCenter(binmax);
4617 maxxx = pentea->GetBinCenter(binmax+2);
4618 x[0] = pentea->GetBinCenter(binmax);
4619 x[1] = pentea->GetBinCenter(binmax+1);
4620 x[2] = pentea->GetBinCenter(binmax+2);
4621 y[0] = pentea->GetBinContent(binmax);
4622 y[1] = pentea->GetBinContent(binmax+1);
4623 y[2] = pentea->GetBinContent(binmax+2);
4624 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4627 minnn = pentea->GetBinCenter(binmax-1);
4628 maxxx = pentea->GetBinCenter(binmax+2);
4629 x[0] = pentea->GetBinCenter(binmax-1);
4630 x[1] = pentea->GetBinCenter(binmax);
4631 x[2] = pentea->GetBinCenter(binmax+1);
4632 x[3] = pentea->GetBinCenter(binmax+2);
4633 y[0] = pentea->GetBinContent(binmax-1);
4634 y[1] = pentea->GetBinContent(binmax);
4635 y[2] = pentea->GetBinContent(binmax+1);
4636 y[3] = pentea->GetBinContent(binmax+2);
4637 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4640 minnn = pentea->GetBinCenter(binmax-2);
4641 maxxx = pentea->GetBinCenter(binmax+2);
4642 x[0] = pentea->GetBinCenter(binmax-2);
4643 x[1] = pentea->GetBinCenter(binmax-1);
4644 x[2] = pentea->GetBinCenter(binmax);
4645 x[3] = pentea->GetBinCenter(binmax+1);
4646 x[4] = pentea->GetBinCenter(binmax+2);
4647 y[0] = pentea->GetBinContent(binmax-2);
4648 y[1] = pentea->GetBinContent(binmax-1);
4649 y[2] = pentea->GetBinContent(binmax);
4650 y[3] = pentea->GetBinContent(binmax+1);
4651 y[4] = pentea->GetBinContent(binmax+2);
4652 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4660 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4661 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4663 Double_t step = (maxxx-minnn)/10000;
4665 Double_t maxvalue = 0.0;
4666 Double_t placemaximum = minnn;
4667 for(Int_t o = 0; o < 10000; o++){
4668 if(o == 0) maxvalue = polynomeb->Eval(l);
4669 if(maxvalue < (polynomeb->Eval(l))){
4670 maxvalue = polynomeb->Eval(l);
4675 fPhd[0] = placemaximum;
4678 // Amplification region
4681 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4682 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4688 Double_t minn = 0.0;
4689 Double_t maxx = 0.0;
4701 Int_t kase1 = nbins - binmax;
4703 //Determination of minn and maxx
4704 //case binmax = nbins
4709 minn = projPH->GetBinCenter(binmax-2);
4710 maxx = projPH->GetBinCenter(binmax);
4711 x[0] = projPH->GetBinCenter(binmax-2);
4712 x[1] = projPH->GetBinCenter(binmax-1);
4713 x[2] = projPH->GetBinCenter(binmax);
4714 y[0] = projPH->GetBinContent(binmax-2);
4715 y[1] = projPH->GetBinContent(binmax-1);
4716 y[2] = projPH->GetBinContent(binmax);
4717 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4718 //AliInfo("At the limit for the drift!");
4721 minn = projPH->GetBinCenter(binmax-2);
4722 maxx = projPH->GetBinCenter(binmax+1);
4723 x[0] = projPH->GetBinCenter(binmax-2);
4724 x[1] = projPH->GetBinCenter(binmax-1);
4725 x[2] = projPH->GetBinCenter(binmax);
4726 x[3] = projPH->GetBinCenter(binmax+1);
4727 y[0] = projPH->GetBinContent(binmax-2);
4728 y[1] = projPH->GetBinContent(binmax-1);
4729 y[2] = projPH->GetBinContent(binmax);
4730 y[3] = projPH->GetBinContent(binmax+1);
4731 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4740 minn = projPH->GetBinCenter(binmax);
4741 maxx = projPH->GetBinCenter(binmax+2);
4742 x[0] = projPH->GetBinCenter(binmax);
4743 x[1] = projPH->GetBinCenter(binmax+1);
4744 x[2] = projPH->GetBinCenter(binmax+2);
4745 y[0] = projPH->GetBinContent(binmax);
4746 y[1] = projPH->GetBinContent(binmax+1);
4747 y[2] = projPH->GetBinContent(binmax+2);
4748 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4751 minn = projPH->GetBinCenter(binmax-1);
4752 maxx = projPH->GetBinCenter(binmax+2);
4753 x[0] = projPH->GetBinCenter(binmax-1);
4754 x[1] = projPH->GetBinCenter(binmax);
4755 x[2] = projPH->GetBinCenter(binmax+1);
4756 x[3] = projPH->GetBinCenter(binmax+2);
4757 y[0] = projPH->GetBinContent(binmax-1);
4758 y[1] = projPH->GetBinContent(binmax);
4759 y[2] = projPH->GetBinContent(binmax+1);
4760 y[3] = projPH->GetBinContent(binmax+2);
4761 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4764 minn = projPH->GetBinCenter(binmax-2);
4765 maxx = projPH->GetBinCenter(binmax+2);
4766 x[0] = projPH->GetBinCenter(binmax-2);
4767 x[1] = projPH->GetBinCenter(binmax-1);
4768 x[2] = projPH->GetBinCenter(binmax);
4769 x[3] = projPH->GetBinCenter(binmax+1);
4770 x[4] = projPH->GetBinCenter(binmax+2);
4771 y[0] = projPH->GetBinContent(binmax-2);
4772 y[1] = projPH->GetBinContent(binmax-1);
4773 y[2] = projPH->GetBinContent(binmax);
4774 y[3] = projPH->GetBinContent(binmax+1);
4775 y[4] = projPH->GetBinContent(binmax+2);
4776 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4783 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4784 polynomea->SetParameters(c0,c1,c2,c3,c4);
4786 Double_t step = (maxx-minn)/1000;
4788 Double_t maxvalue = 0.0;
4789 Double_t placemaximum = minn;
4790 for(Int_t o = 0; o < 1000; o++){
4791 if(o == 0) maxvalue = polynomea->Eval(l);
4792 if(maxvalue < (polynomea->Eval(l))){
4793 maxvalue = polynomea->Eval(l);
4798 fPhd[1] = placemaximum;
4802 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4803 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4804 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4807 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4813 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4817 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4818 AliInfo("Too many fluctuations at the end!");
4821 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4822 AliInfo("Too many fluctuations at the end!");
4825 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4826 AliInfo("No entries for the next bin!");
4827 pente->SetBinContent(binmin,0);
4828 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4844 Bool_t case1 = kFALSE;
4845 Bool_t case2 = kFALSE;
4846 Bool_t case4 = kFALSE;
4848 //Determination of min and max
4849 //case binmin <= nbins-3
4851 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4852 min = pente->GetBinCenter(binmin-2);
4853 max = pente->GetBinCenter(binmin+2);
4854 x[0] = pente->GetBinCenter(binmin-2);
4855 x[1] = pente->GetBinCenter(binmin-1);
4856 x[2] = pente->GetBinCenter(binmin);
4857 x[3] = pente->GetBinCenter(binmin+1);
4858 x[4] = pente->GetBinCenter(binmin+2);
4859 y[0] = pente->GetBinContent(binmin-2);
4860 y[1] = pente->GetBinContent(binmin-1);
4861 y[2] = pente->GetBinContent(binmin);
4862 y[3] = pente->GetBinContent(binmin+1);
4863 y[4] = pente->GetBinContent(binmin+2);
4864 //Calcul the polynome de Lagrange
4865 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4867 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4868 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4869 //AliInfo("polynome 4 false 1");
4872 if(((binmin+3) <= (nbins-1)) &&
4873 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4874 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4875 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4876 AliInfo("polynome 4 false 2");
4880 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4881 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4882 //AliInfo("polynome 4 case 1");
4885 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4886 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4887 //AliInfo("polynome 4 case 4");
4892 //case binmin = nbins-2
4894 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4896 min = pente->GetBinCenter(binmin-2);
4897 max = pente->GetBinCenter(binmin+1);
4898 x[0] = pente->GetBinCenter(binmin-2);
4899 x[1] = pente->GetBinCenter(binmin-1);
4900 x[2] = pente->GetBinCenter(binmin);
4901 x[3] = pente->GetBinCenter(binmin+1);
4902 y[0] = pente->GetBinContent(binmin-2);
4903 y[1] = pente->GetBinContent(binmin-1);
4904 y[2] = pente->GetBinContent(binmin);
4905 y[3] = pente->GetBinContent(binmin+1);
4906 //Calcul the polynome de Lagrange
4907 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4908 //richtung +: nothing
4910 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4911 //AliInfo("polynome 3- case 2");
4916 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4918 min = pente->GetBinCenter(binmin-1);
4919 max = pente->GetBinCenter(binmin+2);
4920 x[0] = pente->GetBinCenter(binmin-1);
4921 x[1] = pente->GetBinCenter(binmin);
4922 x[2] = pente->GetBinCenter(binmin+1);
4923 x[3] = pente->GetBinCenter(binmin+2);
4924 y[0] = pente->GetBinContent(binmin-1);
4925 y[1] = pente->GetBinContent(binmin);
4926 y[2] = pente->GetBinContent(binmin+1);
4927 y[3] = pente->GetBinContent(binmin+2);
4928 //Calcul the polynome de Lagrange
4929 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4931 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4932 //AliInfo("polynome 3+ case 2");
4937 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4938 min = pente->GetBinCenter(binmin);
4939 max = pente->GetBinCenter(binmin+2);
4940 x[0] = pente->GetBinCenter(binmin);
4941 x[1] = pente->GetBinCenter(binmin+1);
4942 x[2] = pente->GetBinCenter(binmin+2);
4943 y[0] = pente->GetBinContent(binmin);
4944 y[1] = pente->GetBinContent(binmin+1);
4945 y[2] = pente->GetBinContent(binmin+2);
4946 //Calcul the polynome de Lagrange
4947 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4949 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4950 //AliInfo("polynome 2+ false");
4955 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4957 min = pente->GetBinCenter(binmin-1);
4958 max = pente->GetBinCenter(binmin+1);
4959 x[0] = pente->GetBinCenter(binmin-1);
4960 x[1] = pente->GetBinCenter(binmin);
4961 x[2] = pente->GetBinCenter(binmin+1);
4962 y[0] = pente->GetBinContent(binmin-1);
4963 y[1] = pente->GetBinContent(binmin);
4964 y[2] = pente->GetBinContent(binmin+1);
4965 //Calcul the polynome de Lagrange
4966 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4967 //richtung +: nothing
4968 //richtung -: nothing
4970 //case binmin = nbins-1
4972 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4973 min = pente->GetBinCenter(binmin-2);
4974 max = pente->GetBinCenter(binmin);
4975 x[0] = pente->GetBinCenter(binmin-2);
4976 x[1] = pente->GetBinCenter(binmin-1);
4977 x[2] = pente->GetBinCenter(binmin);
4978 y[0] = pente->GetBinContent(binmin-2);
4979 y[1] = pente->GetBinContent(binmin-1);
4980 y[2] = pente->GetBinContent(binmin);
4981 //Calcul the polynome de Lagrange
4982 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4983 //AliInfo("At the limit for the drift!");
4984 //fluctuation too big!
4985 //richtung +: nothing
4987 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4988 //AliInfo("polynome 2- false ");
4992 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4994 AliInfo("At the limit for the drift and not usable!");
4998 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5000 AliInfo("For the drift...problem!");
5002 //pass but should not happen
5003 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
5005 AliInfo("For the drift...problem!");
5009 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5010 polynome->SetParameters(c0,c1,c2,c3,c4);
5011 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5012 Double_t step = (max-min)/1000;
5014 Double_t minvalue = 0.0;
5015 Double_t placeminimum = min;
5016 for(Int_t o = 0; o < 1000; o++){
5017 if(o == 0) minvalue = polynome->Eval(l);
5018 if(minvalue > (polynome->Eval(l))){
5019 minvalue = polynome->Eval(l);
5024 fPhd[2] = placeminimum;
5026 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
5027 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
5028 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
5030 Float_t fPhdt0 = 0.0;
5031 Float_t t0Shift = 0.0;
5034 t0Shift = fT0Shift1;
5038 t0Shift = fT0Shift0;
5041 if ((fPhd[2] > fPhd[0]) &&
5042 (fPhd[2] > fPhd[1]) &&
5043 (fPhd[1] > fPhd[0]) &&
5045 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5046 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5047 else fNumberFitSuccess++;
5048 if (fPhdt0 >= 0.0) {
5049 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5050 //printf("Value of timeoffset %f\n",fCurrentCoef2[0]);
5051 if (fCurrentCoef2[0] < -3.0) {
5052 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5056 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5060 ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
5061 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5063 if((fPhd[1] > fPhd[0]) &&
5065 if (fPhdt0 >= 0.0) {
5066 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5067 if (fCurrentCoef2[0] < -3.0) {
5068 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5072 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5076 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5077 //printf("Fit failed!\n");
5081 if (fDebugLevel == 1) {
5082 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5085 line->SetLineColor(2);
5086 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5087 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5088 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5089 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5090 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5091 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5092 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
5093 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5096 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5103 if(polynome) delete polynome;
5104 if(polynomea) delete polynomea;
5105 if(polynomeb) delete polynomeb;
5106 //if(x) delete [] x;
5107 //if(y) delete [] y;
5108 if(line) delete line;
5113 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5114 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5115 //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]);
5116 projPH->SetDirectory(0);
5120 //_____________________________________________________________________________
5121 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5124 // Fit methode for the drift velocity
5128 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5131 TAxis *xpph = projPH->GetXaxis();
5132 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5134 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5135 fPH->SetParameter(0,0.469); // Scaling
5136 fPH->SetParameter(1,0.18); // Start
5137 fPH->SetParameter(2,0.0857325); // AR
5138 fPH->SetParameter(3,1.89); // DR
5139 fPH->SetParameter(4,0.08); // QA/QD
5140 fPH->SetParameter(5,0.0); // Baseline
5142 TLine *line = new TLine();
5144 fCurrentCoef[0] = 0.0;
5145 fCurrentCoef2[0] = 0.0;
5146 fCurrentCoefE = 0.0;
5147 fCurrentCoefE2 = 0.0;
5149 if (idect%fFitPHPeriode == 0) {
5151 AliInfo(Form("The detector %d will be fitted",idect));
5152 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5153 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5154 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5155 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5156 fPH->SetParameter(4,0.225); // QA/QD
5157 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5159 if (fDebugLevel != 1) {
5160 projPH->Fit(fPH,"0M","",0.0,upedge);
5163 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5165 projPH->Fit(fPH,"M+","",0.0,upedge);
5167 line->SetLineColor(4);
5168 line->DrawLine(fPH->GetParameter(1)
5170 ,fPH->GetParameter(1)
5171 ,projPH->GetMaximum());
5172 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5174 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5175 ,projPH->GetMaximum());
5176 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5178 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5179 ,projPH->GetMaximum());
5182 if (fPH->GetParameter(3) != 0) {
5183 fNumberFitSuccess++;
5184 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
5185 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
5186 fCurrentCoef2[0] = fPH->GetParameter(1);
5187 fCurrentCoefE2 = fPH->GetParError(1);
5190 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5191 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5197 // Put the default value
5198 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5199 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5202 if (fDebugLevel != 1) {
5207 //_____________________________________________________________________________
5208 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5211 // Fit methode for the sigma of the pad response function
5216 fCurrentCoef[0] = 0.0;
5217 fCurrentCoefE = 0.0;
5219 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5221 if(TMath::Abs(ret+4) <= 0.000000001){
5222 fCurrentCoef[0] = -fCurrentCoef[1];
5226 fNumberFitSuccess++;
5227 fCurrentCoef[0] = param[2];
5228 fCurrentCoefE = ret;
5232 //_____________________________________________________________________________
5233 Double_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)
5236 // Fit methode for the sigma of the pad response function
5239 //We should have at least 3 points
5240 if(nBins <=3) return -4.0;
5242 TLinearFitter fitter(3,"pol2");
5243 fitter.StoreData(kFALSE);
5244 fitter.ClearPoints();
5246 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5247 Float_t entries = 0;
5248 Int_t nbbinwithentries = 0;
5249 for (Int_t i=0; i<nBins; i++){
5251 if(arraye[i] > 15) nbbinwithentries++;
5252 //printf("entries for i %d: %f\n",i,arraye[i]);
5254 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5255 //printf("entries %f\n",entries);
5256 //printf("nbbinwithentries %d\n",nbbinwithentries);
5259 Float_t errorm = 0.0;
5260 Float_t errorn = 0.0;
5261 Float_t error = 0.0;
5264 for (Int_t ibin=0;ibin<nBins; ibin++){
5265 Float_t entriesI = arraye[ibin];
5266 Float_t valueI = arraym[ibin];
5267 Double_t xcenter = 0.0;
5269 if ((entriesI>15) && (valueI>0.0)){
5270 xcenter = xMin+(ibin+0.5)*binWidth;
5275 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5276 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5277 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5280 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5281 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5282 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5284 if(TMath::Abs(error) < 0.000000001) continue;
5285 val = TMath::Log(Float_t(valueI));
5286 fitter.AddPoint(&xcenter,val,error);
5290 if(fDebugLevel > 1){
5292 if ( !fDebugStreamer ) {
5294 TDirectory *backup = gDirectory;
5295 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5296 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5299 Int_t detector = fCountDet;
5300 Int_t layer = GetLayer(fCountDet);
5303 (* fDebugStreamer) << "FitGausMIFill"<<
5304 "detector="<<detector<<
5308 "entriesI="<<entriesI<<
5311 "xcenter="<<xcenter<<
5321 if(npoints <=3) return -4.0;
5326 fitter.GetParameters(par);
5327 chi2 = fitter.GetChisquare()/Float_t(npoints);
5330 if (!param) param = new TVectorD(3);
5331 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5332 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5333 Double_t deltax = (fitter.GetParError(2))/x;
5334 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5337 (*param)[1] = par[1]/(-2.*par[2]);
5338 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5339 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5340 if ( lnparam0>307 ) return -4;
5341 (*param)[0] = TMath::Exp(lnparam0);
5343 if(fDebugLevel > 1){
5345 if ( !fDebugStreamer ) {
5347 TDirectory *backup = gDirectory;
5348 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5349 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5352 Int_t detector = fCountDet;
5353 Int_t layer = GetLayer(fCountDet);
5356 (* fDebugStreamer) << "FitGausMIFit"<<
5357 "detector="<<detector<<
5360 "errorsigma="<<chi2<<
5361 "mean="<<(*param)[1]<<
5362 "sigma="<<(*param)[2]<<
5363 "constant="<<(*param)[0]<<
5368 if((chi2/(*param)[2]) > 0.1){
5370 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5375 if(fDebugLevel == 1){
5376 TString name("PRF");
5377 name += (Int_t)xMin;
5378 name += (Int_t)xMax;
5379 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5382 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5383 for(Int_t k = 0; k < nBins; k++){
5384 histo->SetBinContent(k+1,arraym[k]);
5385 histo->SetBinError(k+1,arrayme[k]);
5388 name += "functionf";
5389 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5390 f1->SetParameter(0, (*param)[0]);
5391 f1->SetParameter(1, (*param)[1]);
5392 f1->SetParameter(2, (*param)[2]);
5400 //_____________________________________________________________________________
5401 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5404 // Fit methode for the sigma of the pad response function
5407 fCurrentCoef[0] = 0.0;
5408 fCurrentCoefE = 0.0;
5410 if (fDebugLevel != 1) {
5411 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5414 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5416 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5419 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5420 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5421 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5423 fNumberFitSuccess++;
5426 //_____________________________________________________________________________
5427 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5430 // Fit methode for the sigma of the pad response function
5432 fCurrentCoef[0] = 0.0;
5433 fCurrentCoefE = 0.0;
5434 if (fDebugLevel == 1) {
5435 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5439 fCurrentCoef[0] = projPRF->GetRMS();
5440 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5442 fNumberFitSuccess++;
5445 //_____________________________________________________________________________
5446 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5449 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5452 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5455 Int_t nbins = (Int_t)(nybins/(2*nbg));
5456 Float_t lowedge = -3.0*nbg;
5457 Float_t upedge = lowedge + 3.0;
5460 Double_t xvalues = -0.2*nbg+0.1;
5462 Int_t total = 2*nbg;
5465 for(Int_t k = 0; k < total; k++){
5466 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5468 y = fCurrentCoef[0]*fCurrentCoef[0];
5469 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5472 if(fDebugLevel > 1){
5474 if ( !fDebugStreamer ) {
5476 TDirectory *backup = gDirectory;
5477 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5478 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5481 Int_t detector = fCountDet;
5482 Int_t layer = GetLayer(fCountDet);
5483 Int_t nbtotal = total;
5485 Float_t low = lowedge;
5486 Float_t up = upedge;
5487 Float_t tnp = xvalues;
5488 Float_t wid = fCurrentCoef[0];
5489 Float_t widfE = fCurrentCoefE;
5491 (* fDebugStreamer) << "FitTnpRange0"<<
5492 "detector="<<detector<<
5494 "nbtotal="<<nbtotal<<
5512 fCurrentCoefE = 0.0;
5513 fCurrentCoef[0] = 0.0;
5515 //printf("npoints\n",npoints);
5518 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5523 linearfitter.Eval();
5524 linearfitter.GetParameters(pars0);
5525 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5526 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5527 Double_t min0 = 0.0;
5528 Double_t ermin0 = 0.0;
5529 //Double_t prfe0 = 0.0;
5530 Double_t prf0 = 0.0;
5531 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5532 min0 = -pars0[1]/(2*pars0[2]);
5533 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5534 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5537 prfe0 = linearfitter->GetParError(0)*pointError0
5538 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5539 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5540 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5541 fCurrentCoefE = (Float_t) prfe0;
5543 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5546 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5550 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5553 if(fDebugLevel > 1){
5555 if ( !fDebugStreamer ) {
5557 TDirectory *backup = gDirectory;
5558 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5559 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5562 Int_t detector = fCountDet;
5563 Int_t layer = GetLayer(fCountDet);
5564 Int_t nbtotal = total;
5565 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5566 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5568 (* fDebugStreamer) << "FitTnpRange1"<<
5569 "detector="<<detector<<
5571 "nbtotal="<<nbtotal<<
5575 "npoints="<<npoints<<
5578 "sigmaprf="<<fCurrentCoef[0]<<
5579 "sigprf="<<fCurrentCoef[1]<<
5586 //_____________________________________________________________________________
5587 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5590 // Only mean methode for the gain factor
5593 fCurrentCoef[0] = mean;
5594 fCurrentCoefE = 0.0;
5595 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5596 if (fDebugLevel == 1) {
5597 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5601 CalculChargeCoefMean(kTRUE);
5602 fNumberFitSuccess++;
5604 //_____________________________________________________________________________
5605 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5608 // mean w methode for the gain factor
5612 Int_t nybins = projch->GetNbinsX();
5614 //The weight function
5615 Double_t a = 0.00228515;
5616 Double_t b = -0.00231487;
5617 Double_t c = 0.00044298;
5618 Double_t d = -0.00379239;
5619 Double_t e = 0.00338349;
5629 //A arbitrary error for the moment
5630 fCurrentCoefE = 0.0;
5631 fCurrentCoef[0] = 0.0;
5634 Double_t sumw = 0.0;
5636 Float_t sumAll = (Float_t) nentries;
5637 Int_t sumCurrent = 0;
5638 for(Int_t k = 0; k <nybins; k++){
5639 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5640 if (fraction>0.95) break;
5641 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5642 e*fraction*fraction*fraction*fraction;
5643 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5644 sum += weight*projch->GetBinContent(k+1);
5645 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5646 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5648 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5650 if (fDebugLevel == 1) {
5651 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5655 fNumberFitSuccess++;
5656 CalculChargeCoefMean(kTRUE);
5658 //_____________________________________________________________________________
5659 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5662 // mean w methode for the gain factor
5666 Int_t nybins = projch->GetNbinsX();
5668 //The weight function
5669 Double_t a = 0.00228515;
5670 Double_t b = -0.00231487;
5671 Double_t c = 0.00044298;
5672 Double_t d = -0.00379239;
5673 Double_t e = 0.00338349;
5683 //A arbitrary error for the moment
5684 fCurrentCoefE = 0.0;
5685 fCurrentCoef[0] = 0.0;
5688 Double_t sumw = 0.0;
5690 Int_t sumCurrent = 0;
5691 for(Int_t k = 0; k <nybins; k++){
5692 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5693 if (fraction>0.95) break;
5694 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5695 e*fraction*fraction*fraction*fraction;
5696 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5697 sum += weight*projch->GetBinContent(k+1);
5698 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5699 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5701 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5703 if (fDebugLevel == 1) {
5704 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5708 fNumberFitSuccess++;
5710 //_____________________________________________________________________________
5711 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5714 // Fit methode for the gain factor
5717 fCurrentCoef[0] = 0.0;
5718 fCurrentCoefE = 0.0;
5719 Double_t chisqrl = 0.0;
5720 Double_t chisqrg = 0.0;
5721 Double_t chisqr = 0.0;
5722 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5724 projch->Fit("landau","0",""
5725 ,(Double_t) mean/fBeginFitCharge
5726 ,projch->GetBinCenter(projch->GetNbinsX()));
5727 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5728 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5729 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5730 chisqrl = projch->GetFunction("landau")->GetChisquare();
5732 projch->Fit("gaus","0",""
5733 ,(Double_t) mean/fBeginFitCharge
5734 ,projch->GetBinCenter(projch->GetNbinsX()));
5735 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5736 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5737 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5739 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5740 if (fDebugLevel != 1) {
5741 projch->Fit("fLandauGaus","0",""
5742 ,(Double_t) mean/fBeginFitCharge
5743 ,projch->GetBinCenter(projch->GetNbinsX()));
5744 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5747 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5749 projch->Fit("fLandauGaus","+",""
5750 ,(Double_t) mean/fBeginFitCharge
5751 ,projch->GetBinCenter(projch->GetNbinsX()));
5752 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5754 fLandauGaus->Draw("same");
5757 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5758 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5759 fNumberFitSuccess++;
5760 CalculChargeCoefMean(kTRUE);
5761 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5762 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5765 CalculChargeCoefMean(kFALSE);
5766 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5769 if (fDebugLevel != 1) {
5774 //_____________________________________________________________________________
5775 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5778 // Fit methode for the gain factor more time consuming
5782 //Some parameters to initialise
5783 Double_t widthLandau, widthGaus, mPV, integral;
5784 Double_t chisquarel = 0.0;
5785 Double_t chisquareg = 0.0;
5786 projch->Fit("landau","0M+",""
5788 ,projch->GetBinCenter(projch->GetNbinsX()));
5789 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5790 chisquarel = projch->GetFunction("landau")->GetChisquare();
5791 projch->Fit("gaus","0M+",""
5793 ,projch->GetBinCenter(projch->GetNbinsX()));
5794 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5795 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5797 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5798 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5800 // Setting fit range and start values
5802 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5803 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5804 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5805 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5806 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5807 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5808 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5811 fCurrentCoef[0] = 0.0;
5812 fCurrentCoefE = 0.0;
5816 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5821 Double_t projchPeak;
5822 Double_t projchFWHM;
5823 LanGauPro(fp,projchPeak,projchFWHM);
5825 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5826 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5827 fNumberFitSuccess++;
5828 CalculChargeCoefMean(kTRUE);
5829 fCurrentCoef[0] = fp[1];
5830 fCurrentCoefE = fpe[1];
5831 //chargeCoefE2 = chisqr;
5834 CalculChargeCoefMean(kFALSE);
5835 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5837 if (fDebugLevel == 1) {
5838 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5839 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5842 fitsnr->Draw("same");
5848 //_____________________________________________________________________________
5849 void 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
5852 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5854 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5855 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5856 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5861 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5862 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5866 //_____________________________________________________________________________
5867 void 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
5870 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5872 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5873 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5874 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5875 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5879 c2 = -(x0*(x[1]+x[2]+x[3])
5880 +x1*(x[0]+x[2]+x[3])
5881 +x2*(x[0]+x[1]+x[3])
5882 +x3*(x[0]+x[1]+x[2]));
5883 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5884 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5885 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5886 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5888 c0 = -(x0*x[1]*x[2]*x[3]
5891 +x3*x[0]*x[1]*x[2]);
5896 //_____________________________________________________________________________
5897 void 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
5900 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5902 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5903 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5904 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5905 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5906 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5909 c4 = x0+x1+x2+x3+x4;
5910 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5911 +x1*(x[0]+x[2]+x[3]+x[4])
5912 +x2*(x[0]+x[1]+x[3]+x[4])
5913 +x3*(x[0]+x[1]+x[2]+x[4])
5914 +x4*(x[0]+x[1]+x[2]+x[3]));
5915 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])
5916 +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])
5917 +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])
5918 +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])
5919 +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]));
5921 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])
5922 +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])
5923 +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])
5924 +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])
5925 +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]));
5927 c0 = (x0*x[1]*x[2]*x[3]*x[4]
5928 +x1*x[0]*x[2]*x[3]*x[4]
5929 +x2*x[0]*x[1]*x[3]*x[4]
5930 +x3*x[0]*x[1]*x[2]*x[4]
5931 +x4*x[0]*x[1]*x[2]*x[3]);
5934 //_____________________________________________________________________________
5935 void AliTRDCalibraFit::NormierungCharge()
5938 // Normalisation of the gain factor resulting for the fits
5941 // Calcul of the mean of choosen method by fFitChargeNDB
5943 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5944 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5946 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5947 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5948 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5949 if (GetStack(detector) == 2) {
5952 if (GetStack(detector) != 2) {
5955 for (Int_t j = 0; j < total; j++) {
5963 fScaleFitFactor = fScaleFitFactor / sum;
5966 fScaleFitFactor = 1.0;
5969 //methode de boeuf mais bon...
5970 Double_t scalefactor = fScaleFitFactor;
5972 if(fDebugLevel > 1){
5974 if ( !fDebugStreamer ) {
5976 TDirectory *backup = gDirectory;
5977 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5978 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5980 (* fDebugStreamer) << "NormierungCharge"<<
5981 "scalefactor="<<scalefactor<<
5985 //_____________________________________________________________________________
5986 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5989 // Rebin of the 1D histo for the gain calibration if needed.
5990 // you have to choose fRebin, divider of fNumberBinCharge
5993 TAxis *xhist = hist->GetXaxis();
5994 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5995 ,xhist->GetBinLowEdge(1)
5996 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5998 AliInfo(Form("fRebin: %d",fRebin));
6000 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6002 for (Int_t ji = i; ji < i+fRebin; ji++) {
6003 sum += hist->GetBinContent(ji);
6006 rehist->SetBinContent(k,sum);
6014 //_____________________________________________________________________________
6015 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
6018 // Rebin of the 1D histo for the gain calibration if needed
6019 // you have to choose fRebin divider of fNumberBinCharge
6022 TAxis *xhist = hist->GetXaxis();
6023 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6024 ,xhist->GetBinLowEdge(1)
6025 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6027 AliInfo(Form("fRebin: %d",fRebin));
6029 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6031 for (Int_t ji = i; ji < i+fRebin; ji++) {
6032 sum += hist->GetBinContent(ji);
6035 rehist->SetBinContent(k,sum);
6043 //____________Some basic geometry function_____________________________________
6046 //_____________________________________________________________________________
6047 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
6050 // Reconstruct the plane number from the detector number
6053 return ((Int_t) (d % 6));
6057 //_____________________________________________________________________________
6058 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
6061 // Reconstruct the stack number from the detector number
6063 const Int_t kNlayer = 6;
6065 return ((Int_t) (d % 30) / kNlayer);
6069 //_____________________________________________________________________________
6070 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6073 // Reconstruct the sector number from the detector number
6077 return ((Int_t) (d / fg));
6082 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6084 //_______________________________________________________________________________
6085 void AliTRDCalibraFit::ResetVectorFit()
6088 // Reset the VectorFits
6091 fVectorFit.SetOwner();
6093 fVectorFit2.SetOwner();
6094 fVectorFit2.Clear();
6098 //____________Private Functions________________________________________________
6101 //_____________________________________________________________________________
6102 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
6105 // Function for the fit
6108 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6110 //PARAMETERS FOR FIT PH
6112 //fAsymmGauss->SetParameter(0,0.113755);
6113 //fAsymmGauss->SetParameter(1,0.350706);
6114 //fAsymmGauss->SetParameter(2,0.0604244);
6115 //fAsymmGauss->SetParameter(3,7.65596);
6116 //fAsymmGauss->SetParameter(4,1.00124);
6117 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6125 Double_t dx = 0.005;
6126 Double_t xs = par[1];
6128 Double_t paras[2] = { 0.0, 0.0 };
6131 if ((xs >= par[1]) &&
6132 (xs < (par[1]+par[2]))) {
6133 //fAsymmGauss->SetParameter(0,par[0]);
6134 //fAsymmGauss->SetParameter(1,xs);
6135 //ss += fAsymmGauss->Eval(xx);
6138 ss += AsymmGauss(&xx,paras);
6140 if ((xs >= (par[1]+par[2])) &&
6141 (xs < (par[1]+par[2]+par[3]))) {
6142 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6143 //fAsymmGauss->SetParameter(1,xs);
6144 //ss += fAsymmGauss->Eval(xx);
6145 paras[0] = par[0]*par[4];
6147 ss += AsymmGauss(&xx,paras);
6156 //_____________________________________________________________________________
6157 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
6160 // Function for the fit
6163 //par[0] = normalization
6171 Double_t par1save = par[1];
6172 //Double_t par2save = par[2];
6173 Double_t par2save = 0.0604244;
6174 //Double_t par3save = par[3];
6175 Double_t par3save = 7.65596;
6176 //Double_t par5save = par[5];
6177 Double_t par5save = 0.870597;
6178 Double_t dx = x[0] - par1save;
6180 Double_t sigma2 = par2save*par2save;
6181 Double_t sqrt2 = TMath::Sqrt(2.0);
6182 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6183 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6184 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6185 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6187 //return par[0]*(exp1+par[4]*exp2);
6188 return par[0] * (exp1 + 1.00124 * exp2);
6192 //_____________________________________________________________________________
6193 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
6196 // Sum Landau + Gaus with identical mean
6199 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6200 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6201 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6202 Double_t val = valLandau + valGaus;
6208 //_____________________________________________________________________________
6209 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6212 // Function for the fit
6215 // par[0]=Width (scale) parameter of Landau density
6216 // par[1]=Most Probable (MP, location) parameter of Landau density
6217 // par[2]=Total area (integral -inf to inf, normalization constant)
6218 // par[3]=Width (sigma) of convoluted Gaussian function
6220 // In the Landau distribution (represented by the CERNLIB approximation),
6221 // the maximum is located at x=-0.22278298 with the location parameter=0.
6222 // This shift is corrected within this function, so that the actual
6223 // maximum is identical to the MP parameter.
6226 // Numeric constants
6227 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6228 Double_t mpshift = -0.22278298; // Landau maximum location
6230 // Control constants
6231 Double_t np = 100.0; // Number of convolution steps
6232 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6244 // MP shift correction
6245 mpc = par[1] - mpshift * par[0];
6247 // Range of convolution integral
6248 xlow = x[0] - sc * par[3];
6249 xupp = x[0] + sc * par[3];
6251 step = (xupp - xlow) / np;
6253 // Convolution integral of Landau and Gaussian by sum
6254 for (i = 1.0; i <= np/2; i++) {
6256 xx = xlow + (i-.5) * step;
6257 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6258 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6260 xx = xupp - (i-.5) * step;
6261 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6262 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6266 return (par[2] * step * sum * invsq2pi / par[3]);
6269 //_____________________________________________________________________________
6270 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6271 , const Double_t *parlimitslo, const Double_t *parlimitshi
6272 , Double_t *fitparams, Double_t *fiterrors
6273 , Double_t *chiSqr, Int_t *ndf) const
6276 // Function for the fit
6280 Char_t funname[100];
6282 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6287 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6288 ffit->SetParameters(startvalues);
6289 ffit->SetParNames("Width","MP","Area","GSigma");
6291 for (i = 0; i < 4; i++) {
6292 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6295 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6297 ffit->GetParameters(fitparams); // Obtain fit parameters
6298 for (i = 0; i < 4; i++) {
6299 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6301 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6302 ndf[0] = ffit->GetNDF(); // Obtain ndf
6304 return (ffit); // Return fit function
6308 //_____________________________________________________________________________
6309 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6312 // Function for the fit
6325 Int_t maxcalls = 10000;
6327 // Search for maximum
6328 p = params[1] - 0.1 * params[0];
6329 step = 0.05 * params[0];
6333 while ((l != lold) && (i < maxcalls)) {
6337 l = LanGauFun(&x,params);
6339 step = -step / 10.0;
6344 if (i == maxcalls) {
6350 // Search for right x location of fy
6351 p = maxx + params[0];
6357 while ( (l != lold) && (i < maxcalls) ) {
6362 l = TMath::Abs(LanGauFun(&x,params) - fy);
6376 // Search for left x location of fy
6378 p = maxx - 0.5 * params[0];
6384 while ((l != lold) && (i < maxcalls)) {
6388 l = TMath::Abs(LanGauFun(&x,params) - fy);
6390 step = -step / 10.0;
6395 if (i == maxcalls) {
6404 //_____________________________________________________________________________
6405 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6408 // Gaus with identical mean
6411 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];