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>
63 #include <TTreeStream.h>
64 #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)
140 ,fNumberFitSuccess(0)
147 ,fCalibraMode(new AliTRDCalibraMode())
152 ,fScaleFitFactor(0.0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
167 // Default constructor
170 fGeo = new AliTRDgeometry();
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
177 for (Int_t i = 0; i < 3; i++) {
183 //______________________________________________________________________________________
184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
189 ,fBeginFitCharge(c.fBeginFitCharge)
190 ,fFitPHPeriode(c.fFitPHPeriode)
191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
192 ,fT0Shift0(c.fT0Shift0)
193 ,fT0Shift1(c.fT0Shift1)
194 ,fRangeFitPRF(c.fRangeFitPRF)
196 ,fMinEntries(c.fMinEntries)
198 ,fNumberFit(c.fNumberFit)
199 ,fNumberFitSuccess(c.fNumberFitSuccess)
200 ,fNumberEnt(c.fNumberEnt)
201 ,fStatisticMean(c.fStatisticMean)
203 ,fDebugLevel(c.fDebugLevel)
204 ,fFitVoir(c.fFitVoir)
205 ,fMagneticField(c.fMagneticField)
207 ,fCurrentCoefE(c.fCurrentCoefE)
208 ,fCurrentCoefE2(c.fCurrentCoefE2)
211 ,fScaleFitFactor(c.fScaleFitFactor)
212 ,fEntriesCurrent(c.fEntriesCurrent)
213 ,fCountDet(c.fCountDet)
220 ,fCurrentCoefDetector(0x0)
221 ,fCurrentCoefDetector2(0x0)
229 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
231 //Current variables initialised
232 for (Int_t k = 0; k < 2; k++) {
233 fCurrentCoef[k] = 0.0;
234 fCurrentCoef2[k] = 0.0;
236 for (Int_t i = 0; i < 3; i++) {
240 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
241 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
243 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
244 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
246 fVectorFit.SetName(c.fVectorFit.GetName());
247 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
248 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
249 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
251 if (GetStack(detector) == 2) {
257 Float_t *coef = new Float_t[ntotal];
258 for (Int_t i = 0; i < ntotal; i++) {
259 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
261 fitInfo->SetCoef(coef);
262 fitInfo->SetDetector(detector);
263 fVectorFit.Add((TObject *) fitInfo);
265 fVectorFit.SetName(c.fVectorFit.GetName());
266 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
267 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
268 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
270 if (GetStack(detector) == 2) {
276 Float_t *coef = new Float_t[ntotal];
277 for (Int_t i = 0; i < ntotal; i++) {
278 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
280 fitInfo->SetCoef(coef);
281 fitInfo->SetDetector(detector);
282 fVectorFit2.Add((TObject *) fitInfo);
287 fGeo = new AliTRDgeometry();
290 //____________________________________________________________________________________
291 AliTRDCalibraFit::~AliTRDCalibraFit()
294 // AliTRDCalibraFit destructor
296 if ( fDebugStreamer ) delete fDebugStreamer;
297 if ( fCalDet ) delete fCalDet;
298 if ( fCalDet2 ) delete fCalDet2;
299 if ( fCalROC ) delete fCalROC;
300 if ( fCalROC2 ) delete fCalROC2;
301 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
302 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
304 fVectorFit2.Delete();
310 //_____________________________________________________________________________
311 void AliTRDCalibraFit::Destroy()
323 //_____________________________________________________________________________
324 void AliTRDCalibraFit::DestroyDebugStreamer()
327 // Delete DebugStreamer
330 if ( fDebugStreamer ) delete fDebugStreamer;
331 fDebugStreamer = 0x0;
334 //__________________________________________________________________________________
335 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
338 // From the drift velocity and t0
339 // return the position of the peak and maximum negative slope
342 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
343 Double_t widbins = 0.1; // 0.1 mus
345 //peak and maxnegslope in mus
346 Double_t begind = t0*widbins + fT0Shift0;
347 Double_t peakd = t0*widbins + fT0Shift1;
348 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
350 // peak and maxnegslope in timebin
351 begin = TMath::Nint(begind*widbins);
352 peak = TMath::Nint(peakd*widbins);
353 end = TMath::Nint(maxnegslope*widbins);
356 //____________Functions fit Online CH2d________________________________________
357 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
360 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
361 // calibration group normalized the resulted coefficients (to 1 normally)
364 // Set the calibration mode
365 const char *name = ch->GetTitle();
366 if(!SetModeCalibration(name,0)) return kFALSE;
368 // Number of Ybins (detectors or groups of pads)
369 Int_t nbins = ch->GetNbinsX();// charge
370 Int_t nybins = ch->GetNbinsY();// groups number
371 if (!InitFit(nybins,0)) {
377 fStatisticMean = 0.0;
379 fNumberFitSuccess = 0;
381 // Init fCountDet and fCount
382 InitfCountDetAndfCount(0);
383 // Beginning of the loop betwwen dect1 and dect2
384 for (Int_t idect = fDect1; idect < fDect2; idect++) {
385 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
386 UpdatefCountDetAndfCount(idect,0);
387 ReconstructFitRowMinRowMax(idect, 0);
389 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
390 projch->SetDirectory(0);
391 // Number of entries for this calibration group
392 Double_t nentries = 0.0;
394 for (Int_t k = 0; k < nbins; k++) {
395 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
396 nentries += ch->GetBinContent(binnb);
397 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
398 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
400 projch->SetEntries(nentries);
401 //printf("The number of entries for the group %d is %f\n",idect,nentries);
406 // Rebin and statistic stuff
408 projch = ReBin((TH1I *) projch);
410 // This detector has not enough statistics or was off
411 if (nentries <= fMinEntries) {
412 NotEnoughStatisticCH(idect);
413 if (fDebugLevel != 1) {
418 // Statistics of the group fitted
419 fStatisticMean += nentries;
424 case 0: FitMeanW((TH1 *) projch, nentries); break;
425 case 1: FitMean((TH1 *) projch, nentries, mean); break;
426 case 2: FitCH((TH1 *) projch, mean); break;
427 case 3: FitBisCH((TH1 *) projch, mean); break;
428 default: return kFALSE;
431 FillInfosFitCH(idect);
433 if (fDebugLevel != 1) {
438 if (fDebugLevel != 1) {
442 if (fNumberFit > 0) {
443 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));
444 fStatisticMean = fStatisticMean / fNumberFit;
447 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
449 delete fDebugStreamer;
450 fDebugStreamer = 0x0;
454 //____________Functions fit Online CH2d________________________________________
455 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
458 // Reconstruct a 1D histo from the vectorCH for each calibration group,
459 // fit the histo, normalized the resulted coefficients (to 1 normally)
462 // Set the calibraMode
463 const char *name = calvect->GetNameCH();
464 if(!SetModeCalibration(name,0)) return kFALSE;
466 // Number of Xbins (detectors or groups of pads)
467 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
473 fStatisticMean = 0.0;
475 fNumberFitSuccess = 0;
477 // Init fCountDet and fCount
478 InitfCountDetAndfCount(0);
479 // Beginning of the loop between dect1 and dect2
480 for (Int_t idect = fDect1; idect < fDect2; idect++) {
481 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
482 UpdatefCountDetAndfCount(idect,0);
483 ReconstructFitRowMinRowMax(idect,0);
485 Double_t nentries = 0.0;
488 Bool_t something = kTRUE;
489 if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
493 projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
494 projch->SetDirectory(0);
495 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
496 nentries += projch->GetBinContent(k+1);
497 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
503 //printf("The number of entries for the group %d is %f\n",idect,nentries);
506 projch = ReBin((TH1F *) projch);
509 // This detector has not enough statistics or was not found in VectorCH
510 if (nentries <= fMinEntries) {
511 NotEnoughStatisticCH(idect);
512 if (fDebugLevel != 1) {
513 if(projch) delete projch;
517 // Statistic of the histos fitted
518 fStatisticMean += nentries;
523 case 0: FitMeanW((TH1 *) projch, nentries); break;
524 case 1: FitMean((TH1 *) projch, nentries, mean); break;
525 case 2: FitCH((TH1 *) projch, mean); break;
526 case 3: FitBisCH((TH1 *) projch, mean); break;
527 default: return kFALSE;
530 FillInfosFitCH(idect);
532 if (fDebugLevel != 1) {
537 if (fDebugLevel != 1) {
541 if (fNumberFit > 0) {
542 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));
543 fStatisticMean = fStatisticMean / fNumberFit;
546 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
548 delete fDebugStreamer;
549 fDebugStreamer = 0x0;
552 //________________functions fit Online PH2d____________________________________
553 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
556 // Take the 1D profiles (average pulse height), projections of the 2D PH
557 // on the Xaxis, for each calibration group
558 // Reconstruct a drift velocity
559 // A first calibration of T0 is also made using the same method
562 // Set the calibration mode
563 const char *name = ph->GetTitle();
564 if(!SetModeCalibration(name,1)) return kFALSE;
566 // Number of Xbins (detectors or groups of pads)
567 Int_t nbins = ph->GetNbinsX();// time
568 Int_t nybins = ph->GetNbinsY();// calibration group
569 if (!InitFit(nybins,1)) {
575 fStatisticMean = 0.0;
577 fNumberFitSuccess = 0;
579 // Init fCountDet and fCount
580 InitfCountDetAndfCount(1);
581 // Beginning of the loop
582 for (Int_t idect = fDect1; idect < fDect2; idect++) {
583 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
584 UpdatefCountDetAndfCount(idect,1);
585 ReconstructFitRowMinRowMax(idect,1);
587 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
588 projph->SetDirectory(0);
589 // Number of entries for this calibration group
590 Double_t nentries = 0;
591 for (Int_t k = 0; k < nbins; k++) {
592 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
593 nentries += ph->GetBinEntries(binnb);
598 //printf("The number of entries for the group %d is %f\n",idect,nentries);
599 // This detector has not enough statistics or was off
600 if (nentries <= fMinEntries) {
601 //printf("Not enough statistic!\n");
602 NotEnoughStatisticPH(idect,nentries);
603 if (fDebugLevel != 1) {
608 // Statistics of the histos fitted
610 fStatisticMean += nentries;
611 // Calcul of "real" coef
612 CalculVdriftCoefMean();
617 case 0: FitLagrangePoly((TH1 *) projph); break;
618 case 1: FitPente((TH1 *) projph); break;
619 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
620 default: return kFALSE;
622 // Fill the tree if end of a detector or only the pointer to the branch!!!
623 FillInfosFitPH(idect,nentries);
625 if (fDebugLevel != 1) {
630 if (fNumberFit > 0) {
631 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));
632 fStatisticMean = fStatisticMean / fNumberFit;
635 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
637 delete fDebugStreamer;
638 fDebugStreamer = 0x0;
641 //____________Functions fit Online PH2d________________________________________
642 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
645 // Reconstruct the average pulse height from the vectorPH for each
647 // Reconstruct a drift velocity
648 // A first calibration of T0 is also made using the same method (slope method)
651 // Set the calibration mode
652 const char *name = calvect->GetNamePH();
653 if(!SetModeCalibration(name,1)) return kFALSE;
655 // Number of Xbins (detectors or groups of pads)
656 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
662 fStatisticMean = 0.0;
664 fNumberFitSuccess = 0;
666 // Init fCountDet and fCount
667 InitfCountDetAndfCount(1);
668 // Beginning of the loop
669 for (Int_t idect = fDect1; idect < fDect2; idect++) {
670 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
671 UpdatefCountDetAndfCount(idect,1);
672 ReconstructFitRowMinRowMax(idect,1);
676 Bool_t something = kTRUE;
677 if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
681 projph = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)));
682 projph->SetDirectory(0);
684 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
685 // This detector has not enough statistics or was off
686 if (fEntriesCurrent <= fMinEntries) {
687 //printf("Not enough stat!\n");
688 NotEnoughStatisticPH(idect,fEntriesCurrent);
689 if (fDebugLevel != 1) {
690 if(projph) delete projph;
694 // Statistic of the histos fitted
696 fStatisticMean += fEntriesCurrent;
697 // Calcul of "real" coef
698 CalculVdriftCoefMean();
703 case 0: FitLagrangePoly((TH1 *) projph); break;
704 case 1: FitPente((TH1 *) projph); break;
705 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
706 default: return kFALSE;
708 // Fill the tree if end of a detector or only the pointer to the branch!!!
709 FillInfosFitPH(idect,fEntriesCurrent);
711 if (fDebugLevel != 1) {
717 if (fNumberFit > 0) {
718 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));
719 fStatisticMean = fStatisticMean / fNumberFit;
722 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
724 delete fDebugStreamer;
725 fDebugStreamer = 0x0;
728 //____________Functions fit Online PRF2d_______________________________________
729 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
732 // Take the 1D profiles (pad response function), projections of the 2D PRF
733 // on the Xaxis, for each calibration group
734 // Fit with a gaussian to reconstruct the sigma of the pad response function
737 // Set the calibration mode
738 const char *name = prf->GetTitle();
739 if(!SetModeCalibration(name,2)) return kFALSE;
741 // Number of Ybins (detectors or groups of pads)
742 Int_t nybins = prf->GetNbinsY();// calibration groups
743 Int_t nbins = prf->GetNbinsX();// bins
744 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
745 if((nbg > 0) || (nbg == -1)) return kFALSE;
746 if (!InitFit(nybins,2)) {
752 fStatisticMean = 0.0;
754 fNumberFitSuccess = 0;
756 // Init fCountDet and fCount
757 InitfCountDetAndfCount(2);
758 // Beginning of the loop
759 for (Int_t idect = fDect1; idect < fDect2; idect++) {
760 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
761 UpdatefCountDetAndfCount(idect,2);
762 ReconstructFitRowMinRowMax(idect,2);
764 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
765 projprf->SetDirectory(0);
766 // Number of entries for this calibration group
767 Double_t nentries = 0;
768 for (Int_t k = 0; k < nbins; k++) {
769 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
770 nentries += prf->GetBinEntries(binnb);
772 if(nentries > 0) fNumberEnt++;
773 // This detector has not enough statistics or was off
774 if (nentries <= fMinEntries) {
775 NotEnoughStatisticPRF(idect);
776 if (fDebugLevel != 1) {
781 // Statistics of the histos fitted
783 fStatisticMean += nentries;
784 // Calcul of "real" coef
789 case 0: FitPRF((TH1 *) projprf); break;
790 case 1: RmsPRF((TH1 *) projprf); break;
791 default: return kFALSE;
793 // Fill the tree if end of a detector or only the pointer to the branch!!!
794 FillInfosFitPRF(idect);
796 if (fDebugLevel != 1) {
801 if (fNumberFit > 0) {
802 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
803 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
804 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
805 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
806 fStatisticMean = fStatisticMean / fNumberFit;
809 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
811 delete fDebugStreamer;
812 fDebugStreamer = 0x0;
815 //____________Functions fit Online PRF2d_______________________________________
816 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
819 // Take the 1D profiles (pad response function), projections of the 2D PRF
820 // on the Xaxis, for each calibration group
821 // Fit with a gaussian to reconstruct the sigma of the pad response function
824 // Set the calibration mode
825 const char *name = prf->GetTitle();
826 if(!SetModeCalibration(name,2)) return kFALSE;
828 // Number of Ybins (detectors or groups of pads)
829 TAxis *xprf = prf->GetXaxis();
830 TAxis *yprf = prf->GetYaxis();
831 Int_t nybins = yprf->GetNbins();// calibration groups
832 Int_t nbins = xprf->GetNbins();// bins
833 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
834 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
835 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
836 if(nbg == -1) return kFALSE;
837 if(nbg > 0) fMethod = 1;
839 if (!InitFit(nybins,2)) {
845 fStatisticMean = 0.0;
847 fNumberFitSuccess = 0;
849 // Init fCountDet and fCount
850 InitfCountDetAndfCount(2);
851 // Beginning of the loop
852 for (Int_t idect = fDect1; idect < fDect2; idect++) {
853 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
854 UpdatefCountDetAndfCount(idect,2);
855 ReconstructFitRowMinRowMax(idect,2);
856 // Build the array of entries and sum
857 TArrayD arraye = TArrayD(nbins);
858 TArrayD arraym = TArrayD(nbins);
859 TArrayD arrayme = TArrayD(nbins);
860 Double_t nentries = 0;
861 //printf("nbins %d\n",nbins);
862 for (Int_t k = 0; k < nbins; k++) {
863 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
864 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
865 Double_t mean = (Double_t)prf->GetBinContent(binnb);
866 Double_t error = (Double_t)prf->GetBinError(binnb);
867 //printf("for %d we have %f\n",k,entries);
869 arraye.AddAt(entries,k);
870 arraym.AddAt(mean,k);
871 arrayme.AddAt(error,k);
873 if(nentries > 0) fNumberEnt++;
874 //printf("The number of entries for the group %d is %f\n",idect,nentries);
875 // This detector has not enough statistics or was off
876 if (nentries <= fMinEntries) {
877 NotEnoughStatisticPRF(idect);
880 // Statistics of the histos fitted
882 fStatisticMean += nentries;
883 // Calcul of "real" coef
888 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
889 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
890 default: return kFALSE;
892 // Fill the tree if end of a detector or only the pointer to the branch!!!
893 FillInfosFitPRF(idect);
896 if (fNumberFit > 0) {
897 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
898 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
899 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
900 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
901 fStatisticMean = fStatisticMean / fNumberFit;
904 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
906 delete fDebugStreamer;
907 fDebugStreamer = 0x0;
910 //____________Functions fit Online PRF2d_______________________________________
911 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
914 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
915 // each calibration group
916 // Fit with a gaussian to reconstruct the sigma of the pad response function
919 // Set the calibra mode
920 const char *name = calvect->GetNamePRF();
921 if(!SetModeCalibration(name,2)) return kFALSE;
922 //printf("test0 %s\n",name);
924 // Number of Xbins (detectors or groups of pads)
925 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
930 ///printf("test2\n");
933 fStatisticMean = 0.0;
935 fNumberFitSuccess = 0;
937 // Init fCountDet and fCount
938 InitfCountDetAndfCount(2);
939 // Beginning of the loop
940 for (Int_t idect = fDect1; idect < fDect2; idect++) {
941 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
942 UpdatefCountDetAndfCount(idect,2);
943 ReconstructFitRowMinRowMax(idect,2);
947 Bool_t something = kTRUE;
948 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
950 TString tname("PRF");
952 projprf = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)));
953 projprf->SetDirectory(0);
955 // This detector has not enough statistics or was off
956 if (fEntriesCurrent <= fMinEntries) {
957 NotEnoughStatisticPRF(idect);
958 if (fDebugLevel != 1) {
959 if(projprf) delete projprf;
963 // Statistic of the histos fitted
965 fStatisticMean += fEntriesCurrent;
966 // Calcul of "real" coef
971 case 1: FitPRF((TH1 *) projprf); break;
972 case 2: RmsPRF((TH1 *) projprf); break;
973 default: return kFALSE;
975 // Fill the tree if end of a detector or only the pointer to the branch!!!
976 FillInfosFitPRF(idect);
978 if (fDebugLevel != 1) {
983 if (fNumberFit > 0) {
984 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));
987 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
989 delete fDebugStreamer;
990 fDebugStreamer = 0x0;
993 //____________Functions fit Online PRF2d_______________________________________
994 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
997 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
998 // each calibration group
999 // Fit with a gaussian to reconstruct the sigma of the pad response function
1002 // Set the calibra mode
1003 const char *name = calvect->GetNamePRF();
1004 if(!SetModeCalibration(name,2)) return kFALSE;
1005 //printf("test0 %s\n",name);
1006 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1007 //printf("test1 %d\n",nbg);
1008 if(nbg == -1) return kFALSE;
1009 if(nbg > 0) fMethod = 1;
1011 // Number of Xbins (detectors or groups of pads)
1012 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1013 //printf("test2\n");
1016 if (!InitFitPRF()) {
1017 //printf("test3\n");
1020 fStatisticMean = 0.0;
1022 fNumberFitSuccess = 0;
1026 Double_t *arrayx = 0;
1027 Double_t *arraye = 0;
1028 Double_t *arraym = 0;
1029 Double_t *arrayme = 0;
1030 Float_t lowedge = 0.0;
1031 Float_t upedge = 0.0;
1032 // Init fCountDet and fCount
1033 InitfCountDetAndfCount(2);
1034 // Beginning of the loop
1035 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1036 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1037 UpdatefCountDetAndfCount(idect,2);
1038 ReconstructFitRowMinRowMax(idect,2);
1040 TGraphErrors *projprftree = 0x0;
1041 fEntriesCurrent = 0;
1042 Bool_t something = kTRUE;
1043 if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1045 TString tname("PRF");
1047 projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1048 nbins = projprftree->GetN();
1049 arrayx = (Double_t *)projprftree->GetX();
1050 arraye = (Double_t *)projprftree->GetEX();
1051 arraym = (Double_t *)projprftree->GetY();
1052 arrayme = (Double_t *)projprftree->GetEY();
1053 Float_t step = arrayx[1]-arrayx[0];
1054 lowedge = arrayx[0] - step/2.0;
1055 upedge = arrayx[(nbins-1)] + step/2.0;
1056 //printf("nbins est %d\n",nbins);
1057 for(Int_t k = 0; k < nbins; k++){
1058 fEntriesCurrent += (Int_t)arraye[k];
1059 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1060 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1062 if(fEntriesCurrent > 0) fNumberEnt++;
1064 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1065 // This detector has not enough statistics or was off
1066 if (fEntriesCurrent <= fMinEntries) {
1067 NotEnoughStatisticPRF(idect);
1068 if(projprftree) delete projprftree;
1071 // Statistic of the histos fitted
1073 fStatisticMean += fEntriesCurrent;
1074 // Calcul of "real" coef
1075 CalculPRFCoefMean();
1079 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1080 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1081 default: return kFALSE;
1083 // Fill the tree if end of a detector or only the pointer to the branch!!!
1084 FillInfosFitPRF(idect);
1086 if (fDebugLevel != 1) {
1091 if (fNumberFit > 0) {
1092 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));
1095 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1097 delete fDebugStreamer;
1098 fDebugStreamer = 0x0;
1101 //____________Functions fit Online CH2d________________________________________
1102 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1105 // The linear method
1108 fStatisticMean = 0.0;
1110 fNumberFitSuccess = 0;
1112 if(!InitFitLinearFitter()) return kFALSE;
1115 for(Int_t idet = 0; idet < 540; idet++){
1118 //printf("detector number %d\n",idet);
1123 fEntriesCurrent = 0;
1125 Bool_t here = calivdli->GetParam(idet,¶m);
1126 Bool_t heree = calivdli->GetError(idet,&error);
1127 //printf("here %d and heree %d\n",here, heree);
1129 fEntriesCurrent = (Int_t) error[2];
1132 //printf("Number of entries %d\n",fEntriesCurrent);
1133 // Nothing found or not enough statistic
1134 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1135 NotEnoughStatisticLinearFitter();
1142 fStatisticMean += fEntriesCurrent;
1145 if((-(param[1])) <= 0.0) {
1146 NotEnoughStatisticLinearFitter();
1150 // CalculDatabaseVdriftandTan
1151 CalculVdriftLorentzCoef();
1154 fNumberFitSuccess ++;
1156 // Put the fCurrentCoef
1157 fCurrentCoef[0] = -param[1];
1158 // here the database must be the one of the reconstruction for the lorentz angle....
1159 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1160 fCurrentCoefE = error[1];
1161 fCurrentCoefE2 = error[0];
1162 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1163 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1167 FillInfosFitLinearFitter();
1172 if (fNumberFit > 0) {
1173 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));
1176 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1178 delete fDebugStreamer;
1179 fDebugStreamer = 0x0;
1183 //____________Functions for seeing if the pad is really okey___________________
1184 //_____________________________________________________________________________
1185 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1188 // Get numberofgroupsprf
1192 const Char_t *pattern0 = "Ngp0";
1193 const Char_t *pattern1 = "Ngp1";
1194 const Char_t *pattern2 = "Ngp2";
1195 const Char_t *pattern3 = "Ngp3";
1196 const Char_t *pattern4 = "Ngp4";
1197 const Char_t *pattern5 = "Ngp5";
1198 const Char_t *pattern6 = "Ngp6";
1201 if (strstr(nametitle,pattern0)) {
1204 if (strstr(nametitle,pattern1)) {
1207 if (strstr(nametitle,pattern2)) {
1210 if (strstr(nametitle,pattern3)) {
1213 if (strstr(nametitle,pattern4)) {
1216 if (strstr(nametitle,pattern5)) {
1219 if (strstr(nametitle,pattern6)){
1226 //_____________________________________________________________________________
1227 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1230 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1231 // corresponding to the given name
1234 if(!SetNzFromTObject(name,i)) return kFALSE;
1235 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1240 //_____________________________________________________________________________
1241 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1244 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1245 // corresponding to the given TObject
1249 const Char_t *patternrphi0 = "Nrphi0";
1250 const Char_t *patternrphi1 = "Nrphi1";
1251 const Char_t *patternrphi2 = "Nrphi2";
1252 const Char_t *patternrphi3 = "Nrphi3";
1253 const Char_t *patternrphi4 = "Nrphi4";
1254 const Char_t *patternrphi5 = "Nrphi5";
1255 const Char_t *patternrphi6 = "Nrphi6";
1258 const Char_t *patternrphi10 = "Nrphi10";
1259 const Char_t *patternrphi100 = "Nrphi100";
1260 const Char_t *patternz10 = "Nz10";
1261 const Char_t *patternz100 = "Nz100";
1264 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1265 fCalibraMode->SetAllTogether(i);
1267 if (fDebugLevel > 1) {
1268 AliInfo(Form("fNbDet %d and 100",fNbDet));
1272 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1273 fCalibraMode->SetPerSuperModule(i);
1275 if (fDebugLevel > 1) {
1276 AliInfo(Form("fNDet %d and 100",fNbDet));
1281 if (strstr(name,patternrphi0)) {
1282 fCalibraMode->SetNrphi(i ,0);
1283 if (fDebugLevel > 1) {
1284 AliInfo(Form("fNbDet %d and 0",fNbDet));
1288 if (strstr(name,patternrphi1)) {
1289 fCalibraMode->SetNrphi(i, 1);
1290 if (fDebugLevel > 1) {
1291 AliInfo(Form("fNbDet %d and 1",fNbDet));
1295 if (strstr(name,patternrphi2)) {
1296 fCalibraMode->SetNrphi(i, 2);
1297 if (fDebugLevel > 1) {
1298 AliInfo(Form("fNbDet %d and 2",fNbDet));
1302 if (strstr(name,patternrphi3)) {
1303 fCalibraMode->SetNrphi(i, 3);
1304 if (fDebugLevel > 1) {
1305 AliInfo(Form("fNbDet %d and 3",fNbDet));
1309 if (strstr(name,patternrphi4)) {
1310 fCalibraMode->SetNrphi(i, 4);
1311 if (fDebugLevel > 1) {
1312 AliInfo(Form("fNbDet %d and 4",fNbDet));
1316 if (strstr(name,patternrphi5)) {
1317 fCalibraMode->SetNrphi(i, 5);
1318 if (fDebugLevel > 1) {
1319 AliInfo(Form("fNbDet %d and 5",fNbDet));
1323 if (strstr(name,patternrphi6)) {
1324 fCalibraMode->SetNrphi(i, 6);
1325 if (fDebugLevel > 1) {
1326 AliInfo(Form("fNbDet %d and 6",fNbDet));
1331 if (fDebugLevel > 1) {
1332 AliInfo(Form("fNbDet %d and rest",fNbDet));
1334 fCalibraMode->SetNrphi(i ,0);
1338 //_____________________________________________________________________________
1339 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1342 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1343 // corresponding to the given TObject
1347 const Char_t *patternz0 = "Nz0";
1348 const Char_t *patternz1 = "Nz1";
1349 const Char_t *patternz2 = "Nz2";
1350 const Char_t *patternz3 = "Nz3";
1351 const Char_t *patternz4 = "Nz4";
1353 const Char_t *patternrphi10 = "Nrphi10";
1354 const Char_t *patternrphi100 = "Nrphi100";
1355 const Char_t *patternz10 = "Nz10";
1356 const Char_t *patternz100 = "Nz100";
1358 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1359 fCalibraMode->SetAllTogether(i);
1361 if (fDebugLevel > 1) {
1362 AliInfo(Form("fNbDet %d and 100",fNbDet));
1366 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1367 fCalibraMode->SetPerSuperModule(i);
1369 if (fDebugLevel > 1) {
1370 AliInfo(Form("fNbDet %d and 10",fNbDet));
1374 if (strstr(name,patternz0)) {
1375 fCalibraMode->SetNz(i, 0);
1376 if (fDebugLevel > 1) {
1377 AliInfo(Form("fNbDet %d and 0",fNbDet));
1381 if (strstr(name,patternz1)) {
1382 fCalibraMode->SetNz(i ,1);
1383 if (fDebugLevel > 1) {
1384 AliInfo(Form("fNbDet %d and 1",fNbDet));
1388 if (strstr(name,patternz2)) {
1389 fCalibraMode->SetNz(i ,2);
1390 if (fDebugLevel > 1) {
1391 AliInfo(Form("fNbDet %d and 2",fNbDet));
1395 if (strstr(name,patternz3)) {
1396 fCalibraMode->SetNz(i ,3);
1397 if (fDebugLevel > 1) {
1398 AliInfo(Form("fNbDet %d and 3",fNbDet));
1402 if (strstr(name,patternz4)) {
1403 fCalibraMode->SetNz(i ,4);
1404 if (fDebugLevel > 1) {
1405 AliInfo(Form("fNbDet %d and 4",fNbDet));
1410 if (fDebugLevel > 1) {
1411 AliInfo(Form("fNbDet %d and rest",fNbDet));
1413 fCalibraMode->SetNz(i ,0);
1416 //______________________________________________________________________
1417 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1419 // ofwhat is equaled to 0: mean value of all passing detectors
1420 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1423 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1425 AliInfo("The Vector Fit is not complete!");
1428 Int_t detector = -1;
1430 Float_t value = 0.0;
1432 /////////////////////////////////
1433 // Calculate the mean values
1434 ////////////////////////////////
1436 ////////////////////////
1437 Double_t meanAll = 0.0;
1438 Double_t meanSupermodule[18];
1439 Double_t meanDetector[540];
1441 Int_t countSupermodule[18];
1442 Int_t countDetector[540];
1443 for(Int_t sm = 0; sm < 18; sm++){
1444 meanSupermodule[sm] = 0.0;
1445 countSupermodule[sm] = 0;
1447 for(Int_t det = 0; det < 540; det++){
1448 meanDetector[det] = 0.0;
1449 countDetector[det] = 0;
1453 for (Int_t k = 0; k < loop; k++) {
1454 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1455 sector = GetSector(detector);
1457 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1459 meanDetector[detector] += value;
1460 countDetector[detector]++;
1461 meanSupermodule[sector] += value;
1462 countSupermodule[sector]++;
1468 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1469 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1470 for (Int_t row = 0; row < rowMax; row++) {
1471 for (Int_t col = 0; col < colMax; col++) {
1472 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1474 meanDetector[detector] += value;
1475 countDetector[detector]++;
1476 meanSupermodule[sector] += value;
1477 countSupermodule[sector]++;
1486 if(countAll > 0) meanAll = meanAll/countAll;
1487 for(Int_t sm = 0; sm < 18; sm++){
1488 if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1490 for(Int_t det = 0; det < 540; det++){
1491 if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1493 // Put the mean value for the no-fitted
1494 /////////////////////////////////////////////
1495 for (Int_t k = 0; k < loop; k++) {
1496 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1497 sector = GetSector(detector);
1498 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1499 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1500 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1502 for (Int_t row = 0; row < rowMax; row++) {
1503 for (Int_t col = 0; col < colMax; col++) {
1504 value = coef[(Int_t)(col*rowMax+row)];
1506 if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1508 if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1509 else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1510 else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1514 if(fDebugLevel > 1){
1516 if ( !fDebugStreamer ) {
1518 TDirectory *backup = gDirectory;
1519 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1520 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1523 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1525 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1526 "detector="<<detector<<
1539 //______________________________________________________________________
1540 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1542 // ofwhat is equaled to 0: mean value of all passing detectors
1543 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1546 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1548 AliInfo("The Vector Fit is not complete!");
1551 Int_t detector = -1;
1553 Float_t value = 0.0;
1555 /////////////////////////////////
1556 // Calculate the mean values
1557 ////////////////////////////////
1559 ////////////////////////
1560 Double_t meanAll = 0.0;
1561 Double_t meanSupermodule[18];
1562 Double_t meanDetector[540];
1564 Int_t countSupermodule[18];
1565 Int_t countDetector[540];
1566 for(Int_t sm = 0; sm < 18; sm++){
1567 meanSupermodule[sm] = 0.0;
1568 countSupermodule[sm] = 0;
1570 for(Int_t det = 0; det < 540; det++){
1571 meanDetector[det] = 0.0;
1572 countDetector[det] = 0;
1576 for (Int_t k = 0; k < loop; k++) {
1577 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1578 sector = GetSector(detector);
1580 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1582 meanDetector[detector] += value;
1583 countDetector[detector]++;
1584 meanSupermodule[sector] += value;
1585 countSupermodule[sector]++;
1591 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1592 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1593 for (Int_t row = 0; row < rowMax; row++) {
1594 for (Int_t col = 0; col < colMax; col++) {
1595 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1597 meanDetector[detector] += value;
1598 countDetector[detector]++;
1599 meanSupermodule[sector] += value;
1600 countSupermodule[sector]++;
1609 if(countAll > 0) meanAll = meanAll/countAll;
1610 for(Int_t sm = 0; sm < 18; sm++){
1611 if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1613 for(Int_t det = 0; det < 540; det++){
1614 if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1616 // Put the mean value for the no-fitted
1617 /////////////////////////////////////////////
1618 for (Int_t k = 0; k < loop; k++) {
1619 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1620 sector = GetSector(detector);
1621 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1622 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1623 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1625 for (Int_t row = 0; row < rowMax; row++) {
1626 for (Int_t col = 0; col < colMax; col++) {
1627 value = coef[(Int_t)(col*rowMax+row)];
1629 if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1631 if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1632 else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1633 else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1637 if(fDebugLevel > 1){
1639 if ( !fDebugStreamer ) {
1641 TDirectory *backup = gDirectory;
1642 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1643 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1646 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1648 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1649 "detector="<<detector<<
1662 //_____________________________________________________________________________
1663 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1666 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1667 // It takes the mean value of the coefficients per detector
1668 // This object has to be written in the database
1671 // Create the DetObject
1672 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1674 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1675 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1676 Int_t detector = -1;
1677 Float_t value = 0.0;
1680 for (Int_t k = 0; k < loop; k++) {
1681 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1684 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1688 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1689 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1690 for (Int_t row = 0; row < rowMax; row++) {
1691 for (Int_t col = 0; col < colMax; col++) {
1692 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1693 mean += TMath::Abs(value);
1697 if(count > 0) mean = mean/count;
1699 object->SetValue(detector,mean);
1704 //_____________________________________________________________________________
1705 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1708 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1709 // It takes the mean value of the coefficients per detector
1710 // This object has to be written in the database
1713 // Create the DetObject
1714 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1717 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1718 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1719 Int_t detector = -1;
1720 Float_t value = 0.0;
1722 for (Int_t k = 0; k < loop; k++) {
1723 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1726 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1727 if(!meanOtherBefore){
1728 if(value > 0) value = value*scaleFitFactor;
1730 else value = value*scaleFitFactor;
1731 mean = TMath::Abs(value);
1735 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1736 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1737 for (Int_t row = 0; row < rowMax; row++) {
1738 for (Int_t col = 0; col < colMax; col++) {
1739 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1740 if(!meanOtherBefore) {
1741 if(value > 0) value = value*scaleFitFactor;
1743 else value = value*scaleFitFactor;
1744 mean += TMath::Abs(value);
1748 if(count > 0) mean = mean/count;
1750 object->SetValue(detector,mean);
1755 //_____________________________________________________________________________
1756 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1759 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1760 // It takes the min value of the coefficients per detector
1761 // This object has to be written in the database
1764 // Create the DetObject
1765 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1767 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1768 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1769 Int_t detector = -1;
1770 Float_t value = 0.0;
1772 for (Int_t k = 0; k < loop; k++) {
1773 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1774 Float_t min = 100.0;
1776 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1778 if(value > 70.0) value = value-100.0;
1783 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1784 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1785 for (Int_t row = 0; row < rowMax; row++) {
1786 for (Int_t col = 0; col < colMax; col++) {
1787 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1789 if(value > 70.0) value = value-100.0;
1791 if(min > value) min = value;
1795 object->SetValue(detector,min);
1801 //_____________________________________________________________________________
1802 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
1805 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1806 // It takes the min value of the coefficients per detector
1807 // This object has to be written in the database
1810 // Create the DetObject
1811 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1814 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1815 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1816 Int_t detector = -1;
1817 Float_t value = 0.0;
1819 for (Int_t k = 0; k < loop; k++) {
1820 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1822 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1823 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1824 Float_t min = 100.0;
1825 for (Int_t row = 0; row < rowMax; row++) {
1826 for (Int_t col = 0; col < colMax; col++) {
1827 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1828 mean += -TMath::Abs(value);
1832 if(count > 0) mean = mean/count;
1834 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1835 object->SetValue(detector,-TMath::Abs(value));
1841 //_____________________________________________________________________________
1842 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
1845 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1846 // You need first to create the object for the detectors,
1847 // where the mean value is put.
1848 // This object has to be written in the database
1851 // Create the DetObject
1852 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1855 for(Int_t k = 0; k < 540; k++){
1856 AliTRDCalROC *calROC = object->GetCalROC(k);
1857 Int_t nchannels = calROC->GetNchannels();
1858 for(Int_t ch = 0; ch < nchannels; ch++){
1859 calROC->SetValue(ch,1.0);
1865 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1866 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1867 Int_t detector = -1;
1868 Float_t value = 0.0;
1870 for (Int_t k = 0; k < loop; k++) {
1871 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1872 AliTRDCalROC *calROC = object->GetCalROC(detector);
1873 Float_t mean = detobject->GetValue(detector);
1874 if(mean == 0) continue;
1875 Int_t rowMax = calROC->GetNrows();
1876 Int_t colMax = calROC->GetNcols();
1877 for (Int_t row = 0; row < rowMax; row++) {
1878 for (Int_t col = 0; col < colMax; col++) {
1879 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1880 if(value > 0) value = value*scaleFitFactor;
1881 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1889 //_____________________________________________________________________________
1890 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1893 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1894 // You need first to create the object for the detectors,
1895 // where the mean value is put.
1896 // This object has to be written in the database
1899 // Create the DetObject
1900 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1903 for(Int_t k = 0; k < 540; k++){
1904 AliTRDCalROC *calROC = object->GetCalROC(k);
1905 Int_t nchannels = calROC->GetNchannels();
1906 for(Int_t ch = 0; ch < nchannels; ch++){
1907 calROC->SetValue(ch,1.0);
1913 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1914 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1915 Int_t detector = -1;
1916 Float_t value = 0.0;
1918 for (Int_t k = 0; k < loop; k++) {
1919 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1920 AliTRDCalROC *calROC = object->GetCalROC(detector);
1921 Float_t mean = detobject->GetValue(detector);
1922 if(mean == 0) continue;
1923 Int_t rowMax = calROC->GetNrows();
1924 Int_t colMax = calROC->GetNcols();
1925 for (Int_t row = 0; row < rowMax; row++) {
1926 for (Int_t col = 0; col < colMax; col++) {
1927 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1928 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1936 //_____________________________________________________________________________
1937 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1940 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1941 // You need first to create the object for the detectors,
1942 // where the mean value is put.
1943 // This object has to be written in the database
1946 // Create the DetObject
1947 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1950 for(Int_t k = 0; k < 540; k++){
1951 AliTRDCalROC *calROC = object->GetCalROC(k);
1952 Int_t nchannels = calROC->GetNchannels();
1953 for(Int_t ch = 0; ch < nchannels; ch++){
1954 calROC->SetValue(ch,0.0);
1960 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1961 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1962 Int_t detector = -1;
1963 Float_t value = 0.0;
1965 for (Int_t k = 0; k < loop; k++) {
1966 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1967 AliTRDCalROC *calROC = object->GetCalROC(detector);
1968 Float_t min = detobject->GetValue(detector);
1969 Int_t rowMax = calROC->GetNrows();
1970 Int_t colMax = calROC->GetNcols();
1971 for (Int_t row = 0; row < rowMax; row++) {
1972 for (Int_t col = 0; col < colMax; col++) {
1973 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1975 if(value > 70.0) value = value - 100.0;
1977 calROC->SetValue(col,row,value-min);
1985 //_____________________________________________________________________________
1986 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
1989 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1990 // This object has to be written in the database
1993 // Create the DetObject
1994 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1996 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1997 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1998 Int_t detector = -1;
1999 Float_t value = 0.0;
2001 for (Int_t k = 0; k < loop; k++) {
2002 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2003 AliTRDCalROC *calROC = object->GetCalROC(detector);
2004 Int_t rowMax = calROC->GetNrows();
2005 Int_t colMax = calROC->GetNcols();
2006 for (Int_t row = 0; row < rowMax; row++) {
2007 for (Int_t col = 0; col < colMax; col++) {
2008 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2009 calROC->SetValue(col,row,TMath::Abs(value));
2017 //_____________________________________________________________________________
2018 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2021 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2022 // 0 successful fit 1 not successful fit
2023 // mean is the mean value over the successful fit
2024 // do not use it for t0: no meaning
2027 // Create the CalObject
2028 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2032 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2034 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2035 for(Int_t k = 0; k < 540; k++){
2036 object->SetValue(k,1.0);
2039 Int_t detector = -1;
2040 Float_t value = 0.0;
2042 for (Int_t k = 0; k < loop; k++) {
2043 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2044 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2045 if(value <= 0) object->SetValue(detector,1.0);
2047 object->SetValue(detector,0.0);
2052 if(count > 0) mean /= count;
2055 //_____________________________________________________________________________
2056 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2059 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2060 // 0 not successful fit 1 successful fit
2061 // mean mean value over the successful fit
2064 // Create the CalObject
2065 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2069 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2071 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2072 for(Int_t k = 0; k < 540; k++){
2073 AliTRDCalROC *calROC = object->GetCalROC(k);
2074 Int_t nchannels = calROC->GetNchannels();
2075 for(Int_t ch = 0; ch < nchannels; ch++){
2076 calROC->SetValue(ch,1.0);
2080 Int_t detector = -1;
2081 Float_t value = 0.0;
2083 for (Int_t k = 0; k < loop; k++) {
2084 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2085 AliTRDCalROC *calROC = object->GetCalROC(detector);
2086 Int_t nchannels = calROC->GetNchannels();
2087 for (Int_t ch = 0; ch < nchannels; ch++) {
2088 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2089 if(value <= 0) calROC->SetValue(ch,1.0);
2091 calROC->SetValue(ch,0.0);
2097 if(count > 0) mean /= count;
2100 //_____________________________________________________________________________
2101 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2104 // Set FitPH if 1 then each detector will be fitted
2107 if (periodeFitPH > 0) {
2108 fFitPHPeriode = periodeFitPH;
2111 AliInfo("periodeFitPH must be higher than 0!");
2115 //_____________________________________________________________________________
2116 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2119 // The fit of the deposited charge distribution begins at
2120 // histo->Mean()/beginFitCharge
2121 // You can here set beginFitCharge
2124 if (beginFitCharge > 0) {
2125 fBeginFitCharge = beginFitCharge;
2128 AliInfo("beginFitCharge must be strict positif!");
2133 //_____________________________________________________________________________
2134 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2137 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2138 // You can here set t0Shift0
2142 fT0Shift0 = t0Shift;
2145 AliInfo("t0Shift0 must be strict positif!");
2150 //_____________________________________________________________________________
2151 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2154 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2155 // You can here set t0Shift1
2159 fT0Shift1 = t0Shift;
2162 AliInfo("t0Shift must be strict positif!");
2167 //_____________________________________________________________________________
2168 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2171 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2172 // You can here set rangeFitPRF
2175 if ((rangeFitPRF > 0) &&
2176 (rangeFitPRF <= 1.5)) {
2177 fRangeFitPRF = rangeFitPRF;
2180 AliInfo("rangeFitPRF must be between 0 and 1.0");
2185 //_____________________________________________________________________________
2186 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2189 // Minimum entries for fitting
2192 if (minEntries > 0) {
2193 fMinEntries = minEntries;
2196 AliInfo("fMinEntries must be >= 0.");
2201 //_____________________________________________________________________________
2202 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2205 // Rebin with rebin time less bins the Ch histo
2206 // You can set here rebin that should divide the number of bins of CH histo
2211 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2214 AliInfo("You have to choose a positiv value!");
2218 //_____________________________________________________________________________
2219 Bool_t AliTRDCalibraFit::FillVectorFit()
2222 // For the Fit functions fill the vector Fit
2225 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2228 if (GetStack(fCountDet) == 2) {
2235 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2236 Float_t *coef = new Float_t[ntotal];
2237 for (Int_t i = 0; i < ntotal; i++) {
2238 coef[i] = fCurrentCoefDetector[i];
2241 Int_t detector = fCountDet;
2243 fitInfo->SetCoef(coef);
2244 fitInfo->SetDetector(detector);
2245 fVectorFit.Add((TObject *) fitInfo);
2250 //_____________________________________________________________________________
2251 Bool_t AliTRDCalibraFit::FillVectorFit2()
2254 // For the Fit functions fill the vector Fit
2257 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2260 if (GetStack(fCountDet) == 2) {
2267 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2268 Float_t *coef = new Float_t[ntotal];
2269 for (Int_t i = 0; i < ntotal; i++) {
2270 coef[i] = fCurrentCoefDetector2[i];
2273 Int_t detector = fCountDet;
2275 fitInfo->SetCoef(coef);
2276 fitInfo->SetDetector(detector);
2277 fVectorFit2.Add((TObject *) fitInfo);
2282 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2283 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2286 // Init the number of expected bins and fDect1[i] fDect2[i]
2289 gStyle->SetPalette(1);
2290 gStyle->SetOptStat(1111);
2291 gStyle->SetPadBorderMode(0);
2292 gStyle->SetCanvasColor(10);
2293 gStyle->SetPadLeftMargin(0.13);
2294 gStyle->SetPadRightMargin(0.01);
2296 // Mode groups of pads: the total number of bins!
2297 CalculNumberOfBinsExpected(i);
2299 // Quick verification that we have the good pad calibration mode!
2300 if (fNumberOfBinsExpected != nbins) {
2301 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2305 // Security for fDebug 3 and 4
2306 if ((fDebugLevel >= 3) &&
2310 AliInfo("This detector doesn't exit!");
2314 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2315 CalculDect1Dect2(i);
2320 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2321 Bool_t AliTRDCalibraFit::InitFitCH()
2324 // Init the fVectorFitCH for normalisation
2325 // Init the histo for debugging
2330 fScaleFitFactor = 0.0;
2331 fCurrentCoefDetector = new Float_t[2304];
2332 for (Int_t k = 0; k < 2304; k++) {
2333 fCurrentCoefDetector[k] = 0.0;
2335 fVectorFit.SetName("gainfactorscoefficients");
2337 // fDebug == 0 nothing
2338 // fDebug == 1 and fFitVoir no histo
2339 if (fDebugLevel == 1) {
2340 if(!CheckFitVoir()) return kFALSE;
2342 //Get the CalDet object
2344 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2346 AliInfo("Could not get calibDB");
2349 if(fCalDet) delete fCalDet;
2350 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2353 Float_t devalue = 1.0;
2354 if(fCalDet) delete fCalDet;
2355 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2356 for(Int_t k = 0; k < 540; k++){
2357 fCalDet->SetValue(k,devalue);
2363 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2364 Bool_t AliTRDCalibraFit::InitFitPH()
2367 // Init the arrays of results
2368 // Init the histos for debugging
2372 fVectorFit.SetName("driftvelocitycoefficients");
2373 fVectorFit2.SetName("t0coefficients");
2375 fCurrentCoefDetector = new Float_t[2304];
2376 for (Int_t k = 0; k < 2304; k++) {
2377 fCurrentCoefDetector[k] = 0.0;
2380 fCurrentCoefDetector2 = new Float_t[2304];
2381 for (Int_t k = 0; k < 2304; k++) {
2382 fCurrentCoefDetector2[k] = 0.0;
2385 //fDebug == 0 nothing
2386 // fDebug == 1 and fFitVoir no histo
2387 if (fDebugLevel == 1) {
2388 if(!CheckFitVoir()) return kFALSE;
2390 //Get the CalDet object
2392 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2394 AliInfo("Could not get calibDB");
2397 if(fCalDet) delete fCalDet;
2398 if(fCalDet2) delete fCalDet2;
2399 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2400 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2403 Float_t devalue = 1.5;
2404 Float_t devalue2 = 0.0;
2405 if(fCalDet) delete fCalDet;
2406 if(fCalDet2) delete fCalDet2;
2407 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2408 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2409 for(Int_t k = 0; k < 540; k++){
2410 fCalDet->SetValue(k,devalue);
2411 fCalDet2->SetValue(k,devalue2);
2416 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2417 Bool_t AliTRDCalibraFit::InitFitPRF()
2420 // Init the calibration mode (Nz, Nrphi), the histograms for
2421 // debugging the fit methods if fDebug > 0,
2425 fVectorFit.SetName("prfwidthcoefficients");
2427 fCurrentCoefDetector = new Float_t[2304];
2428 for (Int_t k = 0; k < 2304; k++) {
2429 fCurrentCoefDetector[k] = 0.0;
2432 // fDebug == 0 nothing
2433 // fDebug == 1 and fFitVoir no histo
2434 if (fDebugLevel == 1) {
2435 if(!CheckFitVoir()) return kFALSE;
2439 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2440 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2443 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2448 fCurrentCoefDetector = new Float_t[2304];
2449 fCurrentCoefDetector2 = new Float_t[2304];
2450 for (Int_t k = 0; k < 2304; k++) {
2451 fCurrentCoefDetector[k] = 0.0;
2452 fCurrentCoefDetector2[k] = 0.0;
2455 //printf("test0\n");
2457 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2459 AliInfo("Could not get calibDB");
2463 //Get the CalDet object
2465 if(fCalDet) delete fCalDet;
2466 if(fCalDet2) delete fCalDet2;
2467 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2468 //printf("test1\n");
2469 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2470 //printf("test2\n");
2471 for(Int_t k = 0; k < 540; k++){
2472 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2474 //printf("test3\n");
2477 Float_t devalue = 1.5;
2478 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2479 if(fCalDet) delete fCalDet;
2480 if(fCalDet2) delete fCalDet2;
2481 //printf("test1\n");
2482 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2483 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2484 //printf("test2\n");
2485 for(Int_t k = 0; k < 540; k++){
2486 fCalDet->SetValue(k,devalue);
2487 fCalDet2->SetValue(k,devalue2);
2489 //printf("test3\n");
2494 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2495 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2498 // Init the current detector where we are fCountDet and the
2499 // next fCount for the functions Fit...
2502 // Loop on the Xbins of ch!!
2503 fCountDet = -1; // Current detector
2504 fCount = 0; // To find the next detector
2507 if (fDebugLevel >= 3) {
2508 // Set countdet to the detector
2509 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2510 // Set counter to write at the end of the detector
2512 // Get the right calib objects
2515 if(fDebugLevel == 1) {
2517 fCalibraMode->CalculXBins(fCountDet,i);
2518 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2520 fCalibraMode->CalculXBins(fCountDet,i);
2522 fCount = fCalibraMode->GetXbins(i);
2524 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2525 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2526 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2527 ,(Int_t) GetStack(fCountDet)
2528 ,(Int_t) GetSector(fCountDet),i);
2531 //_______________________________________________________________________________
2532 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2535 // Calculate the number of bins expected (calibration groups)
2538 fNumberOfBinsExpected = 0;
2540 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2541 fNumberOfBinsExpected = 1;
2545 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2546 fNumberOfBinsExpected = 18;
2550 fCalibraMode->ModePadCalibration(2,i);
2551 fCalibraMode->ModePadFragmentation(0,2,0,i);
2552 fCalibraMode->SetDetChamb2(i);
2553 if (fDebugLevel > 1) {
2554 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2556 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2557 fCalibraMode->ModePadCalibration(0,i);
2558 fCalibraMode->ModePadFragmentation(0,0,0,i);
2559 fCalibraMode->SetDetChamb0(i);
2560 if (fDebugLevel > 1) {
2561 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2563 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2566 //_______________________________________________________________________________
2567 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2570 // Calculate the range of fits
2575 if (fDebugLevel == 1) {
2579 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2581 fDect2 = fNumberOfBinsExpected;
2583 if (fDebugLevel >= 3) {
2584 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2585 fCalibraMode->CalculXBins(fCountDet,i);
2586 fDect1 = fCalibraMode->GetXbins(i);
2587 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2588 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2589 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2590 ,(Int_t) GetStack(fCountDet)
2591 ,(Int_t) GetSector(fCountDet),i);
2592 // Set for the next detector
2593 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2596 //_______________________________________________________________________________
2597 Bool_t AliTRDCalibraFit::CheckFitVoir()
2600 // Check if fFitVoir is in the range
2603 if (fFitVoir < fNumberOfBinsExpected) {
2604 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2607 AliInfo("fFitVoir is out of range of the histo!");
2612 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2613 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2616 // See if we are in a new detector and update the
2617 // variables fNfragZ and fNfragRphi if yes
2618 // Will never happen for only one detector (3 and 4)
2619 // Doesn't matter for 2
2621 if (fCount == idect) {
2622 // On en est au detector (or first detector in the group)
2624 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2625 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2626 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2627 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2628 ,(Int_t) GetStack(fCountDet)
2629 ,(Int_t) GetSector(fCountDet),i);
2630 // Set for the next detector
2631 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2636 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2637 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2640 // Reconstruct the min pad row, max pad row, min pad col and
2641 // max pad col of the calibration group for the Fit functions
2642 // idect is the calibration group inside the detector
2644 if (fDebugLevel != 1) {
2645 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2647 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2648 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2650 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2651 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2654 // For the case where there are not enough entries in the histograms
2655 // of the calibration group, the value present in the choosen database
2656 // will be put. A negativ sign enables to know that a fit was not possible.
2659 if (fDebugLevel == 1) {
2660 AliInfo("The element has not enough statistic to be fitted");
2662 else if (fNbDet > 0){
2663 Int_t firstdetector = fCountDet;
2664 Int_t lastdetector = fCountDet+fNbDet;
2665 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2666 ,idect,firstdetector,lastdetector));
2667 // loop over detectors
2668 for(Int_t det = firstdetector; det < lastdetector; det++){
2670 //Set the calibration object again
2674 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2676 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2677 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2678 ,(Int_t) GetStack(fCountDet)
2679 ,(Int_t) GetSector(fCountDet),0);
2680 // Reconstruct row min row max
2681 ReconstructFitRowMinRowMax(idect,0);
2683 // Calcul the coef from the database choosen for the detector
2684 CalculChargeCoefMean(kFALSE);
2686 //stack 2, not stack 2
2688 if(GetStack(fCountDet) == 2) factor = 12;
2691 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2692 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2693 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2694 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2698 //Put default value negative
2699 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2700 fCurrentCoefE = 0.0;
2705 if(fDebugLevel > 1){
2707 if ( !fDebugStreamer ) {
2709 TDirectory *backup = gDirectory;
2710 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2711 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2714 Int_t detector = fCountDet;
2715 Int_t caligroup = idect;
2716 Short_t rowmin = fCalibraMode->GetRowMin(0);
2717 Short_t rowmax = fCalibraMode->GetRowMax(0);
2718 Short_t colmin = fCalibraMode->GetColMin(0);
2719 Short_t colmax = fCalibraMode->GetColMax(0);
2720 Float_t gf = fCurrentCoef[0];
2721 Float_t gfs = fCurrentCoef[1];
2722 Float_t gfE = fCurrentCoefE;
2724 (*fDebugStreamer) << "FillFillCH" <<
2725 "detector=" << detector <<
2726 "caligroup=" << caligroup <<
2727 "rowmin=" << rowmin <<
2728 "rowmax=" << rowmax <<
2729 "colmin=" << colmin <<
2730 "colmax=" << colmax <<
2738 for (Int_t k = 0; k < 2304; k++) {
2739 fCurrentCoefDetector[k] = 0.0;
2743 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2747 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2748 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2750 // Calcul the coef from the database choosen
2751 CalculChargeCoefMean(kFALSE);
2753 //stack 2, not stack 2
2755 if(GetStack(fCountDet) == 2) factor = 12;
2758 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2759 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2760 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2761 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2765 //Put default value negative
2766 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2767 fCurrentCoefE = 0.0;
2776 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2777 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2780 // For the case where there are not enough entries in the histograms
2781 // of the calibration group, the value present in the choosen database
2782 // will be put. A negativ sign enables to know that a fit was not possible.
2784 if (fDebugLevel == 1) {
2785 AliInfo("The element has not enough statistic to be fitted");
2787 else if (fNbDet > 0) {
2789 Int_t firstdetector = fCountDet;
2790 Int_t lastdetector = fCountDet+fNbDet;
2791 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2792 ,idect,firstdetector,lastdetector));
2793 // loop over detectors
2794 for(Int_t det = firstdetector; det < lastdetector; det++){
2796 //Set the calibration object again
2800 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2802 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
2803 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2804 ,(Int_t) GetStack(fCountDet)
2805 ,(Int_t) GetSector(fCountDet),1);
2806 // Reconstruct row min row max
2807 ReconstructFitRowMinRowMax(idect,1);
2809 // Calcul the coef from the database choosen for the detector
2810 CalculVdriftCoefMean();
2813 //stack 2, not stack 2
2815 if(GetStack(fCountDet) == 2) factor = 12;
2818 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2819 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2820 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2821 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2822 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2826 //Put default value negative
2827 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2828 fCurrentCoefE = 0.0;
2829 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2830 fCurrentCoefE2 = 0.0;
2836 if(fDebugLevel > 1){
2838 if ( !fDebugStreamer ) {
2840 TDirectory *backup = gDirectory;
2841 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
2842 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2846 Int_t detector = fCountDet;
2847 Int_t caligroup = idect;
2848 Short_t rowmin = fCalibraMode->GetRowMin(1);
2849 Short_t rowmax = fCalibraMode->GetRowMax(1);
2850 Short_t colmin = fCalibraMode->GetColMin(1);
2851 Short_t colmax = fCalibraMode->GetColMax(1);
2852 Float_t vf = fCurrentCoef[0];
2853 Float_t vs = fCurrentCoef[1];
2854 Float_t vfE = fCurrentCoefE;
2855 Float_t t0f = fCurrentCoef2[0];
2856 Float_t t0s = fCurrentCoef2[1];
2857 Float_t t0E = fCurrentCoefE2;
2861 (* fDebugStreamer) << "FillFillPH"<<
2862 "detector="<<detector<<
2863 "nentries="<<nentries<<
2864 "caligroup="<<caligroup<<
2878 for (Int_t k = 0; k < 2304; k++) {
2879 fCurrentCoefDetector[k] = 0.0;
2880 fCurrentCoefDetector2[k] = 0.0;
2884 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2888 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2889 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2891 CalculVdriftCoefMean();
2894 //stack 2 and not stack 2
2896 if(GetStack(fCountDet) == 2) factor = 12;
2900 // Fill the fCurrentCoefDetector 2
2901 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2902 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2903 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2904 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2908 // Put the default value
2909 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2910 fCurrentCoefE = 0.0;
2911 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2912 fCurrentCoefE2 = 0.0;
2914 FillFillPH(idect,nentries);
2923 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2924 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2927 // For the case where there are not enough entries in the histograms
2928 // of the calibration group, the value present in the choosen database
2929 // will be put. A negativ sign enables to know that a fit was not possible.
2932 if (fDebugLevel == 1) {
2933 AliInfo("The element has not enough statistic to be fitted");
2935 else if (fNbDet > 0){
2937 Int_t firstdetector = fCountDet;
2938 Int_t lastdetector = fCountDet+fNbDet;
2939 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2940 ,idect,firstdetector,lastdetector));
2942 // loop over detectors
2943 for(Int_t det = firstdetector; det < lastdetector; det++){
2945 //Set the calibration object again
2949 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2951 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
2952 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2953 ,(Int_t) GetStack(fCountDet)
2954 ,(Int_t) GetSector(fCountDet),2);
2955 // Reconstruct row min row max
2956 ReconstructFitRowMinRowMax(idect,2);
2958 // Calcul the coef from the database choosen for the detector
2959 CalculPRFCoefMean();
2961 //stack 2, not stack 2
2963 if(GetStack(fCountDet) == 2) factor = 12;
2966 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2967 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2968 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2969 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2973 //Put default value negative
2974 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2975 fCurrentCoefE = 0.0;
2980 if(fDebugLevel > 1){
2982 if ( !fDebugStreamer ) {
2984 TDirectory *backup = gDirectory;
2985 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
2986 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2989 Int_t detector = fCountDet;
2990 Int_t layer = GetLayer(fCountDet);
2991 Int_t caligroup = idect;
2992 Short_t rowmin = fCalibraMode->GetRowMin(2);
2993 Short_t rowmax = fCalibraMode->GetRowMax(2);
2994 Short_t colmin = fCalibraMode->GetColMin(2);
2995 Short_t colmax = fCalibraMode->GetColMax(2);
2996 Float_t widf = fCurrentCoef[0];
2997 Float_t wids = fCurrentCoef[1];
2998 Float_t widfE = fCurrentCoefE;
3000 (* fDebugStreamer) << "FillFillPRF"<<
3001 "detector="<<detector<<
3003 "caligroup="<<caligroup<<
3014 for (Int_t k = 0; k < 2304; k++) {
3015 fCurrentCoefDetector[k] = 0.0;
3019 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3023 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3024 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3026 CalculPRFCoefMean();
3028 // stack 2 and not stack 2
3030 if(GetStack(fCountDet) == 2) factor = 12;
3034 // Fill the fCurrentCoefDetector
3035 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3036 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3037 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3041 // Put the default value
3042 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3043 fCurrentCoefE = 0.0;
3051 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3052 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3055 // For the case where there are not enough entries in the histograms
3056 // of the calibration group, the value present in the choosen database
3057 // will be put. A negativ sign enables to know that a fit was not possible.
3060 // Calcul the coef from the database choosen
3061 CalculVdriftLorentzCoef();
3064 if(GetStack(fCountDet) == 2) factor = 1728;
3068 // Fill the fCurrentCoefDetector
3069 for (Int_t k = 0; k < factor; k++) {
3070 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3071 // should be negative
3072 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3076 //Put default opposite sign
3077 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3078 fCurrentCoefE = 0.0;
3079 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3080 fCurrentCoefE2 = 0.0;
3082 FillFillLinearFitter();
3087 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3088 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3091 // Fill the coefficients found with the fits or other
3092 // methods from the Fit functions
3095 if (fDebugLevel != 1) {
3097 Int_t firstdetector = fCountDet;
3098 Int_t lastdetector = fCountDet+fNbDet;
3099 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3100 ,idect,firstdetector,lastdetector));
3101 // loop over detectors
3102 for(Int_t det = firstdetector; det < lastdetector; det++){
3104 //Set the calibration object again
3108 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3110 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3111 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3112 ,(Int_t) GetStack(fCountDet)
3113 ,(Int_t) GetSector(fCountDet),0);
3114 // Reconstruct row min row max
3115 ReconstructFitRowMinRowMax(idect,0);
3117 // Calcul the coef from the database choosen for the detector
3118 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3119 else CalculChargeCoefMean(kTRUE);
3121 //stack 2, not stack 2
3123 if(GetStack(fCountDet) == 2) factor = 12;
3126 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3127 Double_t coeftoput = 1.0;
3128 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3129 else coeftoput = fCurrentCoef[0];
3130 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3131 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3132 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3139 if(fDebugLevel > 1){
3141 if ( !fDebugStreamer ) {
3143 TDirectory *backup = gDirectory;
3144 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3145 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3148 Int_t detector = fCountDet;
3149 Int_t caligroup = idect;
3150 Short_t rowmin = fCalibraMode->GetRowMin(0);
3151 Short_t rowmax = fCalibraMode->GetRowMax(0);
3152 Short_t colmin = fCalibraMode->GetColMin(0);
3153 Short_t colmax = fCalibraMode->GetColMax(0);
3154 Float_t gf = fCurrentCoef[0];
3155 Float_t gfs = fCurrentCoef[1];
3156 Float_t gfE = fCurrentCoefE;
3158 (*fDebugStreamer) << "FillFillCH" <<
3159 "detector=" << detector <<
3160 "caligroup=" << caligroup <<
3161 "rowmin=" << rowmin <<
3162 "rowmax=" << rowmax <<
3163 "colmin=" << colmin <<
3164 "colmax=" << colmax <<
3172 for (Int_t k = 0; k < 2304; k++) {
3173 fCurrentCoefDetector[k] = 0.0;
3177 //printf("Check the count now: fCountDet %d\n",fCountDet);
3182 if(GetStack(fCountDet) == 2) factor = 12;
3185 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3186 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3187 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3198 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3199 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3202 // Fill the coefficients found with the fits or other
3203 // methods from the Fit functions
3206 if (fDebugLevel != 1) {
3209 Int_t firstdetector = fCountDet;
3210 Int_t lastdetector = fCountDet+fNbDet;
3211 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3212 ,idect,firstdetector,lastdetector));
3214 // loop over detectors
3215 for(Int_t det = firstdetector; det < lastdetector; det++){
3217 //Set the calibration object again
3221 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3223 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3224 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3225 ,(Int_t) GetStack(fCountDet)
3226 ,(Int_t) GetSector(fCountDet),1);
3227 // Reconstruct row min row max
3228 ReconstructFitRowMinRowMax(idect,1);
3230 // Calcul the coef from the database choosen for the detector
3231 CalculVdriftCoefMean();
3234 //stack 2, not stack 2
3236 if(GetStack(fCountDet) == 2) factor = 12;
3239 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3240 Double_t coeftoput = 1.5;
3241 Double_t coeftoput2 = 0.0;
3243 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3244 else coeftoput = fCurrentCoef[0];
3246 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3247 else coeftoput2 = fCurrentCoef2[0];
3249 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3250 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3251 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3252 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3260 if(fDebugLevel > 1){
3262 if ( !fDebugStreamer ) {
3264 TDirectory *backup = gDirectory;
3265 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3266 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3270 Int_t detector = fCountDet;
3271 Int_t caligroup = idect;
3272 Short_t rowmin = fCalibraMode->GetRowMin(1);
3273 Short_t rowmax = fCalibraMode->GetRowMax(1);
3274 Short_t colmin = fCalibraMode->GetColMin(1);
3275 Short_t colmax = fCalibraMode->GetColMax(1);
3276 Float_t vf = fCurrentCoef[0];
3277 Float_t vs = fCurrentCoef[1];
3278 Float_t vfE = fCurrentCoefE;
3279 Float_t t0f = fCurrentCoef2[0];
3280 Float_t t0s = fCurrentCoef2[1];
3281 Float_t t0E = fCurrentCoefE2;
3285 (* fDebugStreamer) << "FillFillPH"<<
3286 "detector="<<detector<<
3287 "nentries="<<nentries<<
3288 "caligroup="<<caligroup<<
3302 for (Int_t k = 0; k < 2304; k++) {
3303 fCurrentCoefDetector[k] = 0.0;
3304 fCurrentCoefDetector2[k] = 0.0;
3308 //printf("Check the count now: fCountDet %d\n",fCountDet);
3313 if(GetStack(fCountDet) == 2) factor = 12;
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)] = fCurrentCoef[0];
3319 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3323 FillFillPH(idect,nentries);
3328 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3329 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3332 // Fill the coefficients found with the fits or other
3333 // methods from the Fit functions
3336 if (fDebugLevel != 1) {
3339 Int_t firstdetector = fCountDet;
3340 Int_t lastdetector = fCountDet+fNbDet;
3341 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3342 ,idect,firstdetector,lastdetector));
3344 // loop over detectors
3345 for(Int_t det = firstdetector; det < lastdetector; det++){
3347 //Set the calibration object again
3351 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3353 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3354 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3355 ,(Int_t) GetStack(fCountDet)
3356 ,(Int_t) GetSector(fCountDet),2);
3357 // Reconstruct row min row max
3358 ReconstructFitRowMinRowMax(idect,2);
3360 // Calcul the coef from the database choosen for the detector
3361 CalculPRFCoefMean();
3363 //stack 2, not stack 2
3365 if(GetStack(fCountDet) == 2) factor = 12;
3368 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3369 Double_t coeftoput = 1.0;
3370 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3371 else coeftoput = fCurrentCoef[0];
3372 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3373 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3374 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3381 if(fDebugLevel > 1){
3383 if ( !fDebugStreamer ) {
3385 TDirectory *backup = gDirectory;
3386 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3387 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3390 Int_t detector = fCountDet;
3391 Int_t layer = GetLayer(fCountDet);
3392 Int_t caligroup = idect;
3393 Short_t rowmin = fCalibraMode->GetRowMin(2);
3394 Short_t rowmax = fCalibraMode->GetRowMax(2);
3395 Short_t colmin = fCalibraMode->GetColMin(2);
3396 Short_t colmax = fCalibraMode->GetColMax(2);
3397 Float_t widf = fCurrentCoef[0];
3398 Float_t wids = fCurrentCoef[1];
3399 Float_t widfE = fCurrentCoefE;
3401 (* fDebugStreamer) << "FillFillPRF"<<
3402 "detector="<<detector<<
3404 "caligroup="<<caligroup<<
3415 for (Int_t k = 0; k < 2304; k++) {
3416 fCurrentCoefDetector[k] = 0.0;
3420 //printf("Check the count now: fCountDet %d\n",fCountDet);
3425 if(GetStack(fCountDet) == 2) factor = 12;
3428 // Pointer to the branch
3429 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3430 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3431 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3441 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3442 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3445 // Fill the coefficients found with the fits or other
3446 // methods from the Fit functions
3450 if(GetStack(fCountDet) == 2) factor = 1728;
3453 // Pointer to the branch
3454 for (Int_t k = 0; k < factor; k++) {
3455 fCurrentCoefDetector[k] = fCurrentCoef[0];
3456 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3459 FillFillLinearFitter();
3464 //________________________________________________________________________________
3465 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3468 // DebugStream and fVectorFit
3471 // End of one detector
3472 if ((idect == (fCount-1))) {
3475 for (Int_t k = 0; k < 2304; k++) {
3476 fCurrentCoefDetector[k] = 0.0;
3480 if(fDebugLevel > 1){
3482 if ( !fDebugStreamer ) {
3484 TDirectory *backup = gDirectory;
3485 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3486 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3489 Int_t detector = fCountDet;
3490 Int_t caligroup = idect;
3491 Short_t rowmin = fCalibraMode->GetRowMin(0);
3492 Short_t rowmax = fCalibraMode->GetRowMax(0);
3493 Short_t colmin = fCalibraMode->GetColMin(0);
3494 Short_t colmax = fCalibraMode->GetColMax(0);
3495 Float_t gf = fCurrentCoef[0];
3496 Float_t gfs = fCurrentCoef[1];
3497 Float_t gfE = fCurrentCoefE;
3499 (*fDebugStreamer) << "FillFillCH" <<
3500 "detector=" << detector <<
3501 "caligroup=" << caligroup <<
3502 "rowmin=" << rowmin <<
3503 "rowmax=" << rowmax <<
3504 "colmin=" << colmin <<
3505 "colmax=" << colmax <<
3513 //________________________________________________________________________________
3514 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3517 // DebugStream and fVectorFit and fVectorFit2
3520 // End of one detector
3521 if ((idect == (fCount-1))) {
3525 for (Int_t k = 0; k < 2304; k++) {
3526 fCurrentCoefDetector[k] = 0.0;
3527 fCurrentCoefDetector2[k] = 0.0;
3531 if(fDebugLevel > 1){
3533 if ( !fDebugStreamer ) {
3535 TDirectory *backup = gDirectory;
3536 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3537 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3541 Int_t detector = fCountDet;
3542 Int_t caligroup = idect;
3543 Short_t rowmin = fCalibraMode->GetRowMin(1);
3544 Short_t rowmax = fCalibraMode->GetRowMax(1);
3545 Short_t colmin = fCalibraMode->GetColMin(1);
3546 Short_t colmax = fCalibraMode->GetColMax(1);
3547 Float_t vf = fCurrentCoef[0];
3548 Float_t vs = fCurrentCoef[1];
3549 Float_t vfE = fCurrentCoefE;
3550 Float_t t0f = fCurrentCoef2[0];
3551 Float_t t0s = fCurrentCoef2[1];
3552 Float_t t0E = fCurrentCoefE2;
3556 (* fDebugStreamer) << "FillFillPH"<<
3557 "detector="<<detector<<
3558 "nentries="<<nentries<<
3559 "caligroup="<<caligroup<<
3574 //________________________________________________________________________________
3575 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3578 // DebugStream and fVectorFit
3581 // End of one detector
3582 if ((idect == (fCount-1))) {
3585 for (Int_t k = 0; k < 2304; k++) {
3586 fCurrentCoefDetector[k] = 0.0;
3591 if(fDebugLevel > 1){
3593 if ( !fDebugStreamer ) {
3595 TDirectory *backup = gDirectory;
3596 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3597 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3600 Int_t detector = fCountDet;
3601 Int_t layer = GetLayer(fCountDet);
3602 Int_t caligroup = idect;
3603 Short_t rowmin = fCalibraMode->GetRowMin(2);
3604 Short_t rowmax = fCalibraMode->GetRowMax(2);
3605 Short_t colmin = fCalibraMode->GetColMin(2);
3606 Short_t colmax = fCalibraMode->GetColMax(2);
3607 Float_t widf = fCurrentCoef[0];
3608 Float_t wids = fCurrentCoef[1];
3609 Float_t widfE = fCurrentCoefE;
3611 (* fDebugStreamer) << "FillFillPRF"<<
3612 "detector="<<detector<<
3614 "caligroup="<<caligroup<<
3626 //________________________________________________________________________________
3627 void AliTRDCalibraFit::FillFillLinearFitter()
3630 // DebugStream and fVectorFit
3633 // End of one detector
3639 for (Int_t k = 0; k < 2304; k++) {
3640 fCurrentCoefDetector[k] = 0.0;
3641 fCurrentCoefDetector2[k] = 0.0;
3645 if(fDebugLevel > 1){
3647 if ( !fDebugStreamer ) {
3649 TDirectory *backup = gDirectory;
3650 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3651 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3654 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3655 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3656 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3657 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3658 Float_t tiltangle = padplane->GetTiltingAngle();
3659 Int_t detector = fCountDet;
3660 Int_t stack = GetStack(fCountDet);
3661 Int_t layer = GetLayer(fCountDet);
3662 Float_t vf = fCurrentCoef[0];
3663 Float_t vs = fCurrentCoef[1];
3664 Float_t vfE = fCurrentCoefE;
3665 Float_t lorentzangler = fCurrentCoef2[0];
3666 Float_t elorentzangler = fCurrentCoefE2;
3667 Float_t lorentzangles = fCurrentCoef2[1];
3669 (* fDebugStreamer) << "FillFillLinearFitter"<<
3670 "detector="<<detector<<
3675 "tiltangle="<<tiltangle<<
3679 "lorentzangler="<<lorentzangler<<
3680 "Elorentzangler="<<elorentzangler<<
3681 "lorentzangles="<<lorentzangles<<
3687 //____________Calcul Coef Mean_________________________________________________
3689 //_____________________________________________________________________________
3690 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3693 // For the detector Dect calcul the mean time 0
3694 // for the calibration group idect from the choosen database
3697 fCurrentCoef2[1] = 0.0;
3698 if(fDebugLevel != 1){
3699 if(((fCalibraMode->GetNz(1) > 0) ||
3700 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3702 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3703 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3704 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3708 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3714 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3718 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3719 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3720 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3723 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3731 //_____________________________________________________________________________
3732 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3735 // For the detector Dect calcul the mean gain factor
3736 // for the calibration group idect from the choosen database
3739 fCurrentCoef[1] = 0.0;
3740 if(fDebugLevel != 1){
3741 if (((fCalibraMode->GetNz(0) > 0) ||
3742 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3743 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3744 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3745 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3746 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3749 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3753 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3754 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3759 //_____________________________________________________________________________
3760 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3763 // For the detector Dect calcul the mean sigma of pad response
3764 // function for the calibration group idect from the choosen database
3767 fCurrentCoef[1] = 0.0;
3768 if(fDebugLevel != 1){
3769 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3770 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3771 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3774 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3778 //_____________________________________________________________________________
3779 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3782 // For the detector dect calcul the mean drift velocity for the
3783 // calibration group idect from the choosen database
3786 fCurrentCoef[1] = 0.0;
3787 if(fDebugLevel != 1){
3788 if (((fCalibraMode->GetNz(1) > 0) ||
3789 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3791 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3792 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3793 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3797 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3802 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3807 //_____________________________________________________________________________
3808 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
3811 // For the detector fCountDet, mean drift velocity and tan lorentzangle
3814 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
3815 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3819 //_____________________________________________________________________________
3820 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
3823 // Default width of the PRF if there is no database as reference
3828 //case 0: return 0.515;
3829 //case 1: return 0.502;
3830 //case 2: return 0.491;
3831 //case 3: return 0.481;
3832 //case 4: return 0.471;
3833 //case 5: return 0.463;
3834 //default: return 0.0;
3837 case 0: return 0.538429;
3838 case 1: return 0.524302;
3839 case 2: return 0.511591;
3840 case 3: return 0.500140;
3841 case 4: return 0.489821;
3842 case 5: return 0.480524;
3843 default: return 0.0;
3846 //________________________________________________________________________________
3847 void AliTRDCalibraFit::SetCalROC(Int_t i)
3850 // Set the calib object for fCountDet
3853 Float_t value = 0.0;
3855 //Get the CalDet object
3857 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3859 AliInfo("Could not get calibDB");
3866 fCalROC->~AliTRDCalROC();
3867 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3868 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3872 fCalROC->~AliTRDCalROC();
3873 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3874 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3876 fCalROC2->~AliTRDCalROC();
3877 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3878 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3882 fCalROC->~AliTRDCalROC();
3883 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3884 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3893 if(fCalROC) delete fCalROC;
3894 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3895 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3896 fCalROC->SetValue(k,1.0);
3900 if(fCalROC) delete fCalROC;
3901 if(fCalROC2) delete fCalROC2;
3902 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3903 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3904 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3905 fCalROC->SetValue(k,1.0);
3906 fCalROC2->SetValue(k,0.0);
3910 if(fCalROC) delete fCalROC;
3911 value = GetPRFDefault(GetLayer(fCountDet));
3912 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3913 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3914 fCalROC->SetValue(k,value);
3922 //____________Fit Methods______________________________________________________
3924 //_____________________________________________________________________________
3925 void AliTRDCalibraFit::FitPente(TH1* projPH)
3928 // Slope methode for the drift velocity
3932 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3939 fCurrentCoefE = 0.0;
3940 fCurrentCoefE2 = 0.0;
3941 fCurrentCoef[0] = 0.0;
3942 fCurrentCoef2[0] = 0.0;
3943 TLine *line = new TLine();
3946 TAxis *xpph = projPH->GetXaxis();
3947 Int_t nbins = xpph->GetNbins();
3948 Double_t lowedge = xpph->GetBinLowEdge(1);
3949 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3950 Double_t widbins = (upedge - lowedge) / nbins;
3951 Double_t limit = upedge + 0.5 * widbins;
3954 // Beginning of the signal
3955 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3956 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3957 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3959 binmax = (Int_t) pentea->GetMaximumBin();
3962 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3964 if (binmax >= nbins) {
3967 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3969 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3970 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3971 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3972 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3973 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3975 fPhd[0] = -(l3P1am / (2 * l3P2am));
3978 if((l3P1am != 0.0) && (l3P2am != 0.0)){
3979 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3982 // Amplification region
3985 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3986 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3993 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3995 if (binmax >= nbins) {
3998 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4000 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4001 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4002 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4003 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4004 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4006 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4008 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
4009 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4012 fCurrentCoefE2 = fCurrentCoefE;
4015 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4016 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4017 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4020 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4023 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4025 if (binmin >= nbins) {
4028 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4033 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4034 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4035 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4036 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4037 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4038 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4040 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4042 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
4043 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4045 Float_t fPhdt0 = 0.0;
4046 Float_t t0Shift = 0.0;
4049 t0Shift = fT0Shift1;
4053 t0Shift = fT0Shift0;
4056 if ((fPhd[2] > fPhd[0]) &&
4057 (fPhd[2] > fPhd[1]) &&
4058 (fPhd[1] > fPhd[0]) &&
4060 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4061 fNumberFitSuccess++;
4063 if (fPhdt0 >= 0.0) {
4064 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4065 if (fCurrentCoef2[0] < -1.0) {
4066 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4070 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4075 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4076 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4079 if (fDebugLevel == 1) {
4080 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4083 line->SetLineColor(2);
4084 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4085 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4086 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4087 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4088 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4089 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4090 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4091 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4094 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4103 //_____________________________________________________________________________
4104 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4107 // Slope methode but with polynomes de Lagrange
4111 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4114 Double_t *x = new Double_t[5];
4115 Double_t *y = new Double_t[5];
4130 fCurrentCoefE = 0.0;
4131 fCurrentCoefE2 = 1.0;
4132 fCurrentCoef[0] = 0.0;
4133 fCurrentCoef2[0] = 0.0;
4134 TLine *line = new TLine();
4135 TF1 * polynome = 0x0;
4136 TF1 * polynomea = 0x0;
4137 TF1 * polynomeb = 0x0;
4141 TAxis *xpph = projPH->GetXaxis();
4142 Int_t nbins = xpph->GetNbins();
4143 Double_t lowedge = xpph->GetBinLowEdge(1);
4144 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4145 Double_t widbins = (upedge - lowedge) / nbins;
4146 Double_t limit = upedge + 0.5 * widbins;
4151 // Beginning of the signal
4152 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4153 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4154 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4157 binmax = (Int_t) pentea->GetMaximumBin();
4159 Double_t minnn = 0.0;
4160 Double_t maxxx = 0.0;
4162 Int_t kase = nbins-binmax;
4170 minnn = pentea->GetBinCenter(binmax-2);
4171 maxxx = pentea->GetBinCenter(binmax);
4172 x[0] = pentea->GetBinCenter(binmax-2);
4173 x[1] = pentea->GetBinCenter(binmax-1);
4174 x[2] = pentea->GetBinCenter(binmax);
4175 y[0] = pentea->GetBinContent(binmax-2);
4176 y[1] = pentea->GetBinContent(binmax-1);
4177 y[2] = pentea->GetBinContent(binmax);
4178 c = CalculPolynomeLagrange2(x,y);
4179 AliInfo("At the limit for beginning!");
4182 minnn = pentea->GetBinCenter(binmax-2);
4183 maxxx = pentea->GetBinCenter(binmax+1);
4184 x[0] = pentea->GetBinCenter(binmax-2);
4185 x[1] = pentea->GetBinCenter(binmax-1);
4186 x[2] = pentea->GetBinCenter(binmax);
4187 x[3] = pentea->GetBinCenter(binmax+1);
4188 y[0] = pentea->GetBinContent(binmax-2);
4189 y[1] = pentea->GetBinContent(binmax-1);
4190 y[2] = pentea->GetBinContent(binmax);
4191 y[3] = pentea->GetBinContent(binmax+1);
4192 c = CalculPolynomeLagrange3(x,y);
4200 minnn = pentea->GetBinCenter(binmax);
4201 maxxx = pentea->GetBinCenter(binmax+2);
4202 x[0] = pentea->GetBinCenter(binmax);
4203 x[1] = pentea->GetBinCenter(binmax+1);
4204 x[2] = pentea->GetBinCenter(binmax+2);
4205 y[0] = pentea->GetBinContent(binmax);
4206 y[1] = pentea->GetBinContent(binmax+1);
4207 y[2] = pentea->GetBinContent(binmax+2);
4208 c = CalculPolynomeLagrange2(x,y);
4211 minnn = pentea->GetBinCenter(binmax-1);
4212 maxxx = pentea->GetBinCenter(binmax+2);
4213 x[0] = pentea->GetBinCenter(binmax-1);
4214 x[1] = pentea->GetBinCenter(binmax);
4215 x[2] = pentea->GetBinCenter(binmax+1);
4216 x[3] = pentea->GetBinCenter(binmax+2);
4217 y[0] = pentea->GetBinContent(binmax-1);
4218 y[1] = pentea->GetBinContent(binmax);
4219 y[2] = pentea->GetBinContent(binmax+1);
4220 y[3] = pentea->GetBinContent(binmax+2);
4221 c = CalculPolynomeLagrange3(x,y);
4224 minnn = pentea->GetBinCenter(binmax-2);
4225 maxxx = pentea->GetBinCenter(binmax+2);
4226 x[0] = pentea->GetBinCenter(binmax-2);
4227 x[1] = pentea->GetBinCenter(binmax-1);
4228 x[2] = pentea->GetBinCenter(binmax);
4229 x[3] = pentea->GetBinCenter(binmax+1);
4230 x[4] = pentea->GetBinCenter(binmax+2);
4231 y[0] = pentea->GetBinContent(binmax-2);
4232 y[1] = pentea->GetBinContent(binmax-1);
4233 y[2] = pentea->GetBinContent(binmax);
4234 y[3] = pentea->GetBinContent(binmax+1);
4235 y[4] = pentea->GetBinContent(binmax+2);
4236 c = CalculPolynomeLagrange4(x,y);
4244 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4245 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4247 Double_t step = (maxxx-minnn)/10000;
4249 Double_t maxvalue = 0.0;
4250 Double_t placemaximum = minnn;
4251 for(Int_t o = 0; o < 10000; o++){
4252 if(o == 0) maxvalue = polynomeb->Eval(l);
4253 if(maxvalue < (polynomeb->Eval(l))){
4254 maxvalue = polynomeb->Eval(l);
4259 fPhd[0] = placemaximum;
4262 // Amplification region
4265 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4266 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4272 Double_t minn = 0.0;
4273 Double_t maxx = 0.0;
4285 Int_t kase1 = nbins - binmax;
4287 //Determination of minn and maxx
4288 //case binmax = nbins
4293 minn = projPH->GetBinCenter(binmax-2);
4294 maxx = projPH->GetBinCenter(binmax);
4295 x[0] = projPH->GetBinCenter(binmax-2);
4296 x[1] = projPH->GetBinCenter(binmax-1);
4297 x[2] = projPH->GetBinCenter(binmax);
4298 y[0] = projPH->GetBinContent(binmax-2);
4299 y[1] = projPH->GetBinContent(binmax-1);
4300 y[2] = projPH->GetBinContent(binmax);
4301 c = CalculPolynomeLagrange2(x,y);
4302 //AliInfo("At the limit for the drift!");
4305 minn = projPH->GetBinCenter(binmax-2);
4306 maxx = projPH->GetBinCenter(binmax+1);
4307 x[0] = projPH->GetBinCenter(binmax-2);
4308 x[1] = projPH->GetBinCenter(binmax-1);
4309 x[2] = projPH->GetBinCenter(binmax);
4310 x[3] = projPH->GetBinCenter(binmax+1);
4311 y[0] = projPH->GetBinContent(binmax-2);
4312 y[1] = projPH->GetBinContent(binmax-1);
4313 y[2] = projPH->GetBinContent(binmax);
4314 y[3] = projPH->GetBinContent(binmax+1);
4315 c = CalculPolynomeLagrange3(x,y);
4324 minn = projPH->GetBinCenter(binmax);
4325 maxx = projPH->GetBinCenter(binmax+2);
4326 x[0] = projPH->GetBinCenter(binmax);
4327 x[1] = projPH->GetBinCenter(binmax+1);
4328 x[2] = projPH->GetBinCenter(binmax+2);
4329 y[0] = projPH->GetBinContent(binmax);
4330 y[1] = projPH->GetBinContent(binmax+1);
4331 y[2] = projPH->GetBinContent(binmax+2);
4332 c = CalculPolynomeLagrange2(x,y);
4335 minn = projPH->GetBinCenter(binmax-1);
4336 maxx = projPH->GetBinCenter(binmax+2);
4337 x[0] = projPH->GetBinCenter(binmax-1);
4338 x[1] = projPH->GetBinCenter(binmax);
4339 x[2] = projPH->GetBinCenter(binmax+1);
4340 x[3] = projPH->GetBinCenter(binmax+2);
4341 y[0] = projPH->GetBinContent(binmax-1);
4342 y[1] = projPH->GetBinContent(binmax);
4343 y[2] = projPH->GetBinContent(binmax+1);
4344 y[3] = projPH->GetBinContent(binmax+2);
4345 c = CalculPolynomeLagrange3(x,y);
4348 minn = projPH->GetBinCenter(binmax-2);
4349 maxx = projPH->GetBinCenter(binmax+2);
4350 x[0] = projPH->GetBinCenter(binmax-2);
4351 x[1] = projPH->GetBinCenter(binmax-1);
4352 x[2] = projPH->GetBinCenter(binmax);
4353 x[3] = projPH->GetBinCenter(binmax+1);
4354 x[4] = projPH->GetBinCenter(binmax+2);
4355 y[0] = projPH->GetBinContent(binmax-2);
4356 y[1] = projPH->GetBinContent(binmax-1);
4357 y[2] = projPH->GetBinContent(binmax);
4358 y[3] = projPH->GetBinContent(binmax+1);
4359 y[4] = projPH->GetBinContent(binmax+2);
4360 c = CalculPolynomeLagrange4(x,y);
4367 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4368 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4370 Double_t step = (maxx-minn)/1000;
4372 Double_t maxvalue = 0.0;
4373 Double_t placemaximum = minn;
4374 for(Int_t o = 0; o < 1000; o++){
4375 if(o == 0) maxvalue = polynomea->Eval(l);
4376 if(maxvalue < (polynomea->Eval(l))){
4377 maxvalue = polynomea->Eval(l);
4382 fPhd[1] = placemaximum;
4386 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4387 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4388 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4391 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4397 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4401 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4402 AliInfo("Too many fluctuations at the end!");
4405 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4406 AliInfo("Too many fluctuations at the end!");
4409 if(pente->GetBinContent(binmin+1)==0){
4410 AliInfo("No entries for the next bin!");
4411 pente->SetBinContent(binmin,0);
4412 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4428 Bool_t case1 = kFALSE;
4429 Bool_t case2 = kFALSE;
4430 Bool_t case4 = kFALSE;
4432 //Determination of min and max
4433 //case binmin <= nbins-3
4435 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4436 min = pente->GetBinCenter(binmin-2);
4437 max = pente->GetBinCenter(binmin+2);
4438 x[0] = pente->GetBinCenter(binmin-2);
4439 x[1] = pente->GetBinCenter(binmin-1);
4440 x[2] = pente->GetBinCenter(binmin);
4441 x[3] = pente->GetBinCenter(binmin+1);
4442 x[4] = pente->GetBinCenter(binmin+2);
4443 y[0] = pente->GetBinContent(binmin-2);
4444 y[1] = pente->GetBinContent(binmin-1);
4445 y[2] = pente->GetBinContent(binmin);
4446 y[3] = pente->GetBinContent(binmin+1);
4447 y[4] = pente->GetBinContent(binmin+2);
4448 //Calcul the polynome de Lagrange
4449 c = CalculPolynomeLagrange4(x,y);
4451 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4452 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4453 //AliInfo("polynome 4 false 1");
4456 if(((binmin+3) <= (nbins-1)) &&
4457 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4458 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4459 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4460 AliInfo("polynome 4 false 2");
4464 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4465 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4466 //AliInfo("polynome 4 case 1");
4469 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4470 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4471 //AliInfo("polynome 4 case 4");
4476 //case binmin = nbins-2
4478 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4480 min = pente->GetBinCenter(binmin-2);
4481 max = pente->GetBinCenter(binmin+1);
4482 x[0] = pente->GetBinCenter(binmin-2);
4483 x[1] = pente->GetBinCenter(binmin-1);
4484 x[2] = pente->GetBinCenter(binmin);
4485 x[3] = pente->GetBinCenter(binmin+1);
4486 y[0] = pente->GetBinContent(binmin-2);
4487 y[1] = pente->GetBinContent(binmin-1);
4488 y[2] = pente->GetBinContent(binmin);
4489 y[3] = pente->GetBinContent(binmin+1);
4490 //Calcul the polynome de Lagrange
4491 c = CalculPolynomeLagrange3(x,y);
4492 //richtung +: nothing
4494 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4495 //AliInfo("polynome 3- case 2");
4500 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4502 min = pente->GetBinCenter(binmin-1);
4503 max = pente->GetBinCenter(binmin+2);
4504 x[0] = pente->GetBinCenter(binmin-1);
4505 x[1] = pente->GetBinCenter(binmin);
4506 x[2] = pente->GetBinCenter(binmin+1);
4507 x[3] = pente->GetBinCenter(binmin+2);
4508 y[0] = pente->GetBinContent(binmin-1);
4509 y[1] = pente->GetBinContent(binmin);
4510 y[2] = pente->GetBinContent(binmin+1);
4511 y[3] = pente->GetBinContent(binmin+2);
4512 //Calcul the polynome de Lagrange
4513 c = CalculPolynomeLagrange3(x,y);
4515 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4516 //AliInfo("polynome 3+ case 2");
4521 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4522 min = pente->GetBinCenter(binmin);
4523 max = pente->GetBinCenter(binmin+2);
4524 x[0] = pente->GetBinCenter(binmin);
4525 x[1] = pente->GetBinCenter(binmin+1);
4526 x[2] = pente->GetBinCenter(binmin+2);
4527 y[0] = pente->GetBinContent(binmin);
4528 y[1] = pente->GetBinContent(binmin+1);
4529 y[2] = pente->GetBinContent(binmin+2);
4530 //Calcul the polynome de Lagrange
4531 c = CalculPolynomeLagrange2(x,y);
4533 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4534 //AliInfo("polynome 2+ false");
4539 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4541 min = pente->GetBinCenter(binmin-1);
4542 max = pente->GetBinCenter(binmin+1);
4543 x[0] = pente->GetBinCenter(binmin-1);
4544 x[1] = pente->GetBinCenter(binmin);
4545 x[2] = pente->GetBinCenter(binmin+1);
4546 y[0] = pente->GetBinContent(binmin-1);
4547 y[1] = pente->GetBinContent(binmin);
4548 y[2] = pente->GetBinContent(binmin+1);
4549 //Calcul the polynome de Lagrange
4550 c = CalculPolynomeLagrange2(x,y);
4551 //richtung +: nothing
4552 //richtung -: nothing
4554 //case binmin = nbins-1
4556 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4557 min = pente->GetBinCenter(binmin-2);
4558 max = pente->GetBinCenter(binmin);
4559 x[0] = pente->GetBinCenter(binmin-2);
4560 x[1] = pente->GetBinCenter(binmin-1);
4561 x[2] = pente->GetBinCenter(binmin);
4562 y[0] = pente->GetBinContent(binmin-2);
4563 y[1] = pente->GetBinContent(binmin-1);
4564 y[2] = pente->GetBinContent(binmin);
4565 //Calcul the polynome de Lagrange
4566 c = CalculPolynomeLagrange2(x,y);
4567 //AliInfo("At the limit for the drift!");
4568 //fluctuation too big!
4569 //richtung +: nothing
4571 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4572 //AliInfo("polynome 2- false ");
4576 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4578 AliInfo("At the limit for the drift and not usable!");
4582 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4584 AliInfo("For the drift...problem!");
4586 //pass but should not happen
4587 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4589 AliInfo("For the drift...problem!");
4593 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4594 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4595 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4596 Double_t step = (max-min)/1000;
4598 Double_t minvalue = 0.0;
4599 Double_t placeminimum = min;
4600 for(Int_t o = 0; o < 1000; o++){
4601 if(o == 0) minvalue = polynome->Eval(l);
4602 if(minvalue > (polynome->Eval(l))){
4603 minvalue = polynome->Eval(l);
4608 fPhd[2] = placeminimum;
4610 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4611 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4612 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4614 Float_t fPhdt0 = 0.0;
4615 Float_t t0Shift = 0.0;
4618 t0Shift = fT0Shift1;
4622 t0Shift = fT0Shift0;
4625 if ((fPhd[2] > fPhd[0]) &&
4626 (fPhd[2] > fPhd[1]) &&
4627 (fPhd[1] > fPhd[0]) &&
4629 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4630 fNumberFitSuccess++;
4631 if (fPhdt0 >= 0.0) {
4632 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4633 if (fCurrentCoef2[0] < -1.0) {
4634 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4638 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4642 if((fPhd[1] > fPhd[0]) &&
4644 if (fPhdt0 >= 0.0) {
4645 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4646 if (fCurrentCoef2[0] < -1.0) {
4647 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4651 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4655 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4656 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4657 //printf("Fit failed!\n");
4661 if (fDebugLevel == 1) {
4662 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4665 line->SetLineColor(2);
4666 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4667 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4668 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4669 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4670 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4671 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4672 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4673 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4676 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4681 if(pentea) delete pentea;
4682 if(pente) delete pente;
4683 if(polynome) delete polynome;
4684 if(polynomea) delete polynomea;
4685 if(polynomeb) delete polynomeb;
4689 if(line) delete line;
4694 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4695 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4697 projPH->SetDirectory(0);
4701 //_____________________________________________________________________________
4702 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4705 // Fit methode for the drift velocity
4709 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4712 TAxis *xpph = projPH->GetXaxis();
4713 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4715 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4716 fPH->SetParameter(0,0.469); // Scaling
4717 fPH->SetParameter(1,0.18); // Start
4718 fPH->SetParameter(2,0.0857325); // AR
4719 fPH->SetParameter(3,1.89); // DR
4720 fPH->SetParameter(4,0.08); // QA/QD
4721 fPH->SetParameter(5,0.0); // Baseline
4723 TLine *line = new TLine();
4725 fCurrentCoef[0] = 0.0;
4726 fCurrentCoef2[0] = 0.0;
4727 fCurrentCoefE = 0.0;
4728 fCurrentCoefE2 = 0.0;
4730 if (idect%fFitPHPeriode == 0) {
4732 AliInfo(Form("The detector %d will be fitted",idect));
4733 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4734 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4735 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4736 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4737 fPH->SetParameter(4,0.225); // QA/QD
4738 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4740 if (fDebugLevel != 1) {
4741 projPH->Fit(fPH,"0M","",0.0,upedge);
4744 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4746 projPH->Fit(fPH,"M+","",0.0,upedge);
4748 line->SetLineColor(4);
4749 line->DrawLine(fPH->GetParameter(1)
4751 ,fPH->GetParameter(1)
4752 ,projPH->GetMaximum());
4753 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4755 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4756 ,projPH->GetMaximum());
4757 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4759 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4760 ,projPH->GetMaximum());
4763 if (fPH->GetParameter(3) != 0) {
4764 fNumberFitSuccess++;
4765 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4766 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4767 fCurrentCoef2[0] = fPH->GetParameter(1);
4768 fCurrentCoefE2 = fPH->GetParError(1);
4771 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4772 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4778 // Put the default value
4779 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4780 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4783 if (fDebugLevel != 1) {
4788 //_____________________________________________________________________________
4789 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
4792 // Fit methode for the sigma of the pad response function
4797 fCurrentCoef[0] = 0.0;
4798 fCurrentCoefE = 0.0;
4800 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
4803 fCurrentCoef[0] = -fCurrentCoef[1];
4807 fNumberFitSuccess++;
4808 fCurrentCoef[0] = param[2];
4809 fCurrentCoefE = ret;
4813 //_____________________________________________________________________________
4814 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)
4817 // Fit methode for the sigma of the pad response function
4820 //We should have at least 3 points
4821 if(nBins <=3) return -4.0;
4823 TLinearFitter fitter(3,"pol2");
4824 fitter.StoreData(kFALSE);
4825 fitter.ClearPoints();
4827 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
4828 Float_t entries = 0;
4829 Int_t nbbinwithentries = 0;
4830 for (Int_t i=0; i<nBins; i++){
4832 if(arraye[i] > 15) nbbinwithentries++;
4833 //printf("entries for i %d: %f\n",i,arraye[i]);
4835 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
4836 //printf("entries %f\n",entries);
4837 //printf("nbbinwithentries %d\n",nbbinwithentries);
4840 Float_t errorm = 0.0;
4841 Float_t errorn = 0.0;
4842 Float_t error = 0.0;
4845 for (Int_t ibin=0;ibin<nBins; ibin++){
4846 Float_t entriesI = arraye[ibin];
4847 Float_t valueI = arraym[ibin];
4848 Double_t xcenter = 0.0;
4850 if ((entriesI>15) && (valueI>0.0)){
4851 xcenter = xMin+(ibin+0.5)*binWidth;
4856 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
4857 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
4858 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4861 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
4862 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
4863 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4865 if(error == 0.0) continue;
4866 val = TMath::Log(Float_t(valueI));
4867 fitter.AddPoint(&xcenter,val,error);
4871 if(fDebugLevel > 1){
4873 if ( !fDebugStreamer ) {
4875 TDirectory *backup = gDirectory;
4876 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4877 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4880 Int_t detector = fCountDet;
4881 Int_t layer = GetLayer(fCountDet);
4884 (* fDebugStreamer) << "FitGausMIFill"<<
4885 "detector="<<detector<<
4889 "entriesI="<<entriesI<<
4892 "xcenter="<<xcenter<<
4902 if(npoints <=3) return -4.0;
4907 fitter.GetParameters(par);
4908 chi2 = fitter.GetChisquare()/Float_t(npoints);
4911 if (!param) param = new TVectorD(3);
4912 if(par[2] == 0.0) return -4.0;
4913 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
4914 Double_t deltax = (fitter.GetParError(2))/x;
4915 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
4918 (*param)[1] = par[1]/(-2.*par[2]);
4919 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
4920 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
4921 if ( lnparam0>307 ) return -4;
4922 (*param)[0] = TMath::Exp(lnparam0);
4924 if(fDebugLevel > 1){
4926 if ( !fDebugStreamer ) {
4928 TDirectory *backup = gDirectory;
4929 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4930 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4933 Int_t detector = fCountDet;
4934 Int_t layer = GetLayer(fCountDet);
4937 (* fDebugStreamer) << "FitGausMIFit"<<
4938 "detector="<<detector<<
4941 "errorsigma="<<chi2<<
4942 "mean="<<(*param)[1]<<
4943 "sigma="<<(*param)[2]<<
4944 "constant="<<(*param)[0]<<
4949 if((chi2/(*param)[2]) > 0.1){
4951 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
4956 if(fDebugLevel == 1){
4957 TString name("PRF");
4958 name += (Int_t)xMin;
4959 name += (Int_t)xMax;
4960 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
4963 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
4964 for(Int_t k = 0; k < nBins; k++){
4965 histo->SetBinContent(k+1,arraym[k]);
4966 histo->SetBinError(k+1,arrayme[k]);
4969 name += "functionf";
4970 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
4971 f1->SetParameter(0, (*param)[0]);
4972 f1->SetParameter(1, (*param)[1]);
4973 f1->SetParameter(2, (*param)[2]);
4981 //_____________________________________________________________________________
4982 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
4985 // Fit methode for the sigma of the pad response function
4988 fCurrentCoef[0] = 0.0;
4989 fCurrentCoefE = 0.0;
4991 if (fDebugLevel != 1) {
4992 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
4995 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4997 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5000 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5001 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5002 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5004 fNumberFitSuccess++;
5007 //_____________________________________________________________________________
5008 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5011 // Fit methode for the sigma of the pad response function
5013 fCurrentCoef[0] = 0.0;
5014 fCurrentCoefE = 0.0;
5015 if (fDebugLevel == 1) {
5016 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5020 fCurrentCoef[0] = projPRF->GetRMS();
5021 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5023 fNumberFitSuccess++;
5026 //_____________________________________________________________________________
5027 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5030 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5033 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5036 Int_t nbins = (Int_t)(nybins/(2*nbg));
5037 Float_t lowedge = -3.0*nbg;
5038 Float_t upedge = lowedge + 3.0;
5041 Double_t xvalues = -0.2*nbg+0.1;
5043 Int_t total = 2*nbg;
5046 for(Int_t k = 0; k < total; k++){
5047 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5049 y = fCurrentCoef[0]*fCurrentCoef[0];
5050 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5053 if(fDebugLevel > 1){
5055 if ( !fDebugStreamer ) {
5057 TDirectory *backup = gDirectory;
5058 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5059 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5062 Int_t detector = fCountDet;
5063 Int_t layer = GetLayer(fCountDet);
5064 Int_t nbtotal = total;
5066 Float_t low = lowedge;
5067 Float_t up = upedge;
5068 Float_t tnp = xvalues;
5069 Float_t wid = fCurrentCoef[0];
5070 Float_t widfE = fCurrentCoefE;
5072 (* fDebugStreamer) << "FitTnpRange0"<<
5073 "detector="<<detector<<
5075 "nbtotal="<<nbtotal<<
5093 fCurrentCoefE = 0.0;
5094 fCurrentCoef[0] = 0.0;
5096 //printf("npoints\n",npoints);
5099 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5104 linearfitter.Eval();
5105 linearfitter.GetParameters(pars0);
5106 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5107 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5108 Double_t min0 = 0.0;
5109 Double_t ermin0 = 0.0;
5110 //Double_t prfe0 = 0.0;
5111 Double_t prf0 = 0.0;
5112 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
5113 min0 = -pars0[1]/(2*pars0[2]);
5114 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5115 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5118 prfe0 = linearfitter->GetParError(0)*pointError0
5119 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5120 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5121 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5122 fCurrentCoefE = (Float_t) prfe0;
5124 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5127 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5131 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5134 if(fDebugLevel > 1){
5136 if ( !fDebugStreamer ) {
5138 TDirectory *backup = gDirectory;
5139 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5140 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5143 Int_t detector = fCountDet;
5144 Int_t layer = GetLayer(fCountDet);
5145 Int_t nbtotal = total;
5146 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5147 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5149 (* fDebugStreamer) << "FitTnpRange1"<<
5150 "detector="<<detector<<
5152 "nbtotal="<<nbtotal<<
5156 "npoints="<<npoints<<
5159 "sigmaprf="<<fCurrentCoef[0]<<
5160 "sigprf="<<fCurrentCoef[1]<<
5167 //_____________________________________________________________________________
5168 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5171 // Only mean methode for the gain factor
5174 fCurrentCoef[0] = mean;
5175 fCurrentCoefE = 0.0;
5176 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5177 if (fDebugLevel == 1) {
5178 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5182 CalculChargeCoefMean(kTRUE);
5183 fNumberFitSuccess++;
5185 //_____________________________________________________________________________
5186 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5189 // mean w methode for the gain factor
5193 Int_t nybins = projch->GetNbinsX();
5195 //The weight function
5196 Double_t a = 0.00228515;
5197 Double_t b = -0.00231487;
5198 Double_t c = 0.00044298;
5199 Double_t d = -0.00379239;
5200 Double_t e = 0.00338349;
5210 //A arbitrary error for the moment
5211 fCurrentCoefE = 0.0;
5212 fCurrentCoef[0] = 0.0;
5215 Double_t sumw = 0.0;
5217 Float_t sumAll = (Float_t) nentries;
5218 Int_t sumCurrent = 0;
5219 for(Int_t k = 0; k <nybins; k++){
5220 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5221 if (fraction>0.95) break;
5222 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5223 e*fraction*fraction*fraction*fraction;
5224 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5225 sum += weight*projch->GetBinContent(k+1);
5226 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5227 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5229 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5231 if (fDebugLevel == 1) {
5232 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5236 fNumberFitSuccess++;
5237 CalculChargeCoefMean(kTRUE);
5239 //_____________________________________________________________________________
5240 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5243 // mean w methode for the gain factor
5247 Int_t nybins = projch->GetNbinsX();
5249 //The weight function
5250 Double_t a = 0.00228515;
5251 Double_t b = -0.00231487;
5252 Double_t c = 0.00044298;
5253 Double_t d = -0.00379239;
5254 Double_t e = 0.00338349;
5264 //A arbitrary error for the moment
5265 fCurrentCoefE = 0.0;
5266 fCurrentCoef[0] = 0.0;
5269 Double_t sumw = 0.0;
5271 Int_t sumCurrent = 0;
5272 for(Int_t k = 0; k <nybins; k++){
5273 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5274 if (fraction>0.95) break;
5275 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5276 e*fraction*fraction*fraction*fraction;
5277 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5278 sum += weight*projch->GetBinContent(k+1);
5279 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5280 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5282 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5284 if (fDebugLevel == 1) {
5285 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5289 fNumberFitSuccess++;
5291 //_____________________________________________________________________________
5292 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5295 // Fit methode for the gain factor
5298 fCurrentCoef[0] = 0.0;
5299 fCurrentCoefE = 0.0;
5300 Double_t chisqrl = 0.0;
5301 Double_t chisqrg = 0.0;
5302 Double_t chisqr = 0.0;
5303 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5305 projch->Fit("landau","0",""
5306 ,(Double_t) mean/fBeginFitCharge
5307 ,projch->GetBinCenter(projch->GetNbinsX()));
5308 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5309 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5310 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5311 chisqrl = projch->GetFunction("landau")->GetChisquare();
5313 projch->Fit("gaus","0",""
5314 ,(Double_t) mean/fBeginFitCharge
5315 ,projch->GetBinCenter(projch->GetNbinsX()));
5316 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5317 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5318 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5320 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5321 if (fDebugLevel != 1) {
5322 projch->Fit("fLandauGaus","0",""
5323 ,(Double_t) mean/fBeginFitCharge
5324 ,projch->GetBinCenter(projch->GetNbinsX()));
5325 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5328 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5330 projch->Fit("fLandauGaus","+",""
5331 ,(Double_t) mean/fBeginFitCharge
5332 ,projch->GetBinCenter(projch->GetNbinsX()));
5333 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5335 fLandauGaus->Draw("same");
5338 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5339 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5340 fNumberFitSuccess++;
5341 CalculChargeCoefMean(kTRUE);
5342 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5343 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5346 CalculChargeCoefMean(kFALSE);
5347 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5350 if (fDebugLevel != 1) {
5355 //_____________________________________________________________________________
5356 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5359 // Fit methode for the gain factor more time consuming
5363 //Some parameters to initialise
5364 Double_t widthLandau, widthGaus, mPV, integral;
5365 Double_t chisquarel = 0.0;
5366 Double_t chisquareg = 0.0;
5367 projch->Fit("landau","0M+",""
5369 ,projch->GetBinCenter(projch->GetNbinsX()));
5370 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5371 chisquarel = projch->GetFunction("landau")->GetChisquare();
5372 projch->Fit("gaus","0M+",""
5374 ,projch->GetBinCenter(projch->GetNbinsX()));
5375 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5376 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5378 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5379 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5381 // Setting fit range and start values
5383 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5384 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5385 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5386 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5387 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5388 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5389 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5392 fCurrentCoef[0] = 0.0;
5393 fCurrentCoefE = 0.0;
5397 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5402 Double_t projchPeak;
5403 Double_t projchFWHM;
5404 LanGauPro(fp,projchPeak,projchFWHM);
5406 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5407 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5408 fNumberFitSuccess++;
5409 CalculChargeCoefMean(kTRUE);
5410 fCurrentCoef[0] = fp[1];
5411 fCurrentCoefE = fpe[1];
5412 //chargeCoefE2 = chisqr;
5415 CalculChargeCoefMean(kFALSE);
5416 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5418 if (fDebugLevel == 1) {
5419 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5420 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5423 fitsnr->Draw("same");
5429 //_____________________________________________________________________________
5430 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
5433 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5435 Double_t *c = new Double_t[5];
5436 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5437 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5438 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5443 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5444 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5451 //_____________________________________________________________________________
5452 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
5455 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5457 Double_t *c = new Double_t[5];
5458 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5459 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5460 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5461 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5465 c[2] = -(x0*(x[1]+x[2]+x[3])
5466 +x1*(x[0]+x[2]+x[3])
5467 +x2*(x[0]+x[1]+x[3])
5468 +x3*(x[0]+x[1]+x[2]));
5469 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5470 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5471 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5472 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5474 c[0] = -(x0*x[1]*x[2]*x[3]
5477 +x3*x[0]*x[1]*x[2]);
5485 //_____________________________________________________________________________
5486 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
5489 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5491 Double_t *c = new Double_t[5];
5492 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5493 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5494 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5495 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5496 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5499 c[4] = x0+x1+x2+x3+x4;
5500 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
5501 +x1*(x[0]+x[2]+x[3]+x[4])
5502 +x2*(x[0]+x[1]+x[3]+x[4])
5503 +x3*(x[0]+x[1]+x[2]+x[4])
5504 +x4*(x[0]+x[1]+x[2]+x[3]));
5505 c[2] = (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])
5506 +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])
5507 +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])
5508 +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])
5509 +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]));
5511 c[1] = -(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])
5512 +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])
5513 +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])
5514 +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])
5515 +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]));
5517 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
5518 +x1*x[0]*x[2]*x[3]*x[4]
5519 +x2*x[0]*x[1]*x[3]*x[4]
5520 +x3*x[0]*x[1]*x[2]*x[4]
5521 +x4*x[0]*x[1]*x[2]*x[3]);
5527 //_____________________________________________________________________________
5528 void AliTRDCalibraFit::NormierungCharge()
5531 // Normalisation of the gain factor resulting for the fits
5534 // Calcul of the mean of choosen method by fFitChargeNDB
5536 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5537 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5539 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5540 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5541 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5542 if (GetStack(detector) == 2) {
5545 if (GetStack(detector) != 2) {
5548 for (Int_t j = 0; j < total; j++) {
5556 fScaleFitFactor = fScaleFitFactor / sum;
5559 fScaleFitFactor = 1.0;
5562 //methode de boeuf mais bon...
5563 Double_t scalefactor = fScaleFitFactor;
5565 if(fDebugLevel > 1){
5567 if ( !fDebugStreamer ) {
5569 TDirectory *backup = gDirectory;
5570 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5571 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5573 (* fDebugStreamer) << "NormierungCharge"<<
5574 "scalefactor="<<scalefactor<<
5578 //_____________________________________________________________________________
5579 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5582 // Rebin of the 1D histo for the gain calibration if needed.
5583 // you have to choose fRebin, divider of fNumberBinCharge
5586 TAxis *xhist = hist->GetXaxis();
5587 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5588 ,xhist->GetBinLowEdge(1)
5589 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5591 AliInfo(Form("fRebin: %d",fRebin));
5593 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5595 for (Int_t ji = i; ji < i+fRebin; ji++) {
5596 sum += hist->GetBinContent(ji);
5599 rehist->SetBinContent(k,sum);
5607 //_____________________________________________________________________________
5608 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5611 // Rebin of the 1D histo for the gain calibration if needed
5612 // you have to choose fRebin divider of fNumberBinCharge
5615 TAxis *xhist = hist->GetXaxis();
5616 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5617 ,xhist->GetBinLowEdge(1)
5618 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5620 AliInfo(Form("fRebin: %d",fRebin));
5622 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5624 for (Int_t ji = i; ji < i+fRebin; ji++) {
5625 sum += hist->GetBinContent(ji);
5628 rehist->SetBinContent(k,sum);
5636 //_____________________________________________________________________________
5637 TH1F *AliTRDCalibraFit::CorrectTheError(const TGraphErrors *hist)
5640 // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
5641 // to be able to add them after
5642 // We convert it to a TH1F to be able to applied the same fit function method
5643 // After having called this function you can not add the statistics anymore
5648 Int_t nbins = hist->GetN();
5649 Double_t *x = hist->GetX();
5650 Double_t *entries = hist->GetEX();
5651 Double_t *mean = hist->GetY();
5652 Double_t *square = hist->GetEY();
5653 fEntriesCurrent = 0;
5659 Double_t step = x[1] - x[0];
5660 Double_t minvalue = x[0] - step/2;
5661 Double_t maxvalue = x[(nbins-1)] + step/2;
5663 rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
5665 for (Int_t k = 0; k < nbins; k++) {
5666 rehist->SetBinContent(k+1,mean[k]);
5667 if (entries[k] > 0.0) {
5668 fEntriesCurrent += (Int_t) entries[k];
5669 Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
5670 rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
5673 rehist->SetBinError(k+1,0.0);
5677 if(fEntriesCurrent > 0) fNumberEnt++;
5683 //____________Some basic geometry function_____________________________________
5686 //_____________________________________________________________________________
5687 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5690 // Reconstruct the plane number from the detector number
5693 return ((Int_t) (d % 6));
5697 //_____________________________________________________________________________
5698 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5701 // Reconstruct the stack number from the detector number
5703 const Int_t kNlayer = 6;
5705 return ((Int_t) (d % 30) / kNlayer);
5709 //_____________________________________________________________________________
5710 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5713 // Reconstruct the sector number from the detector number
5717 return ((Int_t) (d / fg));
5722 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5724 //_______________________________________________________________________________
5725 void AliTRDCalibraFit::ResetVectorFit()
5728 // Reset the VectorFits
5731 fVectorFit.SetOwner();
5733 fVectorFit2.SetOwner();
5734 fVectorFit2.Clear();
5738 //____________Private Functions________________________________________________
5741 //_____________________________________________________________________________
5742 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5745 // Function for the fit
5748 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5750 //PARAMETERS FOR FIT PH
5752 //fAsymmGauss->SetParameter(0,0.113755);
5753 //fAsymmGauss->SetParameter(1,0.350706);
5754 //fAsymmGauss->SetParameter(2,0.0604244);
5755 //fAsymmGauss->SetParameter(3,7.65596);
5756 //fAsymmGauss->SetParameter(4,1.00124);
5757 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5765 Double_t dx = 0.005;
5766 Double_t xs = par[1];
5768 Double_t paras[2] = { 0.0, 0.0 };
5771 if ((xs >= par[1]) &&
5772 (xs < (par[1]+par[2]))) {
5773 //fAsymmGauss->SetParameter(0,par[0]);
5774 //fAsymmGauss->SetParameter(1,xs);
5775 //ss += fAsymmGauss->Eval(xx);
5778 ss += AsymmGauss(&xx,paras);
5780 if ((xs >= (par[1]+par[2])) &&
5781 (xs < (par[1]+par[2]+par[3]))) {
5782 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5783 //fAsymmGauss->SetParameter(1,xs);
5784 //ss += fAsymmGauss->Eval(xx);
5785 paras[0] = par[0]*par[4];
5787 ss += AsymmGauss(&xx,paras);
5796 //_____________________________________________________________________________
5797 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5800 // Function for the fit
5803 //par[0] = normalization
5811 Double_t par1save = par[1];
5812 //Double_t par2save = par[2];
5813 Double_t par2save = 0.0604244;
5814 //Double_t par3save = par[3];
5815 Double_t par3save = 7.65596;
5816 //Double_t par5save = par[5];
5817 Double_t par5save = 0.870597;
5818 Double_t dx = x[0] - par1save;
5820 Double_t sigma2 = par2save*par2save;
5821 Double_t sqrt2 = TMath::Sqrt(2.0);
5822 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5823 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5824 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5825 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5827 //return par[0]*(exp1+par[4]*exp2);
5828 return par[0] * (exp1 + 1.00124 * exp2);
5832 //_____________________________________________________________________________
5833 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5836 // Sum Landau + Gaus with identical mean
5839 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5840 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5841 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5842 Double_t val = valLandau + valGaus;
5848 //_____________________________________________________________________________
5849 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
5852 // Function for the fit
5855 // par[0]=Width (scale) parameter of Landau density
5856 // par[1]=Most Probable (MP, location) parameter of Landau density
5857 // par[2]=Total area (integral -inf to inf, normalization constant)
5858 // par[3]=Width (sigma) of convoluted Gaussian function
5860 // In the Landau distribution (represented by the CERNLIB approximation),
5861 // the maximum is located at x=-0.22278298 with the location parameter=0.
5862 // This shift is corrected within this function, so that the actual
5863 // maximum is identical to the MP parameter.
5866 // Numeric constants
5867 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
5868 Double_t mpshift = -0.22278298; // Landau maximum location
5870 // Control constants
5871 Double_t np = 100.0; // Number of convolution steps
5872 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
5884 // MP shift correction
5885 mpc = par[1] - mpshift * par[0];
5887 // Range of convolution integral
5888 xlow = x[0] - sc * par[3];
5889 xupp = x[0] + sc * par[3];
5891 step = (xupp - xlow) / np;
5893 // Convolution integral of Landau and Gaussian by sum
5894 for (i = 1.0; i <= np/2; i++) {
5896 xx = xlow + (i-.5) * step;
5897 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5898 sum += fland * TMath::Gaus(x[0],xx,par[3]);
5900 xx = xupp - (i-.5) * step;
5901 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5902 sum += fland * TMath::Gaus(x[0],xx,par[3]);
5906 return (par[2] * step * sum * invsq2pi / par[3]);
5909 //_____________________________________________________________________________
5910 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
5911 , const Double_t *parlimitslo, const Double_t *parlimitshi
5912 , Double_t *fitparams, Double_t *fiterrors
5913 , Double_t *chiSqr, Int_t *ndf) const
5916 // Function for the fit
5920 Char_t funname[100];
5922 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
5927 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
5928 ffit->SetParameters(startvalues);
5929 ffit->SetParNames("Width","MP","Area","GSigma");
5931 for (i = 0; i < 4; i++) {
5932 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
5935 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
5937 ffit->GetParameters(fitparams); // Obtain fit parameters
5938 for (i = 0; i < 4; i++) {
5939 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
5941 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
5942 ndf[0] = ffit->GetNDF(); // Obtain ndf
5944 return (ffit); // Return fit function
5948 //_____________________________________________________________________________
5949 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
5952 // Function for the fit
5965 Int_t maxcalls = 10000;
5967 // Search for maximum
5968 p = params[1] - 0.1 * params[0];
5969 step = 0.05 * params[0];
5973 while ((l != lold) && (i < maxcalls)) {
5977 l = LanGauFun(&x,params);
5979 step = -step / 10.0;
5984 if (i == maxcalls) {
5990 // Search for right x location of fy
5991 p = maxx + params[0];
5997 while ( (l != lold) && (i < maxcalls) ) {
6002 l = TMath::Abs(LanGauFun(&x,params) - fy);
6016 // Search for left x location of fy
6018 p = maxx - 0.5 * params[0];
6024 while ((l != lold) && (i < maxcalls)) {
6028 l = TMath::Abs(LanGauFun(&x,params) - fy);
6030 step = -step / 10.0;
6035 if (i == maxcalls) {
6044 //_____________________________________________________________________________
6045 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6048 // Gaus with identical mean
6051 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];