1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
62 #include <TTreeStream.h>
63 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
71 #include "AliTRDCalibraFit.h"
72 #include "AliTRDCalibraMode.h"
73 #include "AliTRDCalibraVector.h"
74 #include "AliTRDCalibraVdriftLinearFit.h"
75 #include "AliTRDcalibDB.h"
76 #include "AliTRDgeometry.h"
77 #include "AliTRDpadPlane.h"
78 #include "AliTRDgeometry.h"
79 #include "AliTRDCommonParam.h"
80 #include "./Cal/AliTRDCalROC.h"
81 #include "./Cal/AliTRDCalPad.h"
82 #include "./Cal/AliTRDCalDet.h"
85 ClassImp(AliTRDCalibraFit)
87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
90 //_____________singleton implementation_________________________________________________
91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
94 // Singleton implementation
97 if (fgTerminated != kFALSE) {
101 if (fgInstance == 0) {
102 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
124 //______________________________________________________________________________________
125 AliTRDCalibraFit::AliTRDCalibraFit()
128 ,fNumberOfBinsExpected(0)
130 ,fBeginFitCharge(3.5)
132 ,fTakeTheMaxPH(kTRUE)
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 TString name = ch->GetTitle();
367 if(!SetModeCalibration(name,0)) return kFALSE;
369 // Number of Ybins (detectors or groups of pads)
370 Int_t nbins = ch->GetNbinsX();// charge
371 Int_t nybins = ch->GetNbinsY();// groups number
372 if (!InitFit(nybins,0)) {
378 fStatisticMean = 0.0;
380 fNumberFitSuccess = 0;
382 // Init fCountDet and fCount
383 InitfCountDetAndfCount(0);
384 // Beginning of the loop betwwen dect1 and dect2
385 for (Int_t idect = fDect1; idect < fDect2; idect++) {
386 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
387 UpdatefCountDetAndfCount(idect,0);
388 ReconstructFitRowMinRowMax(idect, 0);
390 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
391 projch->SetDirectory(0);
392 // Number of entries for this calibration group
393 Double_t nentries = 0.0;
395 for (Int_t k = 0; k < nbins; k++) {
396 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
397 nentries += ch->GetBinContent(binnb);
398 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
399 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
401 projch->SetEntries(nentries);
402 //printf("The number of entries for the group %d is %f\n",idect,nentries);
407 // Rebin and statistic stuff
409 projch = ReBin((TH1I *) projch);
411 // This detector has not enough statistics or was off
412 if (nentries <= fMinEntries) {
413 NotEnoughStatisticCH(idect);
414 if (fDebugLevel != 1) {
419 // Statistics of the group fitted
420 fStatisticMean += nentries;
425 case 0: FitMeanW((TH1 *) projch, nentries); break;
426 case 1: FitMean((TH1 *) projch, nentries, mean); break;
427 case 2: FitCH((TH1 *) projch, mean); break;
428 case 3: FitBisCH((TH1 *) projch, mean); break;
429 default: return kFALSE;
432 FillInfosFitCH(idect);
434 if (fDebugLevel != 1) {
439 if (fDebugLevel != 1) {
443 if (fNumberFit > 0) {
444 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
445 fStatisticMean = fStatisticMean / fNumberFit;
448 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
450 delete fDebugStreamer;
451 fDebugStreamer = 0x0;
455 //____________Functions fit Online CH2d________________________________________
456 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
459 // Reconstruct a 1D histo from the vectorCH for each calibration group,
460 // fit the histo, normalized the resulted coefficients (to 1 normally)
463 // Set the calibraMode
464 //const char *name = calvect->GetNameCH();
465 TString name = calvect->GetNameCH();
466 if(!SetModeCalibration(name,0)) return kFALSE;
468 // Number of Xbins (detectors or groups of pads)
469 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
475 fStatisticMean = 0.0;
477 fNumberFitSuccess = 0;
479 // Init fCountDet and fCount
480 InitfCountDetAndfCount(0);
481 // Beginning of the loop between dect1 and dect2
482 for (Int_t idect = fDect1; idect < fDect2; idect++) {
483 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
484 UpdatefCountDetAndfCount(idect,0);
485 ReconstructFitRowMinRowMax(idect,0);
487 Double_t nentries = 0.0;
489 if(!calvect->GetCHEntries(fCountDet)) {
490 NotEnoughStatisticCH(idect);
496 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
497 projch->SetDirectory(0);
498 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
499 nentries += projch->GetBinContent(k+1);
500 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
506 //printf("The number of entries for the group %d is %f\n",idect,nentries);
509 projch = ReBin((TH1F *) projch);
511 // This detector has not enough statistics or was not found in VectorCH
512 if (nentries <= fMinEntries) {
513 NotEnoughStatisticCH(idect);
516 // Statistic of the histos fitted
517 fStatisticMean += nentries;
522 case 0: FitMeanW((TH1 *) projch, nentries); break;
523 case 1: FitMean((TH1 *) projch, nentries, mean); break;
524 case 2: FitCH((TH1 *) projch, mean); break;
525 case 3: FitBisCH((TH1 *) projch, mean); break;
526 default: return kFALSE;
529 FillInfosFitCH(idect);
532 if (fDebugLevel != 1) {
536 if (fNumberFit > 0) {
537 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
538 fStatisticMean = fStatisticMean / fNumberFit;
541 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
543 delete fDebugStreamer;
544 fDebugStreamer = 0x0;
547 //________________functions fit Online PH2d____________________________________
548 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
551 // Take the 1D profiles (average pulse height), projections of the 2D PH
552 // on the Xaxis, for each calibration group
553 // Reconstruct a drift velocity
554 // A first calibration of T0 is also made using the same method
557 // Set the calibration mode
558 //const char *name = ph->GetTitle();
559 TString name = ph->GetTitle();
560 if(!SetModeCalibration(name,1)) return kFALSE;
562 //printf("Mode calibration set\n");
564 // Number of Xbins (detectors or groups of pads)
565 Int_t nbins = ph->GetNbinsX();// time
566 Int_t nybins = ph->GetNbinsY();// calibration group
567 if (!InitFit(nybins,1)) {
571 //printf("Init fit\n");
577 //printf("Init fit PH\n");
579 fStatisticMean = 0.0;
581 fNumberFitSuccess = 0;
583 // Init fCountDet and fCount
584 InitfCountDetAndfCount(1);
585 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
587 // Beginning of the loop
588 for (Int_t idect = fDect1; idect < fDect2; idect++) {
589 //printf("idect = %d\n",idect);
590 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
591 UpdatefCountDetAndfCount(idect,1);
592 ReconstructFitRowMinRowMax(idect,1);
594 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
595 projph->SetDirectory(0);
596 // Number of entries for this calibration group
597 Double_t nentries = 0;
598 for (Int_t k = 0; k < nbins; k++) {
599 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
600 nentries += ph->GetBinEntries(binnb);
605 //printf("The number of entries for the group %d is %f\n",idect,nentries);
606 // This detector has not enough statistics or was off
607 if (nentries <= fMinEntries) {
608 //printf("Not enough statistic!\n");
609 NotEnoughStatisticPH(idect,nentries);
610 if (fDebugLevel != 1) {
615 // Statistics of the histos fitted
617 fStatisticMean += nentries;
618 // Calcul of "real" coef
619 CalculVdriftCoefMean();
622 //printf("Method\n");
625 case 0: FitLagrangePoly((TH1 *) projph); break;
626 case 1: FitPente((TH1 *) projph); break;
627 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
628 default: return kFALSE;
630 // Fill the tree if end of a detector or only the pointer to the branch!!!
631 FillInfosFitPH(idect,nentries);
633 if (fDebugLevel != 1) {
638 if (fNumberFit > 0) {
639 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
640 fStatisticMean = fStatisticMean / fNumberFit;
643 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
645 delete fDebugStreamer;
646 fDebugStreamer = 0x0;
649 //____________Functions fit Online PH2d________________________________________
650 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
653 // Reconstruct the average pulse height from the vectorPH for each
655 // Reconstruct a drift velocity
656 // A first calibration of T0 is also made using the same method (slope method)
659 // Set the calibration mode
660 //const char *name = calvect->GetNamePH();
661 TString name = calvect->GetNamePH();
662 if(!SetModeCalibration(name,1)) return kFALSE;
664 // Number of Xbins (detectors or groups of pads)
665 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
671 fStatisticMean = 0.0;
673 fNumberFitSuccess = 0;
675 // Init fCountDet and fCount
676 InitfCountDetAndfCount(1);
677 // Beginning of the loop
678 for (Int_t idect = fDect1; idect < fDect2; idect++) {
679 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
680 UpdatefCountDetAndfCount(idect,1);
681 ReconstructFitRowMinRowMax(idect,1);
684 if(!calvect->GetPHEntries(fCountDet)) {
685 NotEnoughStatisticPH(idect,fEntriesCurrent);
690 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
691 projph->SetDirectory(0);
692 if(fEntriesCurrent > 0) fNumberEnt++;
693 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
694 // This detector has not enough statistics or was off
695 if (fEntriesCurrent <= fMinEntries) {
696 //printf("Not enough stat!\n");
697 NotEnoughStatisticPH(idect,fEntriesCurrent);
700 // Statistic of the histos fitted
702 fStatisticMean += fEntriesCurrent;
703 // Calcul of "real" coef
704 CalculVdriftCoefMean();
709 case 0: FitLagrangePoly((TH1 *) projph); break;
710 case 1: FitPente((TH1 *) projph); break;
711 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
712 default: return kFALSE;
714 // Fill the tree if end of a detector or only the pointer to the branch!!!
715 FillInfosFitPH(idect,fEntriesCurrent);
719 if (fNumberFit > 0) {
720 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
721 fStatisticMean = fStatisticMean / fNumberFit;
724 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
726 delete fDebugStreamer;
727 fDebugStreamer = 0x0;
730 //____________Functions fit Online PRF2d_______________________________________
731 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
734 // Take the 1D profiles (pad response function), projections of the 2D PRF
735 // on the Xaxis, for each calibration group
736 // Fit with a gaussian to reconstruct the sigma of the pad response function
739 // Set the calibration mode
740 //const char *name = prf->GetTitle();
741 TString name = prf->GetTitle();
742 if(!SetModeCalibration(name,2)) return kFALSE;
744 // Number of Ybins (detectors or groups of pads)
745 Int_t nybins = prf->GetNbinsY();// calibration groups
746 Int_t nbins = prf->GetNbinsX();// bins
747 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
748 if((nbg > 0) || (nbg == -1)) return kFALSE;
749 if (!InitFit(nybins,2)) {
755 fStatisticMean = 0.0;
757 fNumberFitSuccess = 0;
759 // Init fCountDet and fCount
760 InitfCountDetAndfCount(2);
761 // Beginning of the loop
762 for (Int_t idect = fDect1; idect < fDect2; idect++) {
763 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
764 UpdatefCountDetAndfCount(idect,2);
765 ReconstructFitRowMinRowMax(idect,2);
767 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
768 projprf->SetDirectory(0);
769 // Number of entries for this calibration group
770 Double_t nentries = 0;
771 for (Int_t k = 0; k < nbins; k++) {
772 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
773 nentries += prf->GetBinEntries(binnb);
775 if(nentries > 0) fNumberEnt++;
776 // This detector has not enough statistics or was off
777 if (nentries <= fMinEntries) {
778 NotEnoughStatisticPRF(idect);
779 if (fDebugLevel != 1) {
784 // Statistics of the histos fitted
786 fStatisticMean += nentries;
787 // Calcul of "real" coef
792 case 0: FitPRF((TH1 *) projprf); break;
793 case 1: RmsPRF((TH1 *) projprf); break;
794 default: return kFALSE;
796 // Fill the tree if end of a detector or only the pointer to the branch!!!
797 FillInfosFitPRF(idect);
799 if (fDebugLevel != 1) {
804 if (fNumberFit > 0) {
805 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
806 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
807 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
808 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
809 fStatisticMean = fStatisticMean / fNumberFit;
812 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
814 delete fDebugStreamer;
815 fDebugStreamer = 0x0;
818 //____________Functions fit Online PRF2d_______________________________________
819 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
822 // Take the 1D profiles (pad response function), projections of the 2D PRF
823 // on the Xaxis, for each calibration group
824 // Fit with a gaussian to reconstruct the sigma of the pad response function
827 // Set the calibration mode
828 //const char *name = prf->GetTitle();
829 TString name = prf->GetTitle();
830 if(!SetModeCalibration(name,2)) return kFALSE;
832 // Number of Ybins (detectors or groups of pads)
833 TAxis *xprf = prf->GetXaxis();
834 TAxis *yprf = prf->GetYaxis();
835 Int_t nybins = yprf->GetNbins();// calibration groups
836 Int_t nbins = xprf->GetNbins();// bins
837 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
838 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
839 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
840 if(nbg == -1) return kFALSE;
841 if(nbg > 0) fMethod = 1;
843 if (!InitFit(nybins,2)) {
849 fStatisticMean = 0.0;
851 fNumberFitSuccess = 0;
853 // Init fCountDet and fCount
854 InitfCountDetAndfCount(2);
855 // Beginning of the loop
856 for (Int_t idect = fDect1; idect < fDect2; idect++) {
857 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
858 UpdatefCountDetAndfCount(idect,2);
859 ReconstructFitRowMinRowMax(idect,2);
860 // Build the array of entries and sum
861 TArrayD arraye = TArrayD(nbins);
862 TArrayD arraym = TArrayD(nbins);
863 TArrayD arrayme = TArrayD(nbins);
864 Double_t nentries = 0;
865 //printf("nbins %d\n",nbins);
866 for (Int_t k = 0; k < nbins; k++) {
867 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
868 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
869 Double_t mean = (Double_t)prf->GetBinContent(binnb);
870 Double_t error = (Double_t)prf->GetBinError(binnb);
871 //printf("for %d we have %f\n",k,entries);
873 arraye.AddAt(entries,k);
874 arraym.AddAt(mean,k);
875 arrayme.AddAt(error,k);
877 if(nentries > 0) fNumberEnt++;
878 //printf("The number of entries for the group %d is %f\n",idect,nentries);
879 // This detector has not enough statistics or was off
880 if (nentries <= fMinEntries) {
881 NotEnoughStatisticPRF(idect);
884 // Statistics of the histos fitted
886 fStatisticMean += nentries;
887 // Calcul of "real" coef
892 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
893 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
894 default: return kFALSE;
896 // Fill the tree if end of a detector or only the pointer to the branch!!!
897 FillInfosFitPRF(idect);
900 if (fNumberFit > 0) {
901 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
902 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
903 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
904 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
905 fStatisticMean = fStatisticMean / fNumberFit;
908 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
910 delete fDebugStreamer;
911 fDebugStreamer = 0x0;
914 //____________Functions fit Online PRF2d_______________________________________
915 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
918 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
919 // each calibration group
920 // Fit with a gaussian to reconstruct the sigma of the pad response function
923 // Set the calibra mode
924 //const char *name = calvect->GetNamePRF();
925 TString name = calvect->GetNamePRF();
926 if(!SetModeCalibration(name,2)) return kFALSE;
927 //printf("test0 %s\n",name);
929 // Number of Xbins (detectors or groups of pads)
930 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
935 ///printf("test2\n");
938 fStatisticMean = 0.0;
940 fNumberFitSuccess = 0;
942 // Init fCountDet and fCount
943 InitfCountDetAndfCount(2);
944 // Beginning of the loop
945 for (Int_t idect = fDect1; idect < fDect2; idect++) {
946 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
947 UpdatefCountDetAndfCount(idect,2);
948 ReconstructFitRowMinRowMax(idect,2);
951 if(!calvect->GetPRFEntries(fCountDet)) {
952 NotEnoughStatisticPRF(idect);
955 TString tname("PRF");
957 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
958 projprf->SetDirectory(0);
959 if(fEntriesCurrent > 0) fNumberEnt++;
960 // This detector has not enough statistics or was off
961 if (fEntriesCurrent <= fMinEntries) {
962 NotEnoughStatisticPRF(idect);
965 // Statistic of the histos fitted
967 fStatisticMean += fEntriesCurrent;
968 // Calcul of "real" coef
973 case 1: FitPRF((TH1 *) projprf); break;
974 case 2: RmsPRF((TH1 *) projprf); break;
975 default: return kFALSE;
977 // Fill the tree if end of a detector or only the pointer to the branch!!!
978 FillInfosFitPRF(idect);
981 if (fNumberFit > 0) {
982 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
985 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
987 delete fDebugStreamer;
988 fDebugStreamer = 0x0;
991 //____________Functions fit Online PRF2d_______________________________________
992 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
995 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
996 // each calibration group
997 // Fit with a gaussian to reconstruct the sigma of the pad response function
1000 // Set the calibra mode
1001 //const char *name = calvect->GetNamePRF();
1002 TString name = calvect->GetNamePRF();
1003 if(!SetModeCalibration(name,2)) return kFALSE;
1004 //printf("test0 %s\n",name);
1005 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1006 //printf("test1 %d\n",nbg);
1007 if(nbg == -1) return kFALSE;
1008 if(nbg > 0) fMethod = 1;
1010 // Number of Xbins (detectors or groups of pads)
1011 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1012 //printf("test2\n");
1015 if (!InitFitPRF()) {
1016 //printf("test3\n");
1019 fStatisticMean = 0.0;
1021 fNumberFitSuccess = 0;
1025 Double_t *arrayx = 0;
1026 Double_t *arraye = 0;
1027 Double_t *arraym = 0;
1028 Double_t *arrayme = 0;
1029 Float_t lowedge = 0.0;
1030 Float_t upedge = 0.0;
1031 // Init fCountDet and fCount
1032 InitfCountDetAndfCount(2);
1033 // Beginning of the loop
1034 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1035 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1036 UpdatefCountDetAndfCount(idect,2);
1037 ReconstructFitRowMinRowMax(idect,2);
1039 fEntriesCurrent = 0;
1040 if(!calvect->GetPRFEntries(fCountDet)) {
1041 NotEnoughStatisticPRF(idect);
1044 TString tname("PRF");
1046 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1047 nbins = projprftree->GetN();
1048 arrayx = (Double_t *)projprftree->GetX();
1049 arraye = (Double_t *)projprftree->GetEX();
1050 arraym = (Double_t *)projprftree->GetY();
1051 arrayme = (Double_t *)projprftree->GetEY();
1052 Float_t step = arrayx[1]-arrayx[0];
1053 lowedge = arrayx[0] - step/2.0;
1054 upedge = arrayx[(nbins-1)] + step/2.0;
1055 //printf("nbins est %d\n",nbins);
1056 for(Int_t k = 0; k < nbins; k++){
1057 fEntriesCurrent += (Int_t)arraye[k];
1058 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1059 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1061 if(fEntriesCurrent > 0) fNumberEnt++;
1062 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1063 // This detector has not enough statistics or was off
1064 if (fEntriesCurrent <= fMinEntries) {
1065 NotEnoughStatisticPRF(idect);
1068 // Statistic of the histos fitted
1070 fStatisticMean += fEntriesCurrent;
1071 // Calcul of "real" coef
1072 CalculPRFCoefMean();
1076 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1077 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1078 default: return kFALSE;
1080 // Fill the tree if end of a detector or only the pointer to the branch!!!
1081 FillInfosFitPRF(idect);
1084 if (fNumberFit > 0) {
1085 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1088 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1090 delete fDebugStreamer;
1091 fDebugStreamer = 0x0;
1094 //____________Functions fit Online CH2d________________________________________
1095 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1098 // The linear method
1101 fStatisticMean = 0.0;
1103 fNumberFitSuccess = 0;
1105 if(!InitFitLinearFitter()) return kFALSE;
1108 for(Int_t idet = 0; idet < 540; idet++){
1111 //printf("detector number %d\n",idet);
1116 fEntriesCurrent = 0;
1118 Bool_t here = calivdli->GetParam(idet,¶m);
1119 Bool_t heree = calivdli->GetError(idet,&error);
1120 //printf("here %d and heree %d\n",here, heree);
1122 fEntriesCurrent = (Int_t) error[2];
1125 //printf("Number of entries %d\n",fEntriesCurrent);
1126 // Nothing found or not enough statistic
1127 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1128 NotEnoughStatisticLinearFitter();
1135 fStatisticMean += fEntriesCurrent;
1138 if((-(param[1])) <= 0.0) {
1139 NotEnoughStatisticLinearFitter();
1143 // CalculDatabaseVdriftandTan
1144 CalculVdriftLorentzCoef();
1147 fNumberFitSuccess ++;
1149 // Put the fCurrentCoef
1150 fCurrentCoef[0] = -param[1];
1151 // here the database must be the one of the reconstruction for the lorentz angle....
1152 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1153 fCurrentCoefE = error[1];
1154 fCurrentCoefE2 = error[0];
1155 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1156 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1160 FillInfosFitLinearFitter();
1165 if (fNumberFit > 0) {
1166 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1169 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1171 delete fDebugStreamer;
1172 fDebugStreamer = 0x0;
1176 //____________Functions for seeing if the pad is really okey___________________
1177 //_____________________________________________________________________________
1178 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1181 // Get numberofgroupsprf
1185 const Char_t *pattern0 = "Ngp0";
1186 const Char_t *pattern1 = "Ngp1";
1187 const Char_t *pattern2 = "Ngp2";
1188 const Char_t *pattern3 = "Ngp3";
1189 const Char_t *pattern4 = "Ngp4";
1190 const Char_t *pattern5 = "Ngp5";
1191 const Char_t *pattern6 = "Ngp6";
1194 if (strstr(nametitle.Data(),pattern0)) {
1197 if (strstr(nametitle.Data(),pattern1)) {
1200 if (strstr(nametitle.Data(),pattern2)) {
1203 if (strstr(nametitle.Data(),pattern3)) {
1206 if (strstr(nametitle.Data(),pattern4)) {
1209 if (strstr(nametitle.Data(),pattern5)) {
1212 if (strstr(nametitle.Data(),pattern6)){
1219 //_____________________________________________________________________________
1220 Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1223 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1224 // corresponding to the given name
1227 if(!SetNzFromTObject(name,i)) return kFALSE;
1228 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1233 //_____________________________________________________________________________
1234 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1237 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1238 // corresponding to the given TObject
1242 const Char_t *patternrphi0 = "Nrphi0";
1243 const Char_t *patternrphi1 = "Nrphi1";
1244 const Char_t *patternrphi2 = "Nrphi2";
1245 const Char_t *patternrphi3 = "Nrphi3";
1246 const Char_t *patternrphi4 = "Nrphi4";
1247 const Char_t *patternrphi5 = "Nrphi5";
1248 const Char_t *patternrphi6 = "Nrphi6";
1251 const Char_t *patternrphi10 = "Nrphi10";
1252 const Char_t *patternrphi100 = "Nrphi100";
1253 const Char_t *patternz10 = "Nz10";
1254 const Char_t *patternz100 = "Nz100";
1257 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1258 fCalibraMode->SetAllTogether(i);
1260 if (fDebugLevel > 1) {
1261 AliInfo(Form("fNbDet %d and 100",fNbDet));
1265 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1266 fCalibraMode->SetPerSuperModule(i);
1268 if (fDebugLevel > 1) {
1269 AliInfo(Form("fNDet %d and 100",fNbDet));
1274 if (strstr(name.Data(),patternrphi0)) {
1275 fCalibraMode->SetNrphi(i ,0);
1276 if (fDebugLevel > 1) {
1277 AliInfo(Form("fNbDet %d and 0",fNbDet));
1281 if (strstr(name.Data(),patternrphi1)) {
1282 fCalibraMode->SetNrphi(i, 1);
1283 if (fDebugLevel > 1) {
1284 AliInfo(Form("fNbDet %d and 1",fNbDet));
1288 if (strstr(name.Data(),patternrphi2)) {
1289 fCalibraMode->SetNrphi(i, 2);
1290 if (fDebugLevel > 1) {
1291 AliInfo(Form("fNbDet %d and 2",fNbDet));
1295 if (strstr(name.Data(),patternrphi3)) {
1296 fCalibraMode->SetNrphi(i, 3);
1297 if (fDebugLevel > 1) {
1298 AliInfo(Form("fNbDet %d and 3",fNbDet));
1302 if (strstr(name.Data(),patternrphi4)) {
1303 fCalibraMode->SetNrphi(i, 4);
1304 if (fDebugLevel > 1) {
1305 AliInfo(Form("fNbDet %d and 4",fNbDet));
1309 if (strstr(name.Data(),patternrphi5)) {
1310 fCalibraMode->SetNrphi(i, 5);
1311 if (fDebugLevel > 1) {
1312 AliInfo(Form("fNbDet %d and 5",fNbDet));
1316 if (strstr(name.Data(),patternrphi6)) {
1317 fCalibraMode->SetNrphi(i, 6);
1318 if (fDebugLevel > 1) {
1319 AliInfo(Form("fNbDet %d and 6",fNbDet));
1324 if (fDebugLevel > 1) {
1325 AliInfo(Form("fNbDet %d and rest",fNbDet));
1327 fCalibraMode->SetNrphi(i ,0);
1331 //_____________________________________________________________________________
1332 Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1335 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1336 // corresponding to the given TObject
1340 const Char_t *patternz0 = "Nz0";
1341 const Char_t *patternz1 = "Nz1";
1342 const Char_t *patternz2 = "Nz2";
1343 const Char_t *patternz3 = "Nz3";
1344 const Char_t *patternz4 = "Nz4";
1346 const Char_t *patternrphi10 = "Nrphi10";
1347 const Char_t *patternrphi100 = "Nrphi100";
1348 const Char_t *patternz10 = "Nz10";
1349 const Char_t *patternz100 = "Nz100";
1351 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1352 fCalibraMode->SetAllTogether(i);
1354 if (fDebugLevel > 1) {
1355 AliInfo(Form("fNbDet %d and 100",fNbDet));
1359 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1360 fCalibraMode->SetPerSuperModule(i);
1362 if (fDebugLevel > 1) {
1363 AliInfo(Form("fNbDet %d and 10",fNbDet));
1367 if (strstr(name.Data(),patternz0)) {
1368 fCalibraMode->SetNz(i, 0);
1369 if (fDebugLevel > 1) {
1370 AliInfo(Form("fNbDet %d and 0",fNbDet));
1374 if (strstr(name.Data(),patternz1)) {
1375 fCalibraMode->SetNz(i ,1);
1376 if (fDebugLevel > 1) {
1377 AliInfo(Form("fNbDet %d and 1",fNbDet));
1381 if (strstr(name.Data(),patternz2)) {
1382 fCalibraMode->SetNz(i ,2);
1383 if (fDebugLevel > 1) {
1384 AliInfo(Form("fNbDet %d and 2",fNbDet));
1388 if (strstr(name.Data(),patternz3)) {
1389 fCalibraMode->SetNz(i ,3);
1390 if (fDebugLevel > 1) {
1391 AliInfo(Form("fNbDet %d and 3",fNbDet));
1395 if (strstr(name.Data(),patternz4)) {
1396 fCalibraMode->SetNz(i ,4);
1397 if (fDebugLevel > 1) {
1398 AliInfo(Form("fNbDet %d and 4",fNbDet));
1403 if (fDebugLevel > 1) {
1404 AliInfo(Form("fNbDet %d and rest",fNbDet));
1406 fCalibraMode->SetNz(i ,0);
1409 //______________________________________________________________________
1410 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1412 // Remove the results too far from the mean value and rms
1413 // type: 0 gain, 1 vdrift
1417 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1419 AliInfo("The Vector Fit is not complete!");
1422 Int_t detector = -1;
1424 Float_t value = 0.0;
1426 /////////////////////////////////
1427 // Calculate the mean values
1428 ////////////////////////////////
1430 ////////////////////////
1431 Double_t meanAll = 0.0;
1432 Double_t rmsAll = 0.0;
1437 for (Int_t k = 0; k < loop; k++) {
1438 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1439 sector = GetSector(detector);
1441 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1443 rmsAll += value*value;
1449 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1450 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1451 for (Int_t row = 0; row < rowMax; row++) {
1452 for (Int_t col = 0; col < colMax; col++) {
1453 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1455 rmsAll += value*value;
1465 meanAll = meanAll/countAll;
1466 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1468 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1469 /////////////////////////////////////////////////
1471 ////////////////////////////////////////////////
1472 Double_t defaultvalue = -1.0;
1473 if(type==1) defaultvalue = -1.5;
1474 for (Int_t k = 0; k < loop; k++) {
1475 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1476 sector = GetSector(detector);
1477 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1478 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1479 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1481 // remove the results too far away
1482 for (Int_t row = 0; row < rowMax; row++) {
1483 for (Int_t col = 0; col < colMax; col++) {
1484 value = coef[(Int_t)(col*rowMax+row)];
1485 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1486 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1492 //______________________________________________________________________
1493 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1495 // Remove the results too far from the mean and rms
1499 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1501 AliInfo("The Vector Fit is not complete!");
1504 Int_t detector = -1;
1506 Float_t value = 0.0;
1508 /////////////////////////////////
1509 // Calculate the mean values
1510 ////////////////////////////////
1512 ////////////////////////
1513 Double_t meanAll = 0.0;
1514 Double_t rmsAll = 0.0;
1519 for (Int_t k = 0; k < loop; k++) {
1520 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1521 sector = GetSector(detector);
1523 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1526 rmsAll += value*value;
1531 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1532 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1533 for (Int_t row = 0; row < rowMax; row++) {
1534 for (Int_t col = 0; col < colMax; col++) {
1535 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1537 rmsAll += value*value;
1546 meanAll = meanAll/countAll;
1547 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1549 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1550 /////////////////////////////////////////////////
1552 ////////////////////////////////////////////////
1553 for (Int_t k = 0; k < loop; k++) {
1554 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1555 sector = GetSector(detector);
1556 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1557 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1558 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1560 // remove the results too far away
1561 for (Int_t row = 0; row < rowMax; row++) {
1562 for (Int_t col = 0; col < colMax; col++) {
1563 value = coef[(Int_t)(col*rowMax+row)];
1564 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) coef[(Int_t)(col*rowMax+row)] = 100.0;
1569 //______________________________________________________________________
1570 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1572 // ofwhat is equaled to 0: mean value of all passing detectors
1573 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1576 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1578 AliInfo("The Vector Fit is not complete!");
1581 Int_t detector = -1;
1583 Float_t value = 0.0;
1585 /////////////////////////////////
1586 // Calculate the mean values
1587 ////////////////////////////////
1589 ////////////////////////
1590 Double_t meanAll = 0.0;
1591 Double_t meanSupermodule[18];
1592 Double_t meanDetector[540];
1593 Double_t rmsAll = 0.0;
1594 Double_t rmsSupermodule[18];
1595 Double_t rmsDetector[540];
1597 Int_t countSupermodule[18];
1598 Int_t countDetector[540];
1599 for(Int_t sm = 0; sm < 18; sm++){
1600 rmsSupermodule[sm] = 0.0;
1601 meanSupermodule[sm] = 0.0;
1602 countSupermodule[sm] = 0;
1604 for(Int_t det = 0; det < 540; det++){
1605 rmsDetector[det] = 0.0;
1606 meanDetector[det] = 0.0;
1607 countDetector[det] = 0;
1612 for (Int_t k = 0; k < loop; k++) {
1613 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1614 sector = GetSector(detector);
1616 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1618 rmsDetector[detector] += value*value;
1619 meanDetector[detector] += value;
1620 countDetector[detector]++;
1621 rmsSupermodule[sector] += value*value;
1622 meanSupermodule[sector] += value;
1623 countSupermodule[sector]++;
1624 rmsAll += value*value;
1630 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1631 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1632 for (Int_t row = 0; row < rowMax; row++) {
1633 for (Int_t col = 0; col < colMax; col++) {
1634 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1636 rmsDetector[detector] += value*value;
1637 meanDetector[detector] += value;
1638 countDetector[detector]++;
1639 rmsSupermodule[sector] += value*value;
1640 meanSupermodule[sector] += value;
1641 countSupermodule[sector]++;
1642 rmsAll += value*value;
1652 meanAll = meanAll/countAll;
1653 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1655 for(Int_t sm = 0; sm < 18; sm++){
1656 if(countSupermodule[sm] > 0) {
1657 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1658 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1661 for(Int_t det = 0; det < 540; det++){
1662 if(countDetector[det] > 0) {
1663 meanDetector[det] = meanDetector[det]/countDetector[det];
1664 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1667 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1668 ///////////////////////////////////////////////
1669 // Put the mean value for the no-fitted
1670 /////////////////////////////////////////////
1671 for (Int_t k = 0; k < loop; k++) {
1672 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1673 sector = GetSector(detector);
1674 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1675 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1676 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1678 for (Int_t row = 0; row < rowMax; row++) {
1679 for (Int_t col = 0; col < colMax; col++) {
1680 value = coef[(Int_t)(col*rowMax+row)];
1682 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1684 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1685 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1686 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1690 if(fDebugLevel > 1){
1692 if ( !fDebugStreamer ) {
1694 TDirectory *backup = gDirectory;
1695 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1696 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1699 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1701 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1702 "detector="<<detector<<
1714 //______________________________________________________________________
1715 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1717 // ofwhat is equaled to 0: mean value of all passing detectors
1718 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1721 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1723 AliInfo("The Vector Fit is not complete!");
1726 Int_t detector = -1;
1728 Float_t value = 0.0;
1730 /////////////////////////////////
1731 // Calculate the mean values
1732 ////////////////////////////////
1734 ////////////////////////
1735 Double_t meanAll = 0.0;
1736 Double_t rmsAll = 0.0;
1737 Double_t meanSupermodule[18];
1738 Double_t rmsSupermodule[18];
1739 Double_t meanDetector[540];
1740 Double_t rmsDetector[540];
1742 Int_t countSupermodule[18];
1743 Int_t countDetector[540];
1744 for(Int_t sm = 0; sm < 18; sm++){
1745 rmsSupermodule[sm] = 0.0;
1746 meanSupermodule[sm] = 0.0;
1747 countSupermodule[sm] = 0;
1749 for(Int_t det = 0; det < 540; det++){
1750 rmsDetector[det] = 0.0;
1751 meanDetector[det] = 0.0;
1752 countDetector[det] = 0;
1756 for (Int_t k = 0; k < loop; k++) {
1757 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1758 sector = GetSector(detector);
1760 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1762 rmsDetector[detector] += value*value;
1763 meanDetector[detector] += value;
1764 countDetector[detector]++;
1765 rmsSupermodule[sector] += value*value;
1766 meanSupermodule[sector] += value;
1767 countSupermodule[sector]++;
1769 rmsAll += value*value;
1774 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1775 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1776 for (Int_t row = 0; row < rowMax; row++) {
1777 for (Int_t col = 0; col < colMax; col++) {
1778 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1780 rmsDetector[detector] += value*value;
1781 meanDetector[detector] += value;
1782 countDetector[detector]++;
1783 rmsSupermodule[sector] += value*value;
1784 meanSupermodule[sector] += value;
1785 countSupermodule[sector]++;
1786 rmsAll += value*value;
1796 meanAll = meanAll/countAll;
1797 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1799 for(Int_t sm = 0; sm < 18; sm++){
1800 if(countSupermodule[sm] > 0) {
1801 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1802 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1805 for(Int_t det = 0; det < 540; det++){
1806 if(countDetector[det] > 0) {
1807 meanDetector[det] = meanDetector[det]/countDetector[det];
1808 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1811 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1812 ////////////////////////////////////////////
1813 // Put the mean value for the no-fitted
1814 /////////////////////////////////////////////
1815 for (Int_t k = 0; k < loop; k++) {
1816 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1817 sector = GetSector(detector);
1818 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1819 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1820 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1822 for (Int_t row = 0; row < rowMax; row++) {
1823 for (Int_t col = 0; col < colMax; col++) {
1824 value = coef[(Int_t)(col*rowMax+row)];
1826 if((ofwhat == 0) && (meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1828 if((meanDetector[detector] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1829 else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1830 else if((meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1834 if(fDebugLevel > 1){
1836 if ( !fDebugStreamer ) {
1838 TDirectory *backup = gDirectory;
1839 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1840 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1843 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1845 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1846 "detector="<<detector<<
1859 //_____________________________________________________________________________
1860 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1863 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1864 // It takes the mean value of the coefficients per detector
1865 // This object has to be written in the database
1868 // Create the DetObject
1869 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1871 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1872 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1873 Int_t detector = -1;
1874 Float_t value = 0.0;
1877 for (Int_t k = 0; k < loop; k++) {
1878 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1881 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1885 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1886 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1887 for (Int_t row = 0; row < rowMax; row++) {
1888 for (Int_t col = 0; col < colMax; col++) {
1889 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1890 mean += TMath::Abs(value);
1894 if(count > 0) mean = mean/count;
1896 object->SetValue(detector,mean);
1901 //_____________________________________________________________________________
1902 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1905 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1906 // It takes the mean value of the coefficients per detector
1907 // This object has to be written in the database
1910 // Create the DetObject
1911 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1914 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1915 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1916 Int_t detector = -1;
1917 Float_t value = 0.0;
1919 for (Int_t k = 0; k < loop; k++) {
1920 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1923 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1924 if(!meanOtherBefore){
1925 if(value > 0) value = value*scaleFitFactor;
1927 else value = value*scaleFitFactor;
1928 mean = TMath::Abs(value);
1932 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1933 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1934 for (Int_t row = 0; row < rowMax; row++) {
1935 for (Int_t col = 0; col < colMax; col++) {
1936 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1937 if(!meanOtherBefore) {
1938 if(value > 0) value = value*scaleFitFactor;
1940 else value = value*scaleFitFactor;
1941 mean += TMath::Abs(value);
1945 if(count > 0) mean = mean/count;
1947 if(mean < 0.1) mean = 0.1;
1948 object->SetValue(detector,mean);
1953 //_____________________________________________________________________________
1954 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1957 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1958 // It takes the min value of the coefficients per detector
1959 // This object has to be written in the database
1962 // Create the DetObject
1963 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1965 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1966 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1967 Int_t detector = -1;
1968 Float_t value = 0.0;
1970 for (Int_t k = 0; k < loop; k++) {
1971 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1972 Float_t min = 100.0;
1974 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1976 if(value > 70.0) value = value-100.0;
1981 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1982 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1983 for (Int_t row = 0; row < rowMax; row++) {
1984 for (Int_t col = 0; col < colMax; col++) {
1985 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1987 if(value > 70.0) value = value-100.0;
1989 if(min > value) min = value;
1993 object->SetValue(detector,min);
1999 //_____________________________________________________________________________
2000 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2003 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2004 // It takes the min value of the coefficients per detector
2005 // This object has to be written in the database
2008 // Create the DetObject
2009 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2012 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2013 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2014 Int_t detector = -1;
2015 Float_t value = 0.0;
2017 for (Int_t k = 0; k < loop; k++) {
2018 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2020 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2021 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2022 Float_t min = 100.0;
2023 for (Int_t row = 0; row < rowMax; row++) {
2024 for (Int_t col = 0; col < colMax; col++) {
2025 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2026 mean += -TMath::Abs(value);
2030 if(count > 0) mean = mean/count;
2032 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2033 object->SetValue(detector,-TMath::Abs(value));
2039 //_____________________________________________________________________________
2040 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2043 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2044 // You need first to create the object for the detectors,
2045 // where the mean value is put.
2046 // This object has to be written in the database
2049 // Create the DetObject
2050 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2053 for(Int_t k = 0; k < 540; k++){
2054 AliTRDCalROC *calROC = object->GetCalROC(k);
2055 Int_t nchannels = calROC->GetNchannels();
2056 for(Int_t ch = 0; ch < nchannels; ch++){
2057 calROC->SetValue(ch,1.0);
2063 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2064 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2065 Int_t detector = -1;
2066 Float_t value = 0.0;
2068 for (Int_t k = 0; k < loop; k++) {
2069 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2070 AliTRDCalROC *calROC = object->GetCalROC(detector);
2071 Float_t mean = detobject->GetValue(detector);
2072 if(TMath::Abs(mean) <= 0.0000000001) continue;
2073 Int_t rowMax = calROC->GetNrows();
2074 Int_t colMax = calROC->GetNcols();
2075 for (Int_t row = 0; row < rowMax; row++) {
2076 for (Int_t col = 0; col < colMax; col++) {
2077 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2078 if(value > 0) value = value*scaleFitFactor;
2079 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2087 //_____________________________________________________________________________
2088 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2091 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2092 // You need first to create the object for the detectors,
2093 // where the mean value is put.
2094 // This object has to be written in the database
2097 // Create the DetObject
2098 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2101 for(Int_t k = 0; k < 540; k++){
2102 AliTRDCalROC *calROC = object->GetCalROC(k);
2103 Int_t nchannels = calROC->GetNchannels();
2104 for(Int_t ch = 0; ch < nchannels; ch++){
2105 calROC->SetValue(ch,1.0);
2111 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2112 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2113 Int_t detector = -1;
2114 Float_t value = 0.0;
2116 for (Int_t k = 0; k < loop; k++) {
2117 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2118 AliTRDCalROC *calROC = object->GetCalROC(detector);
2119 Float_t mean = detobject->GetValue(detector);
2120 if(mean == 0) continue;
2121 Int_t rowMax = calROC->GetNrows();
2122 Int_t colMax = calROC->GetNcols();
2123 for (Int_t row = 0; row < rowMax; row++) {
2124 for (Int_t col = 0; col < colMax; col++) {
2125 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2126 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2134 //_____________________________________________________________________________
2135 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2138 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2139 // You need first to create the object for the detectors,
2140 // where the mean value is put.
2141 // This object has to be written in the database
2144 // Create the DetObject
2145 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2148 for(Int_t k = 0; k < 540; k++){
2149 AliTRDCalROC *calROC = object->GetCalROC(k);
2150 Int_t nchannels = calROC->GetNchannels();
2151 for(Int_t ch = 0; ch < nchannels; ch++){
2152 calROC->SetValue(ch,0.0);
2158 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2159 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2160 Int_t detector = -1;
2161 Float_t value = 0.0;
2163 for (Int_t k = 0; k < loop; k++) {
2164 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2165 AliTRDCalROC *calROC = object->GetCalROC(detector);
2166 Float_t min = detobject->GetValue(detector);
2167 Int_t rowMax = calROC->GetNrows();
2168 Int_t colMax = calROC->GetNcols();
2169 for (Int_t row = 0; row < rowMax; row++) {
2170 for (Int_t col = 0; col < colMax; col++) {
2171 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2173 if(value > 70.0) value = value - 100.0;
2175 calROC->SetValue(col,row,value-min);
2183 //_____________________________________________________________________________
2184 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2187 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2188 // This object has to be written in the database
2191 // Create the DetObject
2192 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2194 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2195 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2196 Int_t detector = -1;
2197 Float_t value = 0.0;
2199 for (Int_t k = 0; k < loop; k++) {
2200 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2201 AliTRDCalROC *calROC = object->GetCalROC(detector);
2202 Int_t rowMax = calROC->GetNrows();
2203 Int_t colMax = calROC->GetNcols();
2204 for (Int_t row = 0; row < rowMax; row++) {
2205 for (Int_t col = 0; col < colMax; col++) {
2206 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2207 calROC->SetValue(col,row,TMath::Abs(value));
2215 //_____________________________________________________________________________
2216 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2219 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2220 // 0 successful fit 1 not successful fit
2221 // mean is the mean value over the successful fit
2222 // do not use it for t0: no meaning
2225 // Create the CalObject
2226 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2230 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2232 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2233 for(Int_t k = 0; k < 540; k++){
2234 object->SetValue(k,1.0);
2237 Int_t detector = -1;
2238 Float_t value = 0.0;
2240 for (Int_t k = 0; k < loop; k++) {
2241 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2242 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2243 if(value <= 0) object->SetValue(detector,1.0);
2245 object->SetValue(detector,0.0);
2250 if(count > 0) mean /= count;
2253 //_____________________________________________________________________________
2254 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2257 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2258 // 0 not successful fit 1 successful fit
2259 // mean mean value over the successful fit
2262 // Create the CalObject
2263 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2267 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2269 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2270 for(Int_t k = 0; k < 540; k++){
2271 AliTRDCalROC *calROC = object->GetCalROC(k);
2272 Int_t nchannels = calROC->GetNchannels();
2273 for(Int_t ch = 0; ch < nchannels; ch++){
2274 calROC->SetValue(ch,1.0);
2278 Int_t detector = -1;
2279 Float_t value = 0.0;
2281 for (Int_t k = 0; k < loop; k++) {
2282 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2283 AliTRDCalROC *calROC = object->GetCalROC(detector);
2284 Int_t nchannels = calROC->GetNchannels();
2285 for (Int_t ch = 0; ch < nchannels; ch++) {
2286 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2287 if(value <= 0) calROC->SetValue(ch,1.0);
2289 calROC->SetValue(ch,0.0);
2295 if(count > 0) mean /= count;
2298 //_____________________________________________________________________________
2299 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2302 // Set FitPH if 1 then each detector will be fitted
2305 if (periodeFitPH > 0) {
2306 fFitPHPeriode = periodeFitPH;
2309 AliInfo("periodeFitPH must be higher than 0!");
2313 //_____________________________________________________________________________
2314 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2317 // The fit of the deposited charge distribution begins at
2318 // histo->Mean()/beginFitCharge
2319 // You can here set beginFitCharge
2322 if (beginFitCharge > 0) {
2323 fBeginFitCharge = beginFitCharge;
2326 AliInfo("beginFitCharge must be strict positif!");
2331 //_____________________________________________________________________________
2332 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2335 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2336 // You can here set t0Shift0
2340 fT0Shift0 = t0Shift;
2343 AliInfo("t0Shift0 must be strict positif!");
2348 //_____________________________________________________________________________
2349 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2352 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2353 // You can here set t0Shift1
2357 fT0Shift1 = t0Shift;
2360 AliInfo("t0Shift must be strict positif!");
2365 //_____________________________________________________________________________
2366 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2369 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2370 // You can here set rangeFitPRF
2373 if ((rangeFitPRF > 0) &&
2374 (rangeFitPRF <= 1.5)) {
2375 fRangeFitPRF = rangeFitPRF;
2378 AliInfo("rangeFitPRF must be between 0 and 1.0");
2383 //_____________________________________________________________________________
2384 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2387 // Minimum entries for fitting
2390 if (minEntries > 0) {
2391 fMinEntries = minEntries;
2394 AliInfo("fMinEntries must be >= 0.");
2399 //_____________________________________________________________________________
2400 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2403 // Rebin with rebin time less bins the Ch histo
2404 // You can set here rebin that should divide the number of bins of CH histo
2409 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2412 AliInfo("You have to choose a positiv value!");
2416 //_____________________________________________________________________________
2417 Bool_t AliTRDCalibraFit::FillVectorFit()
2420 // For the Fit functions fill the vector Fit
2423 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2426 if (GetStack(fCountDet) == 2) {
2433 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2434 Float_t *coef = new Float_t[ntotal];
2435 for (Int_t i = 0; i < ntotal; i++) {
2436 coef[i] = fCurrentCoefDetector[i];
2439 Int_t detector = fCountDet;
2441 fitInfo->SetCoef(coef);
2442 fitInfo->SetDetector(detector);
2443 fVectorFit.Add((TObject *) fitInfo);
2448 //_____________________________________________________________________________
2449 Bool_t AliTRDCalibraFit::FillVectorFit2()
2452 // For the Fit functions fill the vector Fit
2455 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2458 if (GetStack(fCountDet) == 2) {
2465 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2466 Float_t *coef = new Float_t[ntotal];
2467 for (Int_t i = 0; i < ntotal; i++) {
2468 coef[i] = fCurrentCoefDetector2[i];
2471 Int_t detector = fCountDet;
2473 fitInfo->SetCoef(coef);
2474 fitInfo->SetDetector(detector);
2475 fVectorFit2.Add((TObject *) fitInfo);
2480 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2481 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2484 // Init the number of expected bins and fDect1[i] fDect2[i]
2487 gStyle->SetPalette(1);
2488 gStyle->SetOptStat(1111);
2489 gStyle->SetPadBorderMode(0);
2490 gStyle->SetCanvasColor(10);
2491 gStyle->SetPadLeftMargin(0.13);
2492 gStyle->SetPadRightMargin(0.01);
2494 // Mode groups of pads: the total number of bins!
2495 CalculNumberOfBinsExpected(i);
2497 // Quick verification that we have the good pad calibration mode!
2498 if (fNumberOfBinsExpected != nbins) {
2499 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2503 // Security for fDebug 3 and 4
2504 if ((fDebugLevel >= 3) &&
2508 AliInfo("This detector doesn't exit!");
2512 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2513 CalculDect1Dect2(i);
2518 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2519 Bool_t AliTRDCalibraFit::InitFitCH()
2522 // Init the fVectorFitCH for normalisation
2523 // Init the histo for debugging
2528 fScaleFitFactor = 0.0;
2529 fCurrentCoefDetector = new Float_t[2304];
2530 for (Int_t k = 0; k < 2304; k++) {
2531 fCurrentCoefDetector[k] = 0.0;
2533 fVectorFit.SetName("gainfactorscoefficients");
2535 // fDebug == 0 nothing
2536 // fDebug == 1 and fFitVoir no histo
2537 if (fDebugLevel == 1) {
2538 if(!CheckFitVoir()) return kFALSE;
2540 //Get the CalDet object
2542 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2544 AliInfo("Could not get calibDB");
2547 if(fCalDet) delete fCalDet;
2548 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2551 Float_t devalue = 1.0;
2552 if(fCalDet) delete fCalDet;
2553 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2554 for(Int_t k = 0; k < 540; k++){
2555 fCalDet->SetValue(k,devalue);
2561 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2562 Bool_t AliTRDCalibraFit::InitFitPH()
2565 // Init the arrays of results
2566 // Init the histos for debugging
2570 fVectorFit.SetName("driftvelocitycoefficients");
2571 fVectorFit2.SetName("t0coefficients");
2573 fCurrentCoefDetector = new Float_t[2304];
2574 for (Int_t k = 0; k < 2304; k++) {
2575 fCurrentCoefDetector[k] = 0.0;
2578 fCurrentCoefDetector2 = new Float_t[2304];
2579 for (Int_t k = 0; k < 2304; k++) {
2580 fCurrentCoefDetector2[k] = 0.0;
2583 //fDebug == 0 nothing
2584 // fDebug == 1 and fFitVoir no histo
2585 if (fDebugLevel == 1) {
2586 if(!CheckFitVoir()) return kFALSE;
2588 //Get the CalDet object
2590 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2592 AliInfo("Could not get calibDB");
2595 if(fCalDet) delete fCalDet;
2596 if(fCalDet2) delete fCalDet2;
2597 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2598 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2601 Float_t devalue = 1.5;
2602 Float_t devalue2 = 0.0;
2603 if(fCalDet) delete fCalDet;
2604 if(fCalDet2) delete fCalDet2;
2605 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2606 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2607 for(Int_t k = 0; k < 540; k++){
2608 fCalDet->SetValue(k,devalue);
2609 fCalDet2->SetValue(k,devalue2);
2614 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2615 Bool_t AliTRDCalibraFit::InitFitPRF()
2618 // Init the calibration mode (Nz, Nrphi), the histograms for
2619 // debugging the fit methods if fDebug > 0,
2623 fVectorFit.SetName("prfwidthcoefficients");
2625 fCurrentCoefDetector = new Float_t[2304];
2626 for (Int_t k = 0; k < 2304; k++) {
2627 fCurrentCoefDetector[k] = 0.0;
2630 // fDebug == 0 nothing
2631 // fDebug == 1 and fFitVoir no histo
2632 if (fDebugLevel == 1) {
2633 if(!CheckFitVoir()) return kFALSE;
2637 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2638 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2641 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2646 fCurrentCoefDetector = new Float_t[2304];
2647 fCurrentCoefDetector2 = new Float_t[2304];
2648 for (Int_t k = 0; k < 2304; k++) {
2649 fCurrentCoefDetector[k] = 0.0;
2650 fCurrentCoefDetector2[k] = 0.0;
2653 //printf("test0\n");
2655 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2657 AliInfo("Could not get calibDB");
2661 //Get the CalDet object
2663 if(fCalDet) delete fCalDet;
2664 if(fCalDet2) delete fCalDet2;
2665 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2666 //printf("test1\n");
2667 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2668 //printf("test2\n");
2669 for(Int_t k = 0; k < 540; k++){
2670 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2672 //printf("test3\n");
2675 Float_t devalue = 1.5;
2676 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2677 if(fCalDet) delete fCalDet;
2678 if(fCalDet2) delete fCalDet2;
2679 //printf("test1\n");
2680 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2681 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2682 //printf("test2\n");
2683 for(Int_t k = 0; k < 540; k++){
2684 fCalDet->SetValue(k,devalue);
2685 fCalDet2->SetValue(k,devalue2);
2687 //printf("test3\n");
2692 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2693 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2696 // Init the current detector where we are fCountDet and the
2697 // next fCount for the functions Fit...
2700 // Loop on the Xbins of ch!!
2701 fCountDet = -1; // Current detector
2702 fCount = 0; // To find the next detector
2705 if (fDebugLevel >= 3) {
2706 // Set countdet to the detector
2707 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2708 // Set counter to write at the end of the detector
2710 // Get the right calib objects
2713 if(fDebugLevel == 1) {
2715 fCalibraMode->CalculXBins(fCountDet,i);
2716 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2717 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2719 fCalibraMode->CalculXBins(fCountDet,i);
2720 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2726 fCount = fCalibraMode->GetXbins(i);
2728 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2729 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2730 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2731 ,(Int_t) GetStack(fCountDet)
2732 ,(Int_t) GetSector(fCountDet),i);
2735 //_______________________________________________________________________________
2736 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2739 // Calculate the number of bins expected (calibration groups)
2742 fNumberOfBinsExpected = 0;
2744 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2745 fNumberOfBinsExpected = 1;
2749 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2750 fNumberOfBinsExpected = 18;
2754 fCalibraMode->ModePadCalibration(2,i);
2755 fCalibraMode->ModePadFragmentation(0,2,0,i);
2756 fCalibraMode->SetDetChamb2(i);
2757 if (fDebugLevel > 1) {
2758 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2760 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2761 fCalibraMode->ModePadCalibration(0,i);
2762 fCalibraMode->ModePadFragmentation(0,0,0,i);
2763 fCalibraMode->SetDetChamb0(i);
2764 if (fDebugLevel > 1) {
2765 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2767 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2770 //_______________________________________________________________________________
2771 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2774 // Calculate the range of fits
2779 if (fDebugLevel == 1) {
2783 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2785 fDect2 = fNumberOfBinsExpected;
2787 if (fDebugLevel >= 3) {
2788 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2789 fCalibraMode->CalculXBins(fCountDet,i);
2790 fDect1 = fCalibraMode->GetXbins(i);
2791 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2792 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2793 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2794 ,(Int_t) GetStack(fCountDet)
2795 ,(Int_t) GetSector(fCountDet),i);
2796 // Set for the next detector
2797 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2800 //_______________________________________________________________________________
2801 Bool_t AliTRDCalibraFit::CheckFitVoir()
2804 // Check if fFitVoir is in the range
2807 if (fFitVoir < fNumberOfBinsExpected) {
2808 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2811 AliInfo("fFitVoir is out of range of the histo!");
2816 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2817 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2820 // See if we are in a new detector and update the
2821 // variables fNfragZ and fNfragRphi if yes
2822 // Will never happen for only one detector (3 and 4)
2823 // Doesn't matter for 2
2825 if (fCount == idect) {
2826 // On en est au detector (or first detector in the group)
2828 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2829 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2830 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2831 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2832 ,(Int_t) GetStack(fCountDet)
2833 ,(Int_t) GetSector(fCountDet),i);
2834 // Set for the next detector
2835 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2840 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2841 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2844 // Reconstruct the min pad row, max pad row, min pad col and
2845 // max pad col of the calibration group for the Fit functions
2846 // idect is the calibration group inside the detector
2848 if (fDebugLevel != 1) {
2849 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2851 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2852 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2854 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2855 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2858 // For the case where there are not enough entries in the histograms
2859 // of the calibration group, the value present in the choosen database
2860 // will be put. A negativ sign enables to know that a fit was not possible.
2863 if (fDebugLevel == 1) {
2864 AliInfo("The element has not enough statistic to be fitted");
2866 else if (fNbDet > 0){
2867 Int_t firstdetector = fCountDet;
2868 Int_t lastdetector = fCountDet+fNbDet;
2869 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2870 ,idect,firstdetector,lastdetector));
2871 // loop over detectors
2872 for(Int_t det = firstdetector; det < lastdetector; det++){
2874 //Set the calibration object again
2878 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2880 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2881 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2882 ,(Int_t) GetStack(fCountDet)
2883 ,(Int_t) GetSector(fCountDet),0);
2884 // Reconstruct row min row max
2885 ReconstructFitRowMinRowMax(idect,0);
2887 // Calcul the coef from the database choosen for the detector
2888 CalculChargeCoefMean(kFALSE);
2890 //stack 2, not stack 2
2892 if(GetStack(fCountDet) == 2) factor = 12;
2895 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2896 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2897 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2898 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2902 //Put default value negative
2903 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2904 fCurrentCoefE = 0.0;
2909 if(fDebugLevel > 1){
2911 if ( !fDebugStreamer ) {
2913 TDirectory *backup = gDirectory;
2914 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2915 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2918 Int_t detector = fCountDet;
2919 Int_t caligroup = idect;
2920 Short_t rowmin = fCalibraMode->GetRowMin(0);
2921 Short_t rowmax = fCalibraMode->GetRowMax(0);
2922 Short_t colmin = fCalibraMode->GetColMin(0);
2923 Short_t colmax = fCalibraMode->GetColMax(0);
2924 Float_t gf = fCurrentCoef[0];
2925 Float_t gfs = fCurrentCoef[1];
2926 Float_t gfE = fCurrentCoefE;
2928 (*fDebugStreamer) << "FillFillCH" <<
2929 "detector=" << detector <<
2930 "caligroup=" << caligroup <<
2931 "rowmin=" << rowmin <<
2932 "rowmax=" << rowmax <<
2933 "colmin=" << colmin <<
2934 "colmax=" << colmax <<
2942 for (Int_t k = 0; k < 2304; k++) {
2943 fCurrentCoefDetector[k] = 0.0;
2947 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2951 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2952 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2954 // Calcul the coef from the database choosen
2955 CalculChargeCoefMean(kFALSE);
2957 //stack 2, not stack 2
2959 if(GetStack(fCountDet) == 2) factor = 12;
2962 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2963 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2964 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2965 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2969 //Put default value negative
2970 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2971 fCurrentCoefE = 0.0;
2980 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2981 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2984 // For the case where there are not enough entries in the histograms
2985 // of the calibration group, the value present in the choosen database
2986 // will be put. A negativ sign enables to know that a fit was not possible.
2988 if (fDebugLevel == 1) {
2989 AliInfo("The element has not enough statistic to be fitted");
2991 else if (fNbDet > 0) {
2993 Int_t firstdetector = fCountDet;
2994 Int_t lastdetector = fCountDet+fNbDet;
2995 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2996 ,idect,firstdetector,lastdetector));
2997 // loop over detectors
2998 for(Int_t det = firstdetector; det < lastdetector; det++){
3000 //Set the calibration object again
3004 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3006 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3007 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3008 ,(Int_t) GetStack(fCountDet)
3009 ,(Int_t) GetSector(fCountDet),1);
3010 // Reconstruct row min row max
3011 ReconstructFitRowMinRowMax(idect,1);
3013 // Calcul the coef from the database choosen for the detector
3014 CalculVdriftCoefMean();
3017 //stack 2, not stack 2
3019 if(GetStack(fCountDet) == 2) factor = 12;
3022 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3023 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3024 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3025 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3026 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3030 //Put default value negative
3031 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3032 fCurrentCoefE = 0.0;
3033 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3034 fCurrentCoefE2 = 0.0;
3040 if(fDebugLevel > 1){
3042 if ( !fDebugStreamer ) {
3044 TDirectory *backup = gDirectory;
3045 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3046 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3050 Int_t detector = fCountDet;
3051 Int_t caligroup = idect;
3052 Short_t rowmin = fCalibraMode->GetRowMin(1);
3053 Short_t rowmax = fCalibraMode->GetRowMax(1);
3054 Short_t colmin = fCalibraMode->GetColMin(1);
3055 Short_t colmax = fCalibraMode->GetColMax(1);
3056 Float_t vf = fCurrentCoef[0];
3057 Float_t vs = fCurrentCoef[1];
3058 Float_t vfE = fCurrentCoefE;
3059 Float_t t0f = fCurrentCoef2[0];
3060 Float_t t0s = fCurrentCoef2[1];
3061 Float_t t0E = fCurrentCoefE2;
3065 (* fDebugStreamer) << "FillFillPH"<<
3066 "detector="<<detector<<
3067 "nentries="<<nentries<<
3068 "caligroup="<<caligroup<<
3082 for (Int_t k = 0; k < 2304; k++) {
3083 fCurrentCoefDetector[k] = 0.0;
3084 fCurrentCoefDetector2[k] = 0.0;
3088 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3092 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3093 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3095 CalculVdriftCoefMean();
3098 //stack 2 and not stack 2
3100 if(GetStack(fCountDet) == 2) factor = 12;
3104 // Fill the fCurrentCoefDetector 2
3105 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3106 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3107 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3108 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3112 // Put the default value
3113 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3114 fCurrentCoefE = 0.0;
3115 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3116 fCurrentCoefE2 = 0.0;
3118 FillFillPH(idect,nentries);
3127 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3128 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3131 // For the case where there are not enough entries in the histograms
3132 // of the calibration group, the value present in the choosen database
3133 // will be put. A negativ sign enables to know that a fit was not possible.
3136 if (fDebugLevel == 1) {
3137 AliInfo("The element has not enough statistic to be fitted");
3139 else if (fNbDet > 0){
3141 Int_t firstdetector = fCountDet;
3142 Int_t lastdetector = fCountDet+fNbDet;
3143 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3144 ,idect,firstdetector,lastdetector));
3146 // loop over detectors
3147 for(Int_t det = firstdetector; det < lastdetector; det++){
3149 //Set the calibration object again
3153 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3155 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3156 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3157 ,(Int_t) GetStack(fCountDet)
3158 ,(Int_t) GetSector(fCountDet),2);
3159 // Reconstruct row min row max
3160 ReconstructFitRowMinRowMax(idect,2);
3162 // Calcul the coef from the database choosen for the detector
3163 CalculPRFCoefMean();
3165 //stack 2, not stack 2
3167 if(GetStack(fCountDet) == 2) factor = 12;
3170 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3171 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3172 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3173 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3177 //Put default value negative
3178 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3179 fCurrentCoefE = 0.0;
3184 if(fDebugLevel > 1){
3186 if ( !fDebugStreamer ) {
3188 TDirectory *backup = gDirectory;
3189 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3190 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3193 Int_t detector = fCountDet;
3194 Int_t layer = GetLayer(fCountDet);
3195 Int_t caligroup = idect;
3196 Short_t rowmin = fCalibraMode->GetRowMin(2);
3197 Short_t rowmax = fCalibraMode->GetRowMax(2);
3198 Short_t colmin = fCalibraMode->GetColMin(2);
3199 Short_t colmax = fCalibraMode->GetColMax(2);
3200 Float_t widf = fCurrentCoef[0];
3201 Float_t wids = fCurrentCoef[1];
3202 Float_t widfE = fCurrentCoefE;
3204 (* fDebugStreamer) << "FillFillPRF"<<
3205 "detector="<<detector<<
3207 "caligroup="<<caligroup<<
3218 for (Int_t k = 0; k < 2304; k++) {
3219 fCurrentCoefDetector[k] = 0.0;
3223 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3227 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3228 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3230 CalculPRFCoefMean();
3232 // stack 2 and not stack 2
3234 if(GetStack(fCountDet) == 2) factor = 12;
3238 // Fill the fCurrentCoefDetector
3239 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3240 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3241 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3245 // Put the default value
3246 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3247 fCurrentCoefE = 0.0;
3255 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3256 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3259 // For the case where there are not enough entries in the histograms
3260 // of the calibration group, the value present in the choosen database
3261 // will be put. A negativ sign enables to know that a fit was not possible.
3264 // Calcul the coef from the database choosen
3265 CalculVdriftLorentzCoef();
3268 if(GetStack(fCountDet) == 2) factor = 1728;
3272 // Fill the fCurrentCoefDetector
3273 for (Int_t k = 0; k < factor; k++) {
3274 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3275 // should be negative
3276 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3280 //Put default opposite sign
3281 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3282 fCurrentCoefE = 0.0;
3283 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3284 fCurrentCoefE2 = 0.0;
3286 FillFillLinearFitter();
3291 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3292 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3295 // Fill the coefficients found with the fits or other
3296 // methods from the Fit functions
3299 if (fDebugLevel != 1) {
3301 Int_t firstdetector = fCountDet;
3302 Int_t lastdetector = fCountDet+fNbDet;
3303 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3304 ,idect,firstdetector,lastdetector));
3305 // loop over detectors
3306 for(Int_t det = firstdetector; det < lastdetector; det++){
3308 //Set the calibration object again
3312 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3314 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3315 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3316 ,(Int_t) GetStack(fCountDet)
3317 ,(Int_t) GetSector(fCountDet),0);
3318 // Reconstruct row min row max
3319 ReconstructFitRowMinRowMax(idect,0);
3321 // Calcul the coef from the database choosen for the detector
3322 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3323 else CalculChargeCoefMean(kTRUE);
3325 //stack 2, not stack 2
3327 if(GetStack(fCountDet) == 2) factor = 12;
3330 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3331 Double_t coeftoput = 1.0;
3332 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3333 else coeftoput = fCurrentCoef[0];
3334 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3335 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3336 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3343 if(fDebugLevel > 1){
3345 if ( !fDebugStreamer ) {
3347 TDirectory *backup = gDirectory;
3348 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3349 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3352 Int_t detector = fCountDet;
3353 Int_t caligroup = idect;
3354 Short_t rowmin = fCalibraMode->GetRowMin(0);
3355 Short_t rowmax = fCalibraMode->GetRowMax(0);
3356 Short_t colmin = fCalibraMode->GetColMin(0);
3357 Short_t colmax = fCalibraMode->GetColMax(0);
3358 Float_t gf = fCurrentCoef[0];
3359 Float_t gfs = fCurrentCoef[1];
3360 Float_t gfE = fCurrentCoefE;
3362 (*fDebugStreamer) << "FillFillCH" <<
3363 "detector=" << detector <<
3364 "caligroup=" << caligroup <<
3365 "rowmin=" << rowmin <<
3366 "rowmax=" << rowmax <<
3367 "colmin=" << colmin <<
3368 "colmax=" << colmax <<
3376 for (Int_t k = 0; k < 2304; k++) {
3377 fCurrentCoefDetector[k] = 0.0;
3381 //printf("Check the count now: fCountDet %d\n",fCountDet);
3386 if(GetStack(fCountDet) == 2) factor = 12;
3389 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3390 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3391 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3402 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3403 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3406 // Fill the coefficients found with the fits or other
3407 // methods from the Fit functions
3410 if (fDebugLevel != 1) {
3413 Int_t firstdetector = fCountDet;
3414 Int_t lastdetector = fCountDet+fNbDet;
3415 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3416 ,idect,firstdetector,lastdetector));
3418 // loop over detectors
3419 for(Int_t det = firstdetector; det < lastdetector; det++){
3421 //Set the calibration object again
3425 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3427 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3428 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3429 ,(Int_t) GetStack(fCountDet)
3430 ,(Int_t) GetSector(fCountDet),1);
3431 // Reconstruct row min row max
3432 ReconstructFitRowMinRowMax(idect,1);
3434 // Calcul the coef from the database choosen for the detector
3435 CalculVdriftCoefMean();
3438 //stack 2, not stack 2
3440 if(GetStack(fCountDet) == 2) factor = 12;
3443 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3444 Double_t coeftoput = 1.5;
3445 Double_t coeftoput2 = 0.0;
3447 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3448 else coeftoput = fCurrentCoef[0];
3450 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3451 else coeftoput2 = fCurrentCoef2[0];
3453 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3454 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3455 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3456 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3464 if(fDebugLevel > 1){
3466 if ( !fDebugStreamer ) {
3468 TDirectory *backup = gDirectory;
3469 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3470 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3474 Int_t detector = fCountDet;
3475 Int_t caligroup = idect;
3476 Short_t rowmin = fCalibraMode->GetRowMin(1);
3477 Short_t rowmax = fCalibraMode->GetRowMax(1);
3478 Short_t colmin = fCalibraMode->GetColMin(1);
3479 Short_t colmax = fCalibraMode->GetColMax(1);
3480 Float_t vf = fCurrentCoef[0];
3481 Float_t vs = fCurrentCoef[1];
3482 Float_t vfE = fCurrentCoefE;
3483 Float_t t0f = fCurrentCoef2[0];
3484 Float_t t0s = fCurrentCoef2[1];
3485 Float_t t0E = fCurrentCoefE2;
3489 (* fDebugStreamer) << "FillFillPH"<<
3490 "detector="<<detector<<
3491 "nentries="<<nentries<<
3492 "caligroup="<<caligroup<<
3506 for (Int_t k = 0; k < 2304; k++) {
3507 fCurrentCoefDetector[k] = 0.0;
3508 fCurrentCoefDetector2[k] = 0.0;
3512 //printf("Check the count now: fCountDet %d\n",fCountDet);
3517 if(GetStack(fCountDet) == 2) factor = 12;
3520 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3521 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3522 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3523 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3527 FillFillPH(idect,nentries);
3532 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3533 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3536 // Fill the coefficients found with the fits or other
3537 // methods from the Fit functions
3540 if (fDebugLevel != 1) {
3543 Int_t firstdetector = fCountDet;
3544 Int_t lastdetector = fCountDet+fNbDet;
3545 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3546 ,idect,firstdetector,lastdetector));
3548 // loop over detectors
3549 for(Int_t det = firstdetector; det < lastdetector; det++){
3551 //Set the calibration object again
3555 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3557 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3558 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3559 ,(Int_t) GetStack(fCountDet)
3560 ,(Int_t) GetSector(fCountDet),2);
3561 // Reconstruct row min row max
3562 ReconstructFitRowMinRowMax(idect,2);
3564 // Calcul the coef from the database choosen for the detector
3565 CalculPRFCoefMean();
3567 //stack 2, not stack 2
3569 if(GetStack(fCountDet) == 2) factor = 12;
3572 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3573 Double_t coeftoput = 1.0;
3574 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3575 else coeftoput = fCurrentCoef[0];
3576 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3577 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3578 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3585 if(fDebugLevel > 1){
3587 if ( !fDebugStreamer ) {
3589 TDirectory *backup = gDirectory;
3590 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3591 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3594 Int_t detector = fCountDet;
3595 Int_t layer = GetLayer(fCountDet);
3596 Int_t caligroup = idect;
3597 Short_t rowmin = fCalibraMode->GetRowMin(2);
3598 Short_t rowmax = fCalibraMode->GetRowMax(2);
3599 Short_t colmin = fCalibraMode->GetColMin(2);
3600 Short_t colmax = fCalibraMode->GetColMax(2);
3601 Float_t widf = fCurrentCoef[0];
3602 Float_t wids = fCurrentCoef[1];
3603 Float_t widfE = fCurrentCoefE;
3605 (* fDebugStreamer) << "FillFillPRF"<<
3606 "detector="<<detector<<
3608 "caligroup="<<caligroup<<
3619 for (Int_t k = 0; k < 2304; k++) {
3620 fCurrentCoefDetector[k] = 0.0;
3624 //printf("Check the count now: fCountDet %d\n",fCountDet);
3629 if(GetStack(fCountDet) == 2) factor = 12;
3632 // Pointer to the branch
3633 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3634 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3635 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3645 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3646 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3649 // Fill the coefficients found with the fits or other
3650 // methods from the Fit functions
3654 if(GetStack(fCountDet) == 2) factor = 1728;
3657 // Pointer to the branch
3658 for (Int_t k = 0; k < factor; k++) {
3659 fCurrentCoefDetector[k] = fCurrentCoef[0];
3660 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3663 FillFillLinearFitter();
3668 //________________________________________________________________________________
3669 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3672 // DebugStream and fVectorFit
3675 // End of one detector
3676 if ((idect == (fCount-1))) {
3679 for (Int_t k = 0; k < 2304; k++) {
3680 fCurrentCoefDetector[k] = 0.0;
3684 if(fDebugLevel > 1){
3686 if ( !fDebugStreamer ) {
3688 TDirectory *backup = gDirectory;
3689 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3690 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3693 Int_t detector = fCountDet;
3694 Int_t caligroup = idect;
3695 Short_t rowmin = fCalibraMode->GetRowMin(0);
3696 Short_t rowmax = fCalibraMode->GetRowMax(0);
3697 Short_t colmin = fCalibraMode->GetColMin(0);
3698 Short_t colmax = fCalibraMode->GetColMax(0);
3699 Float_t gf = fCurrentCoef[0];
3700 Float_t gfs = fCurrentCoef[1];
3701 Float_t gfE = fCurrentCoefE;
3703 (*fDebugStreamer) << "FillFillCH" <<
3704 "detector=" << detector <<
3705 "caligroup=" << caligroup <<
3706 "rowmin=" << rowmin <<
3707 "rowmax=" << rowmax <<
3708 "colmin=" << colmin <<
3709 "colmax=" << colmax <<
3717 //________________________________________________________________________________
3718 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3721 // DebugStream and fVectorFit and fVectorFit2
3724 // End of one detector
3725 if ((idect == (fCount-1))) {
3729 for (Int_t k = 0; k < 2304; k++) {
3730 fCurrentCoefDetector[k] = 0.0;
3731 fCurrentCoefDetector2[k] = 0.0;
3735 if(fDebugLevel > 1){
3737 if ( !fDebugStreamer ) {
3739 TDirectory *backup = gDirectory;
3740 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3741 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3745 Int_t detector = fCountDet;
3746 Int_t caligroup = idect;
3747 Short_t rowmin = fCalibraMode->GetRowMin(1);
3748 Short_t rowmax = fCalibraMode->GetRowMax(1);
3749 Short_t colmin = fCalibraMode->GetColMin(1);
3750 Short_t colmax = fCalibraMode->GetColMax(1);
3751 Float_t vf = fCurrentCoef[0];
3752 Float_t vs = fCurrentCoef[1];
3753 Float_t vfE = fCurrentCoefE;
3754 Float_t t0f = fCurrentCoef2[0];
3755 Float_t t0s = fCurrentCoef2[1];
3756 Float_t t0E = fCurrentCoefE2;
3760 (* fDebugStreamer) << "FillFillPH"<<
3761 "detector="<<detector<<
3762 "nentries="<<nentries<<
3763 "caligroup="<<caligroup<<
3778 //________________________________________________________________________________
3779 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3782 // DebugStream and fVectorFit
3785 // End of one detector
3786 if ((idect == (fCount-1))) {
3789 for (Int_t k = 0; k < 2304; k++) {
3790 fCurrentCoefDetector[k] = 0.0;
3795 if(fDebugLevel > 1){
3797 if ( !fDebugStreamer ) {
3799 TDirectory *backup = gDirectory;
3800 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3801 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3804 Int_t detector = fCountDet;
3805 Int_t layer = GetLayer(fCountDet);
3806 Int_t caligroup = idect;
3807 Short_t rowmin = fCalibraMode->GetRowMin(2);
3808 Short_t rowmax = fCalibraMode->GetRowMax(2);
3809 Short_t colmin = fCalibraMode->GetColMin(2);
3810 Short_t colmax = fCalibraMode->GetColMax(2);
3811 Float_t widf = fCurrentCoef[0];
3812 Float_t wids = fCurrentCoef[1];
3813 Float_t widfE = fCurrentCoefE;
3815 (* fDebugStreamer) << "FillFillPRF"<<
3816 "detector="<<detector<<
3818 "caligroup="<<caligroup<<
3830 //________________________________________________________________________________
3831 void AliTRDCalibraFit::FillFillLinearFitter()
3834 // DebugStream and fVectorFit
3837 // End of one detector
3843 for (Int_t k = 0; k < 2304; k++) {
3844 fCurrentCoefDetector[k] = 0.0;
3845 fCurrentCoefDetector2[k] = 0.0;
3849 if(fDebugLevel > 1){
3851 if ( !fDebugStreamer ) {
3853 TDirectory *backup = gDirectory;
3854 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3855 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3858 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3859 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3860 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3861 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3862 Float_t tiltangle = padplane->GetTiltingAngle();
3863 Int_t detector = fCountDet;
3864 Int_t stack = GetStack(fCountDet);
3865 Int_t layer = GetLayer(fCountDet);
3866 Float_t vf = fCurrentCoef[0];
3867 Float_t vs = fCurrentCoef[1];
3868 Float_t vfE = fCurrentCoefE;
3869 Float_t lorentzangler = fCurrentCoef2[0];
3870 Float_t elorentzangler = fCurrentCoefE2;
3871 Float_t lorentzangles = fCurrentCoef2[1];
3873 (* fDebugStreamer) << "FillFillLinearFitter"<<
3874 "detector="<<detector<<
3879 "tiltangle="<<tiltangle<<
3883 "lorentzangler="<<lorentzangler<<
3884 "Elorentzangler="<<elorentzangler<<
3885 "lorentzangles="<<lorentzangles<<
3891 //____________Calcul Coef Mean_________________________________________________
3893 //_____________________________________________________________________________
3894 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3897 // For the detector Dect calcul the mean time 0
3898 // for the calibration group idect from the choosen database
3901 fCurrentCoef2[1] = 0.0;
3902 if(fDebugLevel != 1){
3903 if(((fCalibraMode->GetNz(1) > 0) ||
3904 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3906 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3907 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3908 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3912 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3918 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3922 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3923 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3924 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3927 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3935 //_____________________________________________________________________________
3936 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3939 // For the detector Dect calcul the mean gain factor
3940 // for the calibration group idect from the choosen database
3943 fCurrentCoef[1] = 0.0;
3944 if(fDebugLevel != 1){
3945 if (((fCalibraMode->GetNz(0) > 0) ||
3946 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3947 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3948 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3949 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3950 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3953 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3957 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3958 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3963 //_____________________________________________________________________________
3964 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3967 // For the detector Dect calcul the mean sigma of pad response
3968 // function for the calibration group idect from the choosen database
3971 fCurrentCoef[1] = 0.0;
3972 if(fDebugLevel != 1){
3973 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3974 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3975 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3978 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3982 //_____________________________________________________________________________
3983 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3986 // For the detector dect calcul the mean drift velocity for the
3987 // calibration group idect from the choosen database
3990 fCurrentCoef[1] = 0.0;
3991 if(fDebugLevel != 1){
3992 if (((fCalibraMode->GetNz(1) > 0) ||
3993 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3995 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3996 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3997 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4001 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4006 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4011 //_____________________________________________________________________________
4012 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4015 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4018 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
4019 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4023 //_____________________________________________________________________________
4024 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4027 // Default width of the PRF if there is no database as reference
4032 //case 0: return 0.515;
4033 //case 1: return 0.502;
4034 //case 2: return 0.491;
4035 //case 3: return 0.481;
4036 //case 4: return 0.471;
4037 //case 5: return 0.463;
4038 //default: return 0.0;
4041 case 0: return 0.538429;
4042 case 1: return 0.524302;
4043 case 2: return 0.511591;
4044 case 3: return 0.500140;
4045 case 4: return 0.489821;
4046 case 5: return 0.480524;
4047 default: return 0.0;
4050 //________________________________________________________________________________
4051 void AliTRDCalibraFit::SetCalROC(Int_t i)
4054 // Set the calib object for fCountDet
4057 Float_t value = 0.0;
4059 //Get the CalDet object
4061 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4063 AliInfo("Could not get calibDB");
4070 fCalROC->~AliTRDCalROC();
4071 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4072 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4076 fCalROC->~AliTRDCalROC();
4077 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4078 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4080 fCalROC2->~AliTRDCalROC();
4081 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4082 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4086 fCalROC->~AliTRDCalROC();
4087 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4088 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4097 if(fCalROC) delete fCalROC;
4098 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4099 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4100 fCalROC->SetValue(k,1.0);
4104 if(fCalROC) delete fCalROC;
4105 if(fCalROC2) delete fCalROC2;
4106 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4107 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4108 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4109 fCalROC->SetValue(k,1.0);
4110 fCalROC2->SetValue(k,0.0);
4114 if(fCalROC) delete fCalROC;
4115 value = GetPRFDefault(GetLayer(fCountDet));
4116 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4117 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4118 fCalROC->SetValue(k,value);
4126 //____________Fit Methods______________________________________________________
4128 //_____________________________________________________________________________
4129 void AliTRDCalibraFit::FitPente(TH1* projPH)
4132 // Slope methode for the drift velocity
4136 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4143 fCurrentCoefE = 0.0;
4144 fCurrentCoefE2 = 0.0;
4145 fCurrentCoef[0] = 0.0;
4146 fCurrentCoef2[0] = 0.0;
4147 TLine *line = new TLine();
4150 TAxis *xpph = projPH->GetXaxis();
4151 Int_t nbins = xpph->GetNbins();
4152 Double_t lowedge = xpph->GetBinLowEdge(1);
4153 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4154 Double_t widbins = (upedge - lowedge) / nbins;
4155 Double_t limit = upedge + 0.5 * widbins;
4158 // Beginning of the signal
4159 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4160 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4161 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4163 binmax = (Int_t) pentea->GetMaximumBin();
4166 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4168 if (binmax >= nbins) {
4171 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4173 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4174 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4175 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4176 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4177 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4178 if (TMath::Abs(l3P2am) > 0.00000001) {
4179 fPhd[0] = -(l3P1am / (2 * l3P2am));
4182 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4183 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4186 // Amplification region
4189 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4190 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4197 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4199 if (binmax >= nbins) {
4202 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4204 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4205 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4206 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4207 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4208 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4209 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4210 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4212 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4213 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4216 fCurrentCoefE2 = fCurrentCoefE;
4219 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4220 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4221 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4224 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4227 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4229 if (binmin >= nbins) {
4232 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4237 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4238 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4239 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4240 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4241 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4242 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4243 if (TMath::Abs(l3P2dr) > 0.00000001) {
4244 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4246 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4247 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4249 Float_t fPhdt0 = 0.0;
4250 Float_t t0Shift = 0.0;
4253 t0Shift = fT0Shift1;
4257 t0Shift = fT0Shift0;
4260 if ((fPhd[2] > fPhd[0]) &&
4261 (fPhd[2] > fPhd[1]) &&
4262 (fPhd[1] > fPhd[0]) &&
4264 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4265 fNumberFitSuccess++;
4267 if (fPhdt0 >= 0.0) {
4268 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4269 if (fCurrentCoef2[0] < -1.0) {
4270 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4274 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4279 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4280 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4283 if (fDebugLevel == 1) {
4284 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4287 line->SetLineColor(2);
4288 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4289 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4290 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4291 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4292 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4293 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4294 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4295 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4298 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4307 //_____________________________________________________________________________
4308 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4311 // Slope methode but with polynomes de Lagrange
4315 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4318 //Double_t *x = new Double_t[5];
4319 //Double_t *y = new Double_t[5];
4336 fCurrentCoefE = 0.0;
4337 fCurrentCoefE2 = 1.0;
4338 fCurrentCoef[0] = 0.0;
4339 fCurrentCoef2[0] = 0.0;
4340 TLine *line = new TLine();
4341 TF1 * polynome = 0x0;
4342 TF1 * polynomea = 0x0;
4343 TF1 * polynomeb = 0x0;
4351 TAxis *xpph = projPH->GetXaxis();
4352 Int_t nbins = xpph->GetNbins();
4353 Double_t lowedge = xpph->GetBinLowEdge(1);
4354 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4355 Double_t widbins = (upedge - lowedge) / nbins;
4356 Double_t limit = upedge + 0.5 * widbins;
4361 // Beginning of the signal
4362 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4363 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4364 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4367 binmax = (Int_t) pentea->GetMaximumBin();
4369 Double_t minnn = 0.0;
4370 Double_t maxxx = 0.0;
4372 Int_t kase = nbins-binmax;
4380 minnn = pentea->GetBinCenter(binmax-2);
4381 maxxx = pentea->GetBinCenter(binmax);
4382 x[0] = pentea->GetBinCenter(binmax-2);
4383 x[1] = pentea->GetBinCenter(binmax-1);
4384 x[2] = pentea->GetBinCenter(binmax);
4385 y[0] = pentea->GetBinContent(binmax-2);
4386 y[1] = pentea->GetBinContent(binmax-1);
4387 y[2] = pentea->GetBinContent(binmax);
4388 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4389 AliInfo("At the limit for beginning!");
4392 minnn = pentea->GetBinCenter(binmax-2);
4393 maxxx = pentea->GetBinCenter(binmax+1);
4394 x[0] = pentea->GetBinCenter(binmax-2);
4395 x[1] = pentea->GetBinCenter(binmax-1);
4396 x[2] = pentea->GetBinCenter(binmax);
4397 x[3] = pentea->GetBinCenter(binmax+1);
4398 y[0] = pentea->GetBinContent(binmax-2);
4399 y[1] = pentea->GetBinContent(binmax-1);
4400 y[2] = pentea->GetBinContent(binmax);
4401 y[3] = pentea->GetBinContent(binmax+1);
4402 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4410 minnn = pentea->GetBinCenter(binmax);
4411 maxxx = pentea->GetBinCenter(binmax+2);
4412 x[0] = pentea->GetBinCenter(binmax);
4413 x[1] = pentea->GetBinCenter(binmax+1);
4414 x[2] = pentea->GetBinCenter(binmax+2);
4415 y[0] = pentea->GetBinContent(binmax);
4416 y[1] = pentea->GetBinContent(binmax+1);
4417 y[2] = pentea->GetBinContent(binmax+2);
4418 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4421 minnn = pentea->GetBinCenter(binmax-1);
4422 maxxx = pentea->GetBinCenter(binmax+2);
4423 x[0] = pentea->GetBinCenter(binmax-1);
4424 x[1] = pentea->GetBinCenter(binmax);
4425 x[2] = pentea->GetBinCenter(binmax+1);
4426 x[3] = pentea->GetBinCenter(binmax+2);
4427 y[0] = pentea->GetBinContent(binmax-1);
4428 y[1] = pentea->GetBinContent(binmax);
4429 y[2] = pentea->GetBinContent(binmax+1);
4430 y[3] = pentea->GetBinContent(binmax+2);
4431 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4434 minnn = pentea->GetBinCenter(binmax-2);
4435 maxxx = pentea->GetBinCenter(binmax+2);
4436 x[0] = pentea->GetBinCenter(binmax-2);
4437 x[1] = pentea->GetBinCenter(binmax-1);
4438 x[2] = pentea->GetBinCenter(binmax);
4439 x[3] = pentea->GetBinCenter(binmax+1);
4440 x[4] = pentea->GetBinCenter(binmax+2);
4441 y[0] = pentea->GetBinContent(binmax-2);
4442 y[1] = pentea->GetBinContent(binmax-1);
4443 y[2] = pentea->GetBinContent(binmax);
4444 y[3] = pentea->GetBinContent(binmax+1);
4445 y[4] = pentea->GetBinContent(binmax+2);
4446 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4454 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4455 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4457 Double_t step = (maxxx-minnn)/10000;
4459 Double_t maxvalue = 0.0;
4460 Double_t placemaximum = minnn;
4461 for(Int_t o = 0; o < 10000; o++){
4462 if(o == 0) maxvalue = polynomeb->Eval(l);
4463 if(maxvalue < (polynomeb->Eval(l))){
4464 maxvalue = polynomeb->Eval(l);
4469 fPhd[0] = placemaximum;
4472 // Amplification region
4475 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4476 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4482 Double_t minn = 0.0;
4483 Double_t maxx = 0.0;
4495 Int_t kase1 = nbins - binmax;
4497 //Determination of minn and maxx
4498 //case binmax = nbins
4503 minn = projPH->GetBinCenter(binmax-2);
4504 maxx = projPH->GetBinCenter(binmax);
4505 x[0] = projPH->GetBinCenter(binmax-2);
4506 x[1] = projPH->GetBinCenter(binmax-1);
4507 x[2] = projPH->GetBinCenter(binmax);
4508 y[0] = projPH->GetBinContent(binmax-2);
4509 y[1] = projPH->GetBinContent(binmax-1);
4510 y[2] = projPH->GetBinContent(binmax);
4511 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4512 //AliInfo("At the limit for the drift!");
4515 minn = projPH->GetBinCenter(binmax-2);
4516 maxx = projPH->GetBinCenter(binmax+1);
4517 x[0] = projPH->GetBinCenter(binmax-2);
4518 x[1] = projPH->GetBinCenter(binmax-1);
4519 x[2] = projPH->GetBinCenter(binmax);
4520 x[3] = projPH->GetBinCenter(binmax+1);
4521 y[0] = projPH->GetBinContent(binmax-2);
4522 y[1] = projPH->GetBinContent(binmax-1);
4523 y[2] = projPH->GetBinContent(binmax);
4524 y[3] = projPH->GetBinContent(binmax+1);
4525 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4534 minn = projPH->GetBinCenter(binmax);
4535 maxx = projPH->GetBinCenter(binmax+2);
4536 x[0] = projPH->GetBinCenter(binmax);
4537 x[1] = projPH->GetBinCenter(binmax+1);
4538 x[2] = projPH->GetBinCenter(binmax+2);
4539 y[0] = projPH->GetBinContent(binmax);
4540 y[1] = projPH->GetBinContent(binmax+1);
4541 y[2] = projPH->GetBinContent(binmax+2);
4542 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4545 minn = projPH->GetBinCenter(binmax-1);
4546 maxx = projPH->GetBinCenter(binmax+2);
4547 x[0] = projPH->GetBinCenter(binmax-1);
4548 x[1] = projPH->GetBinCenter(binmax);
4549 x[2] = projPH->GetBinCenter(binmax+1);
4550 x[3] = projPH->GetBinCenter(binmax+2);
4551 y[0] = projPH->GetBinContent(binmax-1);
4552 y[1] = projPH->GetBinContent(binmax);
4553 y[2] = projPH->GetBinContent(binmax+1);
4554 y[3] = projPH->GetBinContent(binmax+2);
4555 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4558 minn = projPH->GetBinCenter(binmax-2);
4559 maxx = projPH->GetBinCenter(binmax+2);
4560 x[0] = projPH->GetBinCenter(binmax-2);
4561 x[1] = projPH->GetBinCenter(binmax-1);
4562 x[2] = projPH->GetBinCenter(binmax);
4563 x[3] = projPH->GetBinCenter(binmax+1);
4564 x[4] = projPH->GetBinCenter(binmax+2);
4565 y[0] = projPH->GetBinContent(binmax-2);
4566 y[1] = projPH->GetBinContent(binmax-1);
4567 y[2] = projPH->GetBinContent(binmax);
4568 y[3] = projPH->GetBinContent(binmax+1);
4569 y[4] = projPH->GetBinContent(binmax+2);
4570 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4577 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4578 polynomea->SetParameters(c0,c1,c2,c3,c4);
4580 Double_t step = (maxx-minn)/1000;
4582 Double_t maxvalue = 0.0;
4583 Double_t placemaximum = minn;
4584 for(Int_t o = 0; o < 1000; o++){
4585 if(o == 0) maxvalue = polynomea->Eval(l);
4586 if(maxvalue < (polynomea->Eval(l))){
4587 maxvalue = polynomea->Eval(l);
4592 fPhd[1] = placemaximum;
4596 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4597 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4598 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4601 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4607 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4611 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4612 AliInfo("Too many fluctuations at the end!");
4615 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4616 AliInfo("Too many fluctuations at the end!");
4619 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4620 AliInfo("No entries for the next bin!");
4621 pente->SetBinContent(binmin,0);
4622 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4638 Bool_t case1 = kFALSE;
4639 Bool_t case2 = kFALSE;
4640 Bool_t case4 = kFALSE;
4642 //Determination of min and max
4643 //case binmin <= nbins-3
4645 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4646 min = pente->GetBinCenter(binmin-2);
4647 max = pente->GetBinCenter(binmin+2);
4648 x[0] = pente->GetBinCenter(binmin-2);
4649 x[1] = pente->GetBinCenter(binmin-1);
4650 x[2] = pente->GetBinCenter(binmin);
4651 x[3] = pente->GetBinCenter(binmin+1);
4652 x[4] = pente->GetBinCenter(binmin+2);
4653 y[0] = pente->GetBinContent(binmin-2);
4654 y[1] = pente->GetBinContent(binmin-1);
4655 y[2] = pente->GetBinContent(binmin);
4656 y[3] = pente->GetBinContent(binmin+1);
4657 y[4] = pente->GetBinContent(binmin+2);
4658 //Calcul the polynome de Lagrange
4659 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4661 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4662 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4663 //AliInfo("polynome 4 false 1");
4666 if(((binmin+3) <= (nbins-1)) &&
4667 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4668 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4669 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4670 AliInfo("polynome 4 false 2");
4674 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4675 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4676 //AliInfo("polynome 4 case 1");
4679 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4680 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4681 //AliInfo("polynome 4 case 4");
4686 //case binmin = nbins-2
4688 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4690 min = pente->GetBinCenter(binmin-2);
4691 max = pente->GetBinCenter(binmin+1);
4692 x[0] = pente->GetBinCenter(binmin-2);
4693 x[1] = pente->GetBinCenter(binmin-1);
4694 x[2] = pente->GetBinCenter(binmin);
4695 x[3] = pente->GetBinCenter(binmin+1);
4696 y[0] = pente->GetBinContent(binmin-2);
4697 y[1] = pente->GetBinContent(binmin-1);
4698 y[2] = pente->GetBinContent(binmin);
4699 y[3] = pente->GetBinContent(binmin+1);
4700 //Calcul the polynome de Lagrange
4701 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4702 //richtung +: nothing
4704 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4705 //AliInfo("polynome 3- case 2");
4710 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4712 min = pente->GetBinCenter(binmin-1);
4713 max = pente->GetBinCenter(binmin+2);
4714 x[0] = pente->GetBinCenter(binmin-1);
4715 x[1] = pente->GetBinCenter(binmin);
4716 x[2] = pente->GetBinCenter(binmin+1);
4717 x[3] = pente->GetBinCenter(binmin+2);
4718 y[0] = pente->GetBinContent(binmin-1);
4719 y[1] = pente->GetBinContent(binmin);
4720 y[2] = pente->GetBinContent(binmin+1);
4721 y[3] = pente->GetBinContent(binmin+2);
4722 //Calcul the polynome de Lagrange
4723 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4725 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4726 //AliInfo("polynome 3+ case 2");
4731 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4732 min = pente->GetBinCenter(binmin);
4733 max = pente->GetBinCenter(binmin+2);
4734 x[0] = pente->GetBinCenter(binmin);
4735 x[1] = pente->GetBinCenter(binmin+1);
4736 x[2] = pente->GetBinCenter(binmin+2);
4737 y[0] = pente->GetBinContent(binmin);
4738 y[1] = pente->GetBinContent(binmin+1);
4739 y[2] = pente->GetBinContent(binmin+2);
4740 //Calcul the polynome de Lagrange
4741 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4743 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4744 //AliInfo("polynome 2+ false");
4749 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4751 min = pente->GetBinCenter(binmin-1);
4752 max = pente->GetBinCenter(binmin+1);
4753 x[0] = pente->GetBinCenter(binmin-1);
4754 x[1] = pente->GetBinCenter(binmin);
4755 x[2] = pente->GetBinCenter(binmin+1);
4756 y[0] = pente->GetBinContent(binmin-1);
4757 y[1] = pente->GetBinContent(binmin);
4758 y[2] = pente->GetBinContent(binmin+1);
4759 //Calcul the polynome de Lagrange
4760 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4761 //richtung +: nothing
4762 //richtung -: nothing
4764 //case binmin = nbins-1
4766 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4767 min = pente->GetBinCenter(binmin-2);
4768 max = pente->GetBinCenter(binmin);
4769 x[0] = pente->GetBinCenter(binmin-2);
4770 x[1] = pente->GetBinCenter(binmin-1);
4771 x[2] = pente->GetBinCenter(binmin);
4772 y[0] = pente->GetBinContent(binmin-2);
4773 y[1] = pente->GetBinContent(binmin-1);
4774 y[2] = pente->GetBinContent(binmin);
4775 //Calcul the polynome de Lagrange
4776 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4777 //AliInfo("At the limit for the drift!");
4778 //fluctuation too big!
4779 //richtung +: nothing
4781 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4782 //AliInfo("polynome 2- false ");
4786 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4788 AliInfo("At the limit for the drift and not usable!");
4792 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4794 AliInfo("For the drift...problem!");
4796 //pass but should not happen
4797 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4799 AliInfo("For the drift...problem!");
4803 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4804 polynome->SetParameters(c0,c1,c2,c3,c4);
4805 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4806 Double_t step = (max-min)/1000;
4808 Double_t minvalue = 0.0;
4809 Double_t placeminimum = min;
4810 for(Int_t o = 0; o < 1000; o++){
4811 if(o == 0) minvalue = polynome->Eval(l);
4812 if(minvalue > (polynome->Eval(l))){
4813 minvalue = polynome->Eval(l);
4818 fPhd[2] = placeminimum;
4820 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4821 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4822 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4824 Float_t fPhdt0 = 0.0;
4825 Float_t t0Shift = 0.0;
4828 t0Shift = fT0Shift1;
4832 t0Shift = fT0Shift0;
4835 if ((fPhd[2] > fPhd[0]) &&
4836 (fPhd[2] > fPhd[1]) &&
4837 (fPhd[1] > fPhd[0]) &&
4839 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4840 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4841 else fNumberFitSuccess++;
4842 if (fPhdt0 >= 0.0) {
4843 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4844 if (fCurrentCoef2[0] < -1.0) {
4845 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4849 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4853 //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4854 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4856 if((fPhd[1] > fPhd[0]) &&
4858 if (fPhdt0 >= 0.0) {
4859 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4860 if (fCurrentCoef2[0] < -1.0) {
4861 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4865 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4869 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4870 //printf("Fit failed!\n");
4874 if (fDebugLevel == 1) {
4875 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4878 line->SetLineColor(2);
4879 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4880 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4881 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4882 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4883 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4884 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4885 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4886 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4889 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4896 if(polynome) delete polynome;
4897 if(polynomea) delete polynomea;
4898 if(polynomeb) delete polynomeb;
4899 //if(x) delete [] x;
4900 //if(y) delete [] y;
4901 if(line) delete line;
4906 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4907 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4909 projPH->SetDirectory(0);
4913 //_____________________________________________________________________________
4914 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4917 // Fit methode for the drift velocity
4921 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4924 TAxis *xpph = projPH->GetXaxis();
4925 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4927 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4928 fPH->SetParameter(0,0.469); // Scaling
4929 fPH->SetParameter(1,0.18); // Start
4930 fPH->SetParameter(2,0.0857325); // AR
4931 fPH->SetParameter(3,1.89); // DR
4932 fPH->SetParameter(4,0.08); // QA/QD
4933 fPH->SetParameter(5,0.0); // Baseline
4935 TLine *line = new TLine();
4937 fCurrentCoef[0] = 0.0;
4938 fCurrentCoef2[0] = 0.0;
4939 fCurrentCoefE = 0.0;
4940 fCurrentCoefE2 = 0.0;
4942 if (idect%fFitPHPeriode == 0) {
4944 AliInfo(Form("The detector %d will be fitted",idect));
4945 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4946 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4947 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4948 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4949 fPH->SetParameter(4,0.225); // QA/QD
4950 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4952 if (fDebugLevel != 1) {
4953 projPH->Fit(fPH,"0M","",0.0,upedge);
4956 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4958 projPH->Fit(fPH,"M+","",0.0,upedge);
4960 line->SetLineColor(4);
4961 line->DrawLine(fPH->GetParameter(1)
4963 ,fPH->GetParameter(1)
4964 ,projPH->GetMaximum());
4965 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4967 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4968 ,projPH->GetMaximum());
4969 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4971 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4972 ,projPH->GetMaximum());
4975 if (fPH->GetParameter(3) != 0) {
4976 fNumberFitSuccess++;
4977 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4978 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4979 fCurrentCoef2[0] = fPH->GetParameter(1);
4980 fCurrentCoefE2 = fPH->GetParError(1);
4983 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4984 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4990 // Put the default value
4991 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4992 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4995 if (fDebugLevel != 1) {
5000 //_____________________________________________________________________________
5001 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5004 // Fit methode for the sigma of the pad response function
5009 fCurrentCoef[0] = 0.0;
5010 fCurrentCoefE = 0.0;
5012 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5014 if(TMath::Abs(ret+4) <= 0.000000001){
5015 fCurrentCoef[0] = -fCurrentCoef[1];
5019 fNumberFitSuccess++;
5020 fCurrentCoef[0] = param[2];
5021 fCurrentCoefE = ret;
5025 //_____________________________________________________________________________
5026 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)
5029 // Fit methode for the sigma of the pad response function
5032 //We should have at least 3 points
5033 if(nBins <=3) return -4.0;
5035 TLinearFitter fitter(3,"pol2");
5036 fitter.StoreData(kFALSE);
5037 fitter.ClearPoints();
5039 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5040 Float_t entries = 0;
5041 Int_t nbbinwithentries = 0;
5042 for (Int_t i=0; i<nBins; i++){
5044 if(arraye[i] > 15) nbbinwithentries++;
5045 //printf("entries for i %d: %f\n",i,arraye[i]);
5047 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5048 //printf("entries %f\n",entries);
5049 //printf("nbbinwithentries %d\n",nbbinwithentries);
5052 Float_t errorm = 0.0;
5053 Float_t errorn = 0.0;
5054 Float_t error = 0.0;
5057 for (Int_t ibin=0;ibin<nBins; ibin++){
5058 Float_t entriesI = arraye[ibin];
5059 Float_t valueI = arraym[ibin];
5060 Double_t xcenter = 0.0;
5062 if ((entriesI>15) && (valueI>0.0)){
5063 xcenter = xMin+(ibin+0.5)*binWidth;
5068 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5069 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5070 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5073 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5074 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5075 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5077 if(TMath::Abs(error) < 0.000000001) continue;
5078 val = TMath::Log(Float_t(valueI));
5079 fitter.AddPoint(&xcenter,val,error);
5083 if(fDebugLevel > 1){
5085 if ( !fDebugStreamer ) {
5087 TDirectory *backup = gDirectory;
5088 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5089 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5092 Int_t detector = fCountDet;
5093 Int_t layer = GetLayer(fCountDet);
5096 (* fDebugStreamer) << "FitGausMIFill"<<
5097 "detector="<<detector<<
5101 "entriesI="<<entriesI<<
5104 "xcenter="<<xcenter<<
5114 if(npoints <=3) return -4.0;
5119 fitter.GetParameters(par);
5120 chi2 = fitter.GetChisquare()/Float_t(npoints);
5123 if (!param) param = new TVectorD(3);
5124 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5125 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5126 Double_t deltax = (fitter.GetParError(2))/x;
5127 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5130 (*param)[1] = par[1]/(-2.*par[2]);
5131 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5132 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5133 if ( lnparam0>307 ) return -4;
5134 (*param)[0] = TMath::Exp(lnparam0);
5136 if(fDebugLevel > 1){
5138 if ( !fDebugStreamer ) {
5140 TDirectory *backup = gDirectory;
5141 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5142 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5145 Int_t detector = fCountDet;
5146 Int_t layer = GetLayer(fCountDet);
5149 (* fDebugStreamer) << "FitGausMIFit"<<
5150 "detector="<<detector<<
5153 "errorsigma="<<chi2<<
5154 "mean="<<(*param)[1]<<
5155 "sigma="<<(*param)[2]<<
5156 "constant="<<(*param)[0]<<
5161 if((chi2/(*param)[2]) > 0.1){
5163 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5168 if(fDebugLevel == 1){
5169 TString name("PRF");
5170 name += (Int_t)xMin;
5171 name += (Int_t)xMax;
5172 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5175 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5176 for(Int_t k = 0; k < nBins; k++){
5177 histo->SetBinContent(k+1,arraym[k]);
5178 histo->SetBinError(k+1,arrayme[k]);
5181 name += "functionf";
5182 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5183 f1->SetParameter(0, (*param)[0]);
5184 f1->SetParameter(1, (*param)[1]);
5185 f1->SetParameter(2, (*param)[2]);
5193 //_____________________________________________________________________________
5194 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5197 // Fit methode for the sigma of the pad response function
5200 fCurrentCoef[0] = 0.0;
5201 fCurrentCoefE = 0.0;
5203 if (fDebugLevel != 1) {
5204 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5207 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5209 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5212 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5213 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5214 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5216 fNumberFitSuccess++;
5219 //_____________________________________________________________________________
5220 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5223 // Fit methode for the sigma of the pad response function
5225 fCurrentCoef[0] = 0.0;
5226 fCurrentCoefE = 0.0;
5227 if (fDebugLevel == 1) {
5228 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5232 fCurrentCoef[0] = projPRF->GetRMS();
5233 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5235 fNumberFitSuccess++;
5238 //_____________________________________________________________________________
5239 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5242 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5245 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5248 Int_t nbins = (Int_t)(nybins/(2*nbg));
5249 Float_t lowedge = -3.0*nbg;
5250 Float_t upedge = lowedge + 3.0;
5253 Double_t xvalues = -0.2*nbg+0.1;
5255 Int_t total = 2*nbg;
5258 for(Int_t k = 0; k < total; k++){
5259 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5261 y = fCurrentCoef[0]*fCurrentCoef[0];
5262 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5265 if(fDebugLevel > 1){
5267 if ( !fDebugStreamer ) {
5269 TDirectory *backup = gDirectory;
5270 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5271 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5274 Int_t detector = fCountDet;
5275 Int_t layer = GetLayer(fCountDet);
5276 Int_t nbtotal = total;
5278 Float_t low = lowedge;
5279 Float_t up = upedge;
5280 Float_t tnp = xvalues;
5281 Float_t wid = fCurrentCoef[0];
5282 Float_t widfE = fCurrentCoefE;
5284 (* fDebugStreamer) << "FitTnpRange0"<<
5285 "detector="<<detector<<
5287 "nbtotal="<<nbtotal<<
5305 fCurrentCoefE = 0.0;
5306 fCurrentCoef[0] = 0.0;
5308 //printf("npoints\n",npoints);
5311 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5316 linearfitter.Eval();
5317 linearfitter.GetParameters(pars0);
5318 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5319 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5320 Double_t min0 = 0.0;
5321 Double_t ermin0 = 0.0;
5322 //Double_t prfe0 = 0.0;
5323 Double_t prf0 = 0.0;
5324 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5325 min0 = -pars0[1]/(2*pars0[2]);
5326 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5327 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5330 prfe0 = linearfitter->GetParError(0)*pointError0
5331 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5332 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5333 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5334 fCurrentCoefE = (Float_t) prfe0;
5336 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5339 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5343 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5346 if(fDebugLevel > 1){
5348 if ( !fDebugStreamer ) {
5350 TDirectory *backup = gDirectory;
5351 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5352 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5355 Int_t detector = fCountDet;
5356 Int_t layer = GetLayer(fCountDet);
5357 Int_t nbtotal = total;
5358 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5359 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5361 (* fDebugStreamer) << "FitTnpRange1"<<
5362 "detector="<<detector<<
5364 "nbtotal="<<nbtotal<<
5368 "npoints="<<npoints<<
5371 "sigmaprf="<<fCurrentCoef[0]<<
5372 "sigprf="<<fCurrentCoef[1]<<
5379 //_____________________________________________________________________________
5380 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5383 // Only mean methode for the gain factor
5386 fCurrentCoef[0] = mean;
5387 fCurrentCoefE = 0.0;
5388 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5389 if (fDebugLevel == 1) {
5390 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5394 CalculChargeCoefMean(kTRUE);
5395 fNumberFitSuccess++;
5397 //_____________________________________________________________________________
5398 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5401 // mean w methode for the gain factor
5405 Int_t nybins = projch->GetNbinsX();
5407 //The weight function
5408 Double_t a = 0.00228515;
5409 Double_t b = -0.00231487;
5410 Double_t c = 0.00044298;
5411 Double_t d = -0.00379239;
5412 Double_t e = 0.00338349;
5422 //A arbitrary error for the moment
5423 fCurrentCoefE = 0.0;
5424 fCurrentCoef[0] = 0.0;
5427 Double_t sumw = 0.0;
5429 Float_t sumAll = (Float_t) nentries;
5430 Int_t sumCurrent = 0;
5431 for(Int_t k = 0; k <nybins; k++){
5432 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5433 if (fraction>0.95) break;
5434 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5435 e*fraction*fraction*fraction*fraction;
5436 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5437 sum += weight*projch->GetBinContent(k+1);
5438 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5439 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5441 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5443 if (fDebugLevel == 1) {
5444 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5448 fNumberFitSuccess++;
5449 CalculChargeCoefMean(kTRUE);
5451 //_____________________________________________________________________________
5452 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5455 // mean w methode for the gain factor
5459 Int_t nybins = projch->GetNbinsX();
5461 //The weight function
5462 Double_t a = 0.00228515;
5463 Double_t b = -0.00231487;
5464 Double_t c = 0.00044298;
5465 Double_t d = -0.00379239;
5466 Double_t e = 0.00338349;
5476 //A arbitrary error for the moment
5477 fCurrentCoefE = 0.0;
5478 fCurrentCoef[0] = 0.0;
5481 Double_t sumw = 0.0;
5483 Int_t sumCurrent = 0;
5484 for(Int_t k = 0; k <nybins; k++){
5485 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5486 if (fraction>0.95) break;
5487 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5488 e*fraction*fraction*fraction*fraction;
5489 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5490 sum += weight*projch->GetBinContent(k+1);
5491 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5492 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5494 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5496 if (fDebugLevel == 1) {
5497 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5501 fNumberFitSuccess++;
5503 //_____________________________________________________________________________
5504 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5507 // Fit methode for the gain factor
5510 fCurrentCoef[0] = 0.0;
5511 fCurrentCoefE = 0.0;
5512 Double_t chisqrl = 0.0;
5513 Double_t chisqrg = 0.0;
5514 Double_t chisqr = 0.0;
5515 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5517 projch->Fit("landau","0",""
5518 ,(Double_t) mean/fBeginFitCharge
5519 ,projch->GetBinCenter(projch->GetNbinsX()));
5520 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5521 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5522 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5523 chisqrl = projch->GetFunction("landau")->GetChisquare();
5525 projch->Fit("gaus","0",""
5526 ,(Double_t) mean/fBeginFitCharge
5527 ,projch->GetBinCenter(projch->GetNbinsX()));
5528 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5529 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5530 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5532 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5533 if (fDebugLevel != 1) {
5534 projch->Fit("fLandauGaus","0",""
5535 ,(Double_t) mean/fBeginFitCharge
5536 ,projch->GetBinCenter(projch->GetNbinsX()));
5537 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5540 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5542 projch->Fit("fLandauGaus","+",""
5543 ,(Double_t) mean/fBeginFitCharge
5544 ,projch->GetBinCenter(projch->GetNbinsX()));
5545 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5547 fLandauGaus->Draw("same");
5550 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5551 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5552 fNumberFitSuccess++;
5553 CalculChargeCoefMean(kTRUE);
5554 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5555 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5558 CalculChargeCoefMean(kFALSE);
5559 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5562 if (fDebugLevel != 1) {
5567 //_____________________________________________________________________________
5568 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5571 // Fit methode for the gain factor more time consuming
5575 //Some parameters to initialise
5576 Double_t widthLandau, widthGaus, mPV, integral;
5577 Double_t chisquarel = 0.0;
5578 Double_t chisquareg = 0.0;
5579 projch->Fit("landau","0M+",""
5581 ,projch->GetBinCenter(projch->GetNbinsX()));
5582 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5583 chisquarel = projch->GetFunction("landau")->GetChisquare();
5584 projch->Fit("gaus","0M+",""
5586 ,projch->GetBinCenter(projch->GetNbinsX()));
5587 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5588 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5590 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5591 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5593 // Setting fit range and start values
5595 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5596 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5597 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5598 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5599 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5600 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5601 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5604 fCurrentCoef[0] = 0.0;
5605 fCurrentCoefE = 0.0;
5609 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5614 Double_t projchPeak;
5615 Double_t projchFWHM;
5616 LanGauPro(fp,projchPeak,projchFWHM);
5618 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5619 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5620 fNumberFitSuccess++;
5621 CalculChargeCoefMean(kTRUE);
5622 fCurrentCoef[0] = fp[1];
5623 fCurrentCoefE = fpe[1];
5624 //chargeCoefE2 = chisqr;
5627 CalculChargeCoefMean(kFALSE);
5628 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5630 if (fDebugLevel == 1) {
5631 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5632 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5635 fitsnr->Draw("same");
5641 //_____________________________________________________________________________
5642 void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5645 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5647 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5648 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5649 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5654 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5655 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5659 //_____________________________________________________________________________
5660 void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5663 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5665 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5666 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5667 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5668 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5672 c2 = -(x0*(x[1]+x[2]+x[3])
5673 +x1*(x[0]+x[2]+x[3])
5674 +x2*(x[0]+x[1]+x[3])
5675 +x3*(x[0]+x[1]+x[2]));
5676 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5677 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5678 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5679 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5681 c0 = -(x0*x[1]*x[2]*x[3]
5684 +x3*x[0]*x[1]*x[2]);
5689 //_____________________________________________________________________________
5690 void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
5693 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5695 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5696 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5697 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5698 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5699 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5702 c4 = x0+x1+x2+x3+x4;
5703 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5704 +x1*(x[0]+x[2]+x[3]+x[4])
5705 +x2*(x[0]+x[1]+x[3]+x[4])
5706 +x3*(x[0]+x[1]+x[2]+x[4])
5707 +x4*(x[0]+x[1]+x[2]+x[3]));
5708 c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5709 +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])
5710 +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])
5711 +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])
5712 +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]));
5714 c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5715 +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])
5716 +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])
5717 +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])
5718 +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]));
5720 c0 = (x0*x[1]*x[2]*x[3]*x[4]
5721 +x1*x[0]*x[2]*x[3]*x[4]
5722 +x2*x[0]*x[1]*x[3]*x[4]
5723 +x3*x[0]*x[1]*x[2]*x[4]
5724 +x4*x[0]*x[1]*x[2]*x[3]);
5727 //_____________________________________________________________________________
5728 void AliTRDCalibraFit::NormierungCharge()
5731 // Normalisation of the gain factor resulting for the fits
5734 // Calcul of the mean of choosen method by fFitChargeNDB
5736 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5737 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5739 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5740 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5741 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5742 if (GetStack(detector) == 2) {
5745 if (GetStack(detector) != 2) {
5748 for (Int_t j = 0; j < total; j++) {
5756 fScaleFitFactor = fScaleFitFactor / sum;
5759 fScaleFitFactor = 1.0;
5762 //methode de boeuf mais bon...
5763 Double_t scalefactor = fScaleFitFactor;
5765 if(fDebugLevel > 1){
5767 if ( !fDebugStreamer ) {
5769 TDirectory *backup = gDirectory;
5770 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5771 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5773 (* fDebugStreamer) << "NormierungCharge"<<
5774 "scalefactor="<<scalefactor<<
5778 //_____________________________________________________________________________
5779 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5782 // Rebin of the 1D histo for the gain calibration if needed.
5783 // you have to choose fRebin, divider of fNumberBinCharge
5786 TAxis *xhist = hist->GetXaxis();
5787 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5788 ,xhist->GetBinLowEdge(1)
5789 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5791 AliInfo(Form("fRebin: %d",fRebin));
5793 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5795 for (Int_t ji = i; ji < i+fRebin; ji++) {
5796 sum += hist->GetBinContent(ji);
5799 rehist->SetBinContent(k,sum);
5807 //_____________________________________________________________________________
5808 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5811 // Rebin of the 1D histo for the gain calibration if needed
5812 // you have to choose fRebin divider of fNumberBinCharge
5815 TAxis *xhist = hist->GetXaxis();
5816 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5817 ,xhist->GetBinLowEdge(1)
5818 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5820 AliInfo(Form("fRebin: %d",fRebin));
5822 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5824 for (Int_t ji = i; ji < i+fRebin; ji++) {
5825 sum += hist->GetBinContent(ji);
5828 rehist->SetBinContent(k,sum);
5836 //____________Some basic geometry function_____________________________________
5839 //_____________________________________________________________________________
5840 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5843 // Reconstruct the plane number from the detector number
5846 return ((Int_t) (d % 6));
5850 //_____________________________________________________________________________
5851 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5854 // Reconstruct the stack number from the detector number
5856 const Int_t kNlayer = 6;
5858 return ((Int_t) (d % 30) / kNlayer);
5862 //_____________________________________________________________________________
5863 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5866 // Reconstruct the sector number from the detector number
5870 return ((Int_t) (d / fg));
5875 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5877 //_______________________________________________________________________________
5878 void AliTRDCalibraFit::ResetVectorFit()
5881 // Reset the VectorFits
5884 fVectorFit.SetOwner();
5886 fVectorFit2.SetOwner();
5887 fVectorFit2.Clear();
5891 //____________Private Functions________________________________________________
5894 //_____________________________________________________________________________
5895 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5898 // Function for the fit
5901 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5903 //PARAMETERS FOR FIT PH
5905 //fAsymmGauss->SetParameter(0,0.113755);
5906 //fAsymmGauss->SetParameter(1,0.350706);
5907 //fAsymmGauss->SetParameter(2,0.0604244);
5908 //fAsymmGauss->SetParameter(3,7.65596);
5909 //fAsymmGauss->SetParameter(4,1.00124);
5910 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5918 Double_t dx = 0.005;
5919 Double_t xs = par[1];
5921 Double_t paras[2] = { 0.0, 0.0 };
5924 if ((xs >= par[1]) &&
5925 (xs < (par[1]+par[2]))) {
5926 //fAsymmGauss->SetParameter(0,par[0]);
5927 //fAsymmGauss->SetParameter(1,xs);
5928 //ss += fAsymmGauss->Eval(xx);
5931 ss += AsymmGauss(&xx,paras);
5933 if ((xs >= (par[1]+par[2])) &&
5934 (xs < (par[1]+par[2]+par[3]))) {
5935 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5936 //fAsymmGauss->SetParameter(1,xs);
5937 //ss += fAsymmGauss->Eval(xx);
5938 paras[0] = par[0]*par[4];
5940 ss += AsymmGauss(&xx,paras);
5949 //_____________________________________________________________________________
5950 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5953 // Function for the fit
5956 //par[0] = normalization
5964 Double_t par1save = par[1];
5965 //Double_t par2save = par[2];
5966 Double_t par2save = 0.0604244;
5967 //Double_t par3save = par[3];
5968 Double_t par3save = 7.65596;
5969 //Double_t par5save = par[5];
5970 Double_t par5save = 0.870597;
5971 Double_t dx = x[0] - par1save;
5973 Double_t sigma2 = par2save*par2save;
5974 Double_t sqrt2 = TMath::Sqrt(2.0);
5975 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5976 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5977 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5978 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5980 //return par[0]*(exp1+par[4]*exp2);
5981 return par[0] * (exp1 + 1.00124 * exp2);
5985 //_____________________________________________________________________________
5986 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5989 // Sum Landau + Gaus with identical mean
5992 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5993 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5994 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5995 Double_t val = valLandau + valGaus;
6001 //_____________________________________________________________________________
6002 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6005 // Function for the fit
6008 // par[0]=Width (scale) parameter of Landau density
6009 // par[1]=Most Probable (MP, location) parameter of Landau density
6010 // par[2]=Total area (integral -inf to inf, normalization constant)
6011 // par[3]=Width (sigma) of convoluted Gaussian function
6013 // In the Landau distribution (represented by the CERNLIB approximation),
6014 // the maximum is located at x=-0.22278298 with the location parameter=0.
6015 // This shift is corrected within this function, so that the actual
6016 // maximum is identical to the MP parameter.
6019 // Numeric constants
6020 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6021 Double_t mpshift = -0.22278298; // Landau maximum location
6023 // Control constants
6024 Double_t np = 100.0; // Number of convolution steps
6025 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6037 // MP shift correction
6038 mpc = par[1] - mpshift * par[0];
6040 // Range of convolution integral
6041 xlow = x[0] - sc * par[3];
6042 xupp = x[0] + sc * par[3];
6044 step = (xupp - xlow) / np;
6046 // Convolution integral of Landau and Gaussian by sum
6047 for (i = 1.0; i <= np/2; i++) {
6049 xx = xlow + (i-.5) * step;
6050 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6051 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6053 xx = xupp - (i-.5) * step;
6054 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6055 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6059 return (par[2] * step * sum * invsq2pi / par[3]);
6062 //_____________________________________________________________________________
6063 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6064 , const Double_t *parlimitslo, const Double_t *parlimitshi
6065 , Double_t *fitparams, Double_t *fiterrors
6066 , Double_t *chiSqr, Int_t *ndf) const
6069 // Function for the fit
6073 Char_t funname[100];
6075 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6080 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6081 ffit->SetParameters(startvalues);
6082 ffit->SetParNames("Width","MP","Area","GSigma");
6084 for (i = 0; i < 4; i++) {
6085 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6088 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6090 ffit->GetParameters(fitparams); // Obtain fit parameters
6091 for (i = 0; i < 4; i++) {
6092 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6094 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6095 ndf[0] = ffit->GetNDF(); // Obtain ndf
6097 return (ffit); // Return fit function
6101 //_____________________________________________________________________________
6102 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6105 // Function for the fit
6118 Int_t maxcalls = 10000;
6120 // Search for maximum
6121 p = params[1] - 0.1 * params[0];
6122 step = 0.05 * params[0];
6126 while ((l != lold) && (i < maxcalls)) {
6130 l = LanGauFun(&x,params);
6132 step = -step / 10.0;
6137 if (i == maxcalls) {
6143 // Search for right x location of fy
6144 p = maxx + params[0];
6150 while ( (l != lold) && (i < maxcalls) ) {
6155 l = TMath::Abs(LanGauFun(&x,params) - fy);
6169 // Search for left x location of fy
6171 p = maxx - 0.5 * params[0];
6177 while ((l != lold) && (i < maxcalls)) {
6181 l = TMath::Abs(LanGauFun(&x,params) - fy);
6183 step = -step / 10.0;
6188 if (i == maxcalls) {
6197 //_____________________________________________________________________________
6198 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6201 // Gaus with identical mean
6204 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];