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 object->SetValue(detector,mean);
1952 //_____________________________________________________________________________
1953 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1956 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1957 // It takes the min value of the coefficients per detector
1958 // This object has to be written in the database
1961 // Create the DetObject
1962 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1964 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1965 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1966 Int_t detector = -1;
1967 Float_t value = 0.0;
1969 for (Int_t k = 0; k < loop; k++) {
1970 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1971 Float_t min = 100.0;
1973 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1975 if(value > 70.0) value = value-100.0;
1980 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1981 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1982 for (Int_t row = 0; row < rowMax; row++) {
1983 for (Int_t col = 0; col < colMax; col++) {
1984 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1986 if(value > 70.0) value = value-100.0;
1988 if(min > value) min = value;
1992 object->SetValue(detector,min);
1998 //_____________________________________________________________________________
1999 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2002 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2003 // It takes the min value of the coefficients per detector
2004 // This object has to be written in the database
2007 // Create the DetObject
2008 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2011 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2012 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2013 Int_t detector = -1;
2014 Float_t value = 0.0;
2016 for (Int_t k = 0; k < loop; k++) {
2017 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2019 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2020 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2021 Float_t min = 100.0;
2022 for (Int_t row = 0; row < rowMax; row++) {
2023 for (Int_t col = 0; col < colMax; col++) {
2024 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2025 mean += -TMath::Abs(value);
2029 if(count > 0) mean = mean/count;
2031 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2032 object->SetValue(detector,-TMath::Abs(value));
2038 //_____________________________________________________________________________
2039 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2042 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2043 // You need first to create the object for the detectors,
2044 // where the mean value is put.
2045 // This object has to be written in the database
2048 // Create the DetObject
2049 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2052 for(Int_t k = 0; k < 540; k++){
2053 AliTRDCalROC *calROC = object->GetCalROC(k);
2054 Int_t nchannels = calROC->GetNchannels();
2055 for(Int_t ch = 0; ch < nchannels; ch++){
2056 calROC->SetValue(ch,1.0);
2062 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2063 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2064 Int_t detector = -1;
2065 Float_t value = 0.0;
2067 for (Int_t k = 0; k < loop; k++) {
2068 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2069 AliTRDCalROC *calROC = object->GetCalROC(detector);
2070 Float_t mean = detobject->GetValue(detector);
2071 if(TMath::Abs(mean) <= 0.0000000001) continue;
2072 Int_t rowMax = calROC->GetNrows();
2073 Int_t colMax = calROC->GetNcols();
2074 for (Int_t row = 0; row < rowMax; row++) {
2075 for (Int_t col = 0; col < colMax; col++) {
2076 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2077 if(value > 0) value = value*scaleFitFactor;
2078 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2086 //_____________________________________________________________________________
2087 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2090 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2091 // You need first to create the object for the detectors,
2092 // where the mean value is put.
2093 // This object has to be written in the database
2096 // Create the DetObject
2097 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2100 for(Int_t k = 0; k < 540; k++){
2101 AliTRDCalROC *calROC = object->GetCalROC(k);
2102 Int_t nchannels = calROC->GetNchannels();
2103 for(Int_t ch = 0; ch < nchannels; ch++){
2104 calROC->SetValue(ch,1.0);
2110 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2111 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2112 Int_t detector = -1;
2113 Float_t value = 0.0;
2115 for (Int_t k = 0; k < loop; k++) {
2116 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2117 AliTRDCalROC *calROC = object->GetCalROC(detector);
2118 Float_t mean = detobject->GetValue(detector);
2119 if(mean == 0) continue;
2120 Int_t rowMax = calROC->GetNrows();
2121 Int_t colMax = calROC->GetNcols();
2122 for (Int_t row = 0; row < rowMax; row++) {
2123 for (Int_t col = 0; col < colMax; col++) {
2124 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2125 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2133 //_____________________________________________________________________________
2134 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2137 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2138 // You need first to create the object for the detectors,
2139 // where the mean value is put.
2140 // This object has to be written in the database
2143 // Create the DetObject
2144 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2147 for(Int_t k = 0; k < 540; k++){
2148 AliTRDCalROC *calROC = object->GetCalROC(k);
2149 Int_t nchannels = calROC->GetNchannels();
2150 for(Int_t ch = 0; ch < nchannels; ch++){
2151 calROC->SetValue(ch,0.0);
2157 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2158 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2159 Int_t detector = -1;
2160 Float_t value = 0.0;
2162 for (Int_t k = 0; k < loop; k++) {
2163 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2164 AliTRDCalROC *calROC = object->GetCalROC(detector);
2165 Float_t min = detobject->GetValue(detector);
2166 Int_t rowMax = calROC->GetNrows();
2167 Int_t colMax = calROC->GetNcols();
2168 for (Int_t row = 0; row < rowMax; row++) {
2169 for (Int_t col = 0; col < colMax; col++) {
2170 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2172 if(value > 70.0) value = value - 100.0;
2174 calROC->SetValue(col,row,value-min);
2182 //_____________________________________________________________________________
2183 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2186 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2187 // This object has to be written in the database
2190 // Create the DetObject
2191 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2193 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2194 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2195 Int_t detector = -1;
2196 Float_t value = 0.0;
2198 for (Int_t k = 0; k < loop; k++) {
2199 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2200 AliTRDCalROC *calROC = object->GetCalROC(detector);
2201 Int_t rowMax = calROC->GetNrows();
2202 Int_t colMax = calROC->GetNcols();
2203 for (Int_t row = 0; row < rowMax; row++) {
2204 for (Int_t col = 0; col < colMax; col++) {
2205 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2206 calROC->SetValue(col,row,TMath::Abs(value));
2214 //_____________________________________________________________________________
2215 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2218 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2219 // 0 successful fit 1 not successful fit
2220 // mean is the mean value over the successful fit
2221 // do not use it for t0: no meaning
2224 // Create the CalObject
2225 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2229 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2231 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2232 for(Int_t k = 0; k < 540; k++){
2233 object->SetValue(k,1.0);
2236 Int_t detector = -1;
2237 Float_t value = 0.0;
2239 for (Int_t k = 0; k < loop; k++) {
2240 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2241 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2242 if(value <= 0) object->SetValue(detector,1.0);
2244 object->SetValue(detector,0.0);
2249 if(count > 0) mean /= count;
2252 //_____________________________________________________________________________
2253 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2256 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2257 // 0 not successful fit 1 successful fit
2258 // mean mean value over the successful fit
2261 // Create the CalObject
2262 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2266 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2268 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2269 for(Int_t k = 0; k < 540; k++){
2270 AliTRDCalROC *calROC = object->GetCalROC(k);
2271 Int_t nchannels = calROC->GetNchannels();
2272 for(Int_t ch = 0; ch < nchannels; ch++){
2273 calROC->SetValue(ch,1.0);
2277 Int_t detector = -1;
2278 Float_t value = 0.0;
2280 for (Int_t k = 0; k < loop; k++) {
2281 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2282 AliTRDCalROC *calROC = object->GetCalROC(detector);
2283 Int_t nchannels = calROC->GetNchannels();
2284 for (Int_t ch = 0; ch < nchannels; ch++) {
2285 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2286 if(value <= 0) calROC->SetValue(ch,1.0);
2288 calROC->SetValue(ch,0.0);
2294 if(count > 0) mean /= count;
2297 //_____________________________________________________________________________
2298 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2301 // Set FitPH if 1 then each detector will be fitted
2304 if (periodeFitPH > 0) {
2305 fFitPHPeriode = periodeFitPH;
2308 AliInfo("periodeFitPH must be higher than 0!");
2312 //_____________________________________________________________________________
2313 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2316 // The fit of the deposited charge distribution begins at
2317 // histo->Mean()/beginFitCharge
2318 // You can here set beginFitCharge
2321 if (beginFitCharge > 0) {
2322 fBeginFitCharge = beginFitCharge;
2325 AliInfo("beginFitCharge must be strict positif!");
2330 //_____________________________________________________________________________
2331 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2334 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2335 // You can here set t0Shift0
2339 fT0Shift0 = t0Shift;
2342 AliInfo("t0Shift0 must be strict positif!");
2347 //_____________________________________________________________________________
2348 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2351 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2352 // You can here set t0Shift1
2356 fT0Shift1 = t0Shift;
2359 AliInfo("t0Shift must be strict positif!");
2364 //_____________________________________________________________________________
2365 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2368 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2369 // You can here set rangeFitPRF
2372 if ((rangeFitPRF > 0) &&
2373 (rangeFitPRF <= 1.5)) {
2374 fRangeFitPRF = rangeFitPRF;
2377 AliInfo("rangeFitPRF must be between 0 and 1.0");
2382 //_____________________________________________________________________________
2383 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2386 // Minimum entries for fitting
2389 if (minEntries > 0) {
2390 fMinEntries = minEntries;
2393 AliInfo("fMinEntries must be >= 0.");
2398 //_____________________________________________________________________________
2399 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2402 // Rebin with rebin time less bins the Ch histo
2403 // You can set here rebin that should divide the number of bins of CH histo
2408 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2411 AliInfo("You have to choose a positiv value!");
2415 //_____________________________________________________________________________
2416 Bool_t AliTRDCalibraFit::FillVectorFit()
2419 // For the Fit functions fill the vector Fit
2422 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2425 if (GetStack(fCountDet) == 2) {
2432 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2433 Float_t *coef = new Float_t[ntotal];
2434 for (Int_t i = 0; i < ntotal; i++) {
2435 coef[i] = fCurrentCoefDetector[i];
2438 Int_t detector = fCountDet;
2440 fitInfo->SetCoef(coef);
2441 fitInfo->SetDetector(detector);
2442 fVectorFit.Add((TObject *) fitInfo);
2447 //_____________________________________________________________________________
2448 Bool_t AliTRDCalibraFit::FillVectorFit2()
2451 // For the Fit functions fill the vector Fit
2454 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2457 if (GetStack(fCountDet) == 2) {
2464 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2465 Float_t *coef = new Float_t[ntotal];
2466 for (Int_t i = 0; i < ntotal; i++) {
2467 coef[i] = fCurrentCoefDetector2[i];
2470 Int_t detector = fCountDet;
2472 fitInfo->SetCoef(coef);
2473 fitInfo->SetDetector(detector);
2474 fVectorFit2.Add((TObject *) fitInfo);
2479 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2480 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2483 // Init the number of expected bins and fDect1[i] fDect2[i]
2486 gStyle->SetPalette(1);
2487 gStyle->SetOptStat(1111);
2488 gStyle->SetPadBorderMode(0);
2489 gStyle->SetCanvasColor(10);
2490 gStyle->SetPadLeftMargin(0.13);
2491 gStyle->SetPadRightMargin(0.01);
2493 // Mode groups of pads: the total number of bins!
2494 CalculNumberOfBinsExpected(i);
2496 // Quick verification that we have the good pad calibration mode!
2497 if (fNumberOfBinsExpected != nbins) {
2498 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2502 // Security for fDebug 3 and 4
2503 if ((fDebugLevel >= 3) &&
2507 AliInfo("This detector doesn't exit!");
2511 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2512 CalculDect1Dect2(i);
2517 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2518 Bool_t AliTRDCalibraFit::InitFitCH()
2521 // Init the fVectorFitCH for normalisation
2522 // Init the histo for debugging
2527 fScaleFitFactor = 0.0;
2528 fCurrentCoefDetector = new Float_t[2304];
2529 for (Int_t k = 0; k < 2304; k++) {
2530 fCurrentCoefDetector[k] = 0.0;
2532 fVectorFit.SetName("gainfactorscoefficients");
2534 // fDebug == 0 nothing
2535 // fDebug == 1 and fFitVoir no histo
2536 if (fDebugLevel == 1) {
2537 if(!CheckFitVoir()) return kFALSE;
2539 //Get the CalDet object
2541 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2543 AliInfo("Could not get calibDB");
2546 if(fCalDet) delete fCalDet;
2547 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2550 Float_t devalue = 1.0;
2551 if(fCalDet) delete fCalDet;
2552 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2553 for(Int_t k = 0; k < 540; k++){
2554 fCalDet->SetValue(k,devalue);
2560 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2561 Bool_t AliTRDCalibraFit::InitFitPH()
2564 // Init the arrays of results
2565 // Init the histos for debugging
2569 fVectorFit.SetName("driftvelocitycoefficients");
2570 fVectorFit2.SetName("t0coefficients");
2572 fCurrentCoefDetector = new Float_t[2304];
2573 for (Int_t k = 0; k < 2304; k++) {
2574 fCurrentCoefDetector[k] = 0.0;
2577 fCurrentCoefDetector2 = new Float_t[2304];
2578 for (Int_t k = 0; k < 2304; k++) {
2579 fCurrentCoefDetector2[k] = 0.0;
2582 //fDebug == 0 nothing
2583 // fDebug == 1 and fFitVoir no histo
2584 if (fDebugLevel == 1) {
2585 if(!CheckFitVoir()) return kFALSE;
2587 //Get the CalDet object
2589 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2591 AliInfo("Could not get calibDB");
2594 if(fCalDet) delete fCalDet;
2595 if(fCalDet2) delete fCalDet2;
2596 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2597 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2600 Float_t devalue = 1.5;
2601 Float_t devalue2 = 0.0;
2602 if(fCalDet) delete fCalDet;
2603 if(fCalDet2) delete fCalDet2;
2604 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2605 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2606 for(Int_t k = 0; k < 540; k++){
2607 fCalDet->SetValue(k,devalue);
2608 fCalDet2->SetValue(k,devalue2);
2613 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2614 Bool_t AliTRDCalibraFit::InitFitPRF()
2617 // Init the calibration mode (Nz, Nrphi), the histograms for
2618 // debugging the fit methods if fDebug > 0,
2622 fVectorFit.SetName("prfwidthcoefficients");
2624 fCurrentCoefDetector = new Float_t[2304];
2625 for (Int_t k = 0; k < 2304; k++) {
2626 fCurrentCoefDetector[k] = 0.0;
2629 // fDebug == 0 nothing
2630 // fDebug == 1 and fFitVoir no histo
2631 if (fDebugLevel == 1) {
2632 if(!CheckFitVoir()) return kFALSE;
2636 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2637 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2640 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2645 fCurrentCoefDetector = new Float_t[2304];
2646 fCurrentCoefDetector2 = new Float_t[2304];
2647 for (Int_t k = 0; k < 2304; k++) {
2648 fCurrentCoefDetector[k] = 0.0;
2649 fCurrentCoefDetector2[k] = 0.0;
2652 //printf("test0\n");
2654 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2656 AliInfo("Could not get calibDB");
2660 //Get the CalDet object
2662 if(fCalDet) delete fCalDet;
2663 if(fCalDet2) delete fCalDet2;
2664 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2665 //printf("test1\n");
2666 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2667 //printf("test2\n");
2668 for(Int_t k = 0; k < 540; k++){
2669 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2671 //printf("test3\n");
2674 Float_t devalue = 1.5;
2675 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2676 if(fCalDet) delete fCalDet;
2677 if(fCalDet2) delete fCalDet2;
2678 //printf("test1\n");
2679 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2680 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2681 //printf("test2\n");
2682 for(Int_t k = 0; k < 540; k++){
2683 fCalDet->SetValue(k,devalue);
2684 fCalDet2->SetValue(k,devalue2);
2686 //printf("test3\n");
2691 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2692 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2695 // Init the current detector where we are fCountDet and the
2696 // next fCount for the functions Fit...
2699 // Loop on the Xbins of ch!!
2700 fCountDet = -1; // Current detector
2701 fCount = 0; // To find the next detector
2704 if (fDebugLevel >= 3) {
2705 // Set countdet to the detector
2706 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2707 // Set counter to write at the end of the detector
2709 // Get the right calib objects
2712 if(fDebugLevel == 1) {
2714 fCalibraMode->CalculXBins(fCountDet,i);
2715 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2716 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2718 fCalibraMode->CalculXBins(fCountDet,i);
2719 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2725 fCount = fCalibraMode->GetXbins(i);
2727 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2728 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2729 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2730 ,(Int_t) GetStack(fCountDet)
2731 ,(Int_t) GetSector(fCountDet),i);
2734 //_______________________________________________________________________________
2735 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2738 // Calculate the number of bins expected (calibration groups)
2741 fNumberOfBinsExpected = 0;
2743 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2744 fNumberOfBinsExpected = 1;
2748 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2749 fNumberOfBinsExpected = 18;
2753 fCalibraMode->ModePadCalibration(2,i);
2754 fCalibraMode->ModePadFragmentation(0,2,0,i);
2755 fCalibraMode->SetDetChamb2(i);
2756 if (fDebugLevel > 1) {
2757 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2759 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2760 fCalibraMode->ModePadCalibration(0,i);
2761 fCalibraMode->ModePadFragmentation(0,0,0,i);
2762 fCalibraMode->SetDetChamb0(i);
2763 if (fDebugLevel > 1) {
2764 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2766 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2769 //_______________________________________________________________________________
2770 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2773 // Calculate the range of fits
2778 if (fDebugLevel == 1) {
2782 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2784 fDect2 = fNumberOfBinsExpected;
2786 if (fDebugLevel >= 3) {
2787 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2788 fCalibraMode->CalculXBins(fCountDet,i);
2789 fDect1 = fCalibraMode->GetXbins(i);
2790 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2791 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2792 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2793 ,(Int_t) GetStack(fCountDet)
2794 ,(Int_t) GetSector(fCountDet),i);
2795 // Set for the next detector
2796 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2799 //_______________________________________________________________________________
2800 Bool_t AliTRDCalibraFit::CheckFitVoir()
2803 // Check if fFitVoir is in the range
2806 if (fFitVoir < fNumberOfBinsExpected) {
2807 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2810 AliInfo("fFitVoir is out of range of the histo!");
2815 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2816 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2819 // See if we are in a new detector and update the
2820 // variables fNfragZ and fNfragRphi if yes
2821 // Will never happen for only one detector (3 and 4)
2822 // Doesn't matter for 2
2824 if (fCount == idect) {
2825 // On en est au detector (or first detector in the group)
2827 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2828 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2829 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2830 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2831 ,(Int_t) GetStack(fCountDet)
2832 ,(Int_t) GetSector(fCountDet),i);
2833 // Set for the next detector
2834 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2839 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2840 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2843 // Reconstruct the min pad row, max pad row, min pad col and
2844 // max pad col of the calibration group for the Fit functions
2845 // idect is the calibration group inside the detector
2847 if (fDebugLevel != 1) {
2848 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2850 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2851 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2853 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2854 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2857 // For the case where there are not enough entries in the histograms
2858 // of the calibration group, the value present in the choosen database
2859 // will be put. A negativ sign enables to know that a fit was not possible.
2862 if (fDebugLevel == 1) {
2863 AliInfo("The element has not enough statistic to be fitted");
2865 else if (fNbDet > 0){
2866 Int_t firstdetector = fCountDet;
2867 Int_t lastdetector = fCountDet+fNbDet;
2868 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2869 ,idect,firstdetector,lastdetector));
2870 // loop over detectors
2871 for(Int_t det = firstdetector; det < lastdetector; det++){
2873 //Set the calibration object again
2877 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2879 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2880 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2881 ,(Int_t) GetStack(fCountDet)
2882 ,(Int_t) GetSector(fCountDet),0);
2883 // Reconstruct row min row max
2884 ReconstructFitRowMinRowMax(idect,0);
2886 // Calcul the coef from the database choosen for the detector
2887 CalculChargeCoefMean(kFALSE);
2889 //stack 2, not stack 2
2891 if(GetStack(fCountDet) == 2) factor = 12;
2894 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2895 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2896 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2897 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2901 //Put default value negative
2902 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2903 fCurrentCoefE = 0.0;
2908 if(fDebugLevel > 1){
2910 if ( !fDebugStreamer ) {
2912 TDirectory *backup = gDirectory;
2913 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2914 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2917 Int_t detector = fCountDet;
2918 Int_t caligroup = idect;
2919 Short_t rowmin = fCalibraMode->GetRowMin(0);
2920 Short_t rowmax = fCalibraMode->GetRowMax(0);
2921 Short_t colmin = fCalibraMode->GetColMin(0);
2922 Short_t colmax = fCalibraMode->GetColMax(0);
2923 Float_t gf = fCurrentCoef[0];
2924 Float_t gfs = fCurrentCoef[1];
2925 Float_t gfE = fCurrentCoefE;
2927 (*fDebugStreamer) << "FillFillCH" <<
2928 "detector=" << detector <<
2929 "caligroup=" << caligroup <<
2930 "rowmin=" << rowmin <<
2931 "rowmax=" << rowmax <<
2932 "colmin=" << colmin <<
2933 "colmax=" << colmax <<
2941 for (Int_t k = 0; k < 2304; k++) {
2942 fCurrentCoefDetector[k] = 0.0;
2946 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2950 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2951 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2953 // Calcul the coef from the database choosen
2954 CalculChargeCoefMean(kFALSE);
2956 //stack 2, not stack 2
2958 if(GetStack(fCountDet) == 2) factor = 12;
2961 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2962 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2963 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2964 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2968 //Put default value negative
2969 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2970 fCurrentCoefE = 0.0;
2979 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2980 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2983 // For the case where there are not enough entries in the histograms
2984 // of the calibration group, the value present in the choosen database
2985 // will be put. A negativ sign enables to know that a fit was not possible.
2987 if (fDebugLevel == 1) {
2988 AliInfo("The element has not enough statistic to be fitted");
2990 else if (fNbDet > 0) {
2992 Int_t firstdetector = fCountDet;
2993 Int_t lastdetector = fCountDet+fNbDet;
2994 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2995 ,idect,firstdetector,lastdetector));
2996 // loop over detectors
2997 for(Int_t det = firstdetector; det < lastdetector; det++){
2999 //Set the calibration object again
3003 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3005 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3006 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3007 ,(Int_t) GetStack(fCountDet)
3008 ,(Int_t) GetSector(fCountDet),1);
3009 // Reconstruct row min row max
3010 ReconstructFitRowMinRowMax(idect,1);
3012 // Calcul the coef from the database choosen for the detector
3013 CalculVdriftCoefMean();
3016 //stack 2, not stack 2
3018 if(GetStack(fCountDet) == 2) factor = 12;
3021 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3022 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3023 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3024 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3025 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3029 //Put default value negative
3030 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3031 fCurrentCoefE = 0.0;
3032 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3033 fCurrentCoefE2 = 0.0;
3039 if(fDebugLevel > 1){
3041 if ( !fDebugStreamer ) {
3043 TDirectory *backup = gDirectory;
3044 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3045 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3049 Int_t detector = fCountDet;
3050 Int_t caligroup = idect;
3051 Short_t rowmin = fCalibraMode->GetRowMin(1);
3052 Short_t rowmax = fCalibraMode->GetRowMax(1);
3053 Short_t colmin = fCalibraMode->GetColMin(1);
3054 Short_t colmax = fCalibraMode->GetColMax(1);
3055 Float_t vf = fCurrentCoef[0];
3056 Float_t vs = fCurrentCoef[1];
3057 Float_t vfE = fCurrentCoefE;
3058 Float_t t0f = fCurrentCoef2[0];
3059 Float_t t0s = fCurrentCoef2[1];
3060 Float_t t0E = fCurrentCoefE2;
3064 (* fDebugStreamer) << "FillFillPH"<<
3065 "detector="<<detector<<
3066 "nentries="<<nentries<<
3067 "caligroup="<<caligroup<<
3081 for (Int_t k = 0; k < 2304; k++) {
3082 fCurrentCoefDetector[k] = 0.0;
3083 fCurrentCoefDetector2[k] = 0.0;
3087 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3091 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3092 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3094 CalculVdriftCoefMean();
3097 //stack 2 and not stack 2
3099 if(GetStack(fCountDet) == 2) factor = 12;
3103 // Fill the fCurrentCoefDetector 2
3104 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3105 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3106 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3107 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3111 // Put the default value
3112 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3113 fCurrentCoefE = 0.0;
3114 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3115 fCurrentCoefE2 = 0.0;
3117 FillFillPH(idect,nentries);
3126 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3127 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3130 // For the case where there are not enough entries in the histograms
3131 // of the calibration group, the value present in the choosen database
3132 // will be put. A negativ sign enables to know that a fit was not possible.
3135 if (fDebugLevel == 1) {
3136 AliInfo("The element has not enough statistic to be fitted");
3138 else if (fNbDet > 0){
3140 Int_t firstdetector = fCountDet;
3141 Int_t lastdetector = fCountDet+fNbDet;
3142 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3143 ,idect,firstdetector,lastdetector));
3145 // loop over detectors
3146 for(Int_t det = firstdetector; det < lastdetector; det++){
3148 //Set the calibration object again
3152 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3154 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3155 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3156 ,(Int_t) GetStack(fCountDet)
3157 ,(Int_t) GetSector(fCountDet),2);
3158 // Reconstruct row min row max
3159 ReconstructFitRowMinRowMax(idect,2);
3161 // Calcul the coef from the database choosen for the detector
3162 CalculPRFCoefMean();
3164 //stack 2, not stack 2
3166 if(GetStack(fCountDet) == 2) factor = 12;
3169 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3170 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3171 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3172 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3176 //Put default value negative
3177 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3178 fCurrentCoefE = 0.0;
3183 if(fDebugLevel > 1){
3185 if ( !fDebugStreamer ) {
3187 TDirectory *backup = gDirectory;
3188 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3189 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3192 Int_t detector = fCountDet;
3193 Int_t layer = GetLayer(fCountDet);
3194 Int_t caligroup = idect;
3195 Short_t rowmin = fCalibraMode->GetRowMin(2);
3196 Short_t rowmax = fCalibraMode->GetRowMax(2);
3197 Short_t colmin = fCalibraMode->GetColMin(2);
3198 Short_t colmax = fCalibraMode->GetColMax(2);
3199 Float_t widf = fCurrentCoef[0];
3200 Float_t wids = fCurrentCoef[1];
3201 Float_t widfE = fCurrentCoefE;
3203 (* fDebugStreamer) << "FillFillPRF"<<
3204 "detector="<<detector<<
3206 "caligroup="<<caligroup<<
3217 for (Int_t k = 0; k < 2304; k++) {
3218 fCurrentCoefDetector[k] = 0.0;
3222 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3226 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3227 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3229 CalculPRFCoefMean();
3231 // stack 2 and not stack 2
3233 if(GetStack(fCountDet) == 2) factor = 12;
3237 // Fill the fCurrentCoefDetector
3238 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3239 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3240 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3244 // Put the default value
3245 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3246 fCurrentCoefE = 0.0;
3254 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3255 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3258 // For the case where there are not enough entries in the histograms
3259 // of the calibration group, the value present in the choosen database
3260 // will be put. A negativ sign enables to know that a fit was not possible.
3263 // Calcul the coef from the database choosen
3264 CalculVdriftLorentzCoef();
3267 if(GetStack(fCountDet) == 2) factor = 1728;
3271 // Fill the fCurrentCoefDetector
3272 for (Int_t k = 0; k < factor; k++) {
3273 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3274 // should be negative
3275 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3279 //Put default opposite sign
3280 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3281 fCurrentCoefE = 0.0;
3282 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3283 fCurrentCoefE2 = 0.0;
3285 FillFillLinearFitter();
3290 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3291 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3294 // Fill the coefficients found with the fits or other
3295 // methods from the Fit functions
3298 if (fDebugLevel != 1) {
3300 Int_t firstdetector = fCountDet;
3301 Int_t lastdetector = fCountDet+fNbDet;
3302 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3303 ,idect,firstdetector,lastdetector));
3304 // loop over detectors
3305 for(Int_t det = firstdetector; det < lastdetector; det++){
3307 //Set the calibration object again
3311 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3313 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3314 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3315 ,(Int_t) GetStack(fCountDet)
3316 ,(Int_t) GetSector(fCountDet),0);
3317 // Reconstruct row min row max
3318 ReconstructFitRowMinRowMax(idect,0);
3320 // Calcul the coef from the database choosen for the detector
3321 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3322 else CalculChargeCoefMean(kTRUE);
3324 //stack 2, not stack 2
3326 if(GetStack(fCountDet) == 2) factor = 12;
3329 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3330 Double_t coeftoput = 1.0;
3331 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3332 else coeftoput = fCurrentCoef[0];
3333 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3334 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3335 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3342 if(fDebugLevel > 1){
3344 if ( !fDebugStreamer ) {
3346 TDirectory *backup = gDirectory;
3347 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3348 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3351 Int_t detector = fCountDet;
3352 Int_t caligroup = idect;
3353 Short_t rowmin = fCalibraMode->GetRowMin(0);
3354 Short_t rowmax = fCalibraMode->GetRowMax(0);
3355 Short_t colmin = fCalibraMode->GetColMin(0);
3356 Short_t colmax = fCalibraMode->GetColMax(0);
3357 Float_t gf = fCurrentCoef[0];
3358 Float_t gfs = fCurrentCoef[1];
3359 Float_t gfE = fCurrentCoefE;
3361 (*fDebugStreamer) << "FillFillCH" <<
3362 "detector=" << detector <<
3363 "caligroup=" << caligroup <<
3364 "rowmin=" << rowmin <<
3365 "rowmax=" << rowmax <<
3366 "colmin=" << colmin <<
3367 "colmax=" << colmax <<
3375 for (Int_t k = 0; k < 2304; k++) {
3376 fCurrentCoefDetector[k] = 0.0;
3380 //printf("Check the count now: fCountDet %d\n",fCountDet);
3385 if(GetStack(fCountDet) == 2) factor = 12;
3388 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3389 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3390 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3401 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3402 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3405 // Fill the coefficients found with the fits or other
3406 // methods from the Fit functions
3409 if (fDebugLevel != 1) {
3412 Int_t firstdetector = fCountDet;
3413 Int_t lastdetector = fCountDet+fNbDet;
3414 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3415 ,idect,firstdetector,lastdetector));
3417 // loop over detectors
3418 for(Int_t det = firstdetector; det < lastdetector; det++){
3420 //Set the calibration object again
3424 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3426 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3427 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3428 ,(Int_t) GetStack(fCountDet)
3429 ,(Int_t) GetSector(fCountDet),1);
3430 // Reconstruct row min row max
3431 ReconstructFitRowMinRowMax(idect,1);
3433 // Calcul the coef from the database choosen for the detector
3434 CalculVdriftCoefMean();
3437 //stack 2, not stack 2
3439 if(GetStack(fCountDet) == 2) factor = 12;
3442 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3443 Double_t coeftoput = 1.5;
3444 Double_t coeftoput2 = 0.0;
3446 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3447 else coeftoput = fCurrentCoef[0];
3449 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3450 else coeftoput2 = fCurrentCoef2[0];
3452 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3453 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3454 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3455 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3463 if(fDebugLevel > 1){
3465 if ( !fDebugStreamer ) {
3467 TDirectory *backup = gDirectory;
3468 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3469 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3473 Int_t detector = fCountDet;
3474 Int_t caligroup = idect;
3475 Short_t rowmin = fCalibraMode->GetRowMin(1);
3476 Short_t rowmax = fCalibraMode->GetRowMax(1);
3477 Short_t colmin = fCalibraMode->GetColMin(1);
3478 Short_t colmax = fCalibraMode->GetColMax(1);
3479 Float_t vf = fCurrentCoef[0];
3480 Float_t vs = fCurrentCoef[1];
3481 Float_t vfE = fCurrentCoefE;
3482 Float_t t0f = fCurrentCoef2[0];
3483 Float_t t0s = fCurrentCoef2[1];
3484 Float_t t0E = fCurrentCoefE2;
3488 (* fDebugStreamer) << "FillFillPH"<<
3489 "detector="<<detector<<
3490 "nentries="<<nentries<<
3491 "caligroup="<<caligroup<<
3505 for (Int_t k = 0; k < 2304; k++) {
3506 fCurrentCoefDetector[k] = 0.0;
3507 fCurrentCoefDetector2[k] = 0.0;
3511 //printf("Check the count now: fCountDet %d\n",fCountDet);
3516 if(GetStack(fCountDet) == 2) factor = 12;
3519 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3520 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3521 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3522 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3526 FillFillPH(idect,nentries);
3531 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3532 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3535 // Fill the coefficients found with the fits or other
3536 // methods from the Fit functions
3539 if (fDebugLevel != 1) {
3542 Int_t firstdetector = fCountDet;
3543 Int_t lastdetector = fCountDet+fNbDet;
3544 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3545 ,idect,firstdetector,lastdetector));
3547 // loop over detectors
3548 for(Int_t det = firstdetector; det < lastdetector; det++){
3550 //Set the calibration object again
3554 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3556 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3557 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3558 ,(Int_t) GetStack(fCountDet)
3559 ,(Int_t) GetSector(fCountDet),2);
3560 // Reconstruct row min row max
3561 ReconstructFitRowMinRowMax(idect,2);
3563 // Calcul the coef from the database choosen for the detector
3564 CalculPRFCoefMean();
3566 //stack 2, not stack 2
3568 if(GetStack(fCountDet) == 2) factor = 12;
3571 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3572 Double_t coeftoput = 1.0;
3573 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3574 else coeftoput = fCurrentCoef[0];
3575 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3576 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3577 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3584 if(fDebugLevel > 1){
3586 if ( !fDebugStreamer ) {
3588 TDirectory *backup = gDirectory;
3589 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3590 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3593 Int_t detector = fCountDet;
3594 Int_t layer = GetLayer(fCountDet);
3595 Int_t caligroup = idect;
3596 Short_t rowmin = fCalibraMode->GetRowMin(2);
3597 Short_t rowmax = fCalibraMode->GetRowMax(2);
3598 Short_t colmin = fCalibraMode->GetColMin(2);
3599 Short_t colmax = fCalibraMode->GetColMax(2);
3600 Float_t widf = fCurrentCoef[0];
3601 Float_t wids = fCurrentCoef[1];
3602 Float_t widfE = fCurrentCoefE;
3604 (* fDebugStreamer) << "FillFillPRF"<<
3605 "detector="<<detector<<
3607 "caligroup="<<caligroup<<
3618 for (Int_t k = 0; k < 2304; k++) {
3619 fCurrentCoefDetector[k] = 0.0;
3623 //printf("Check the count now: fCountDet %d\n",fCountDet);
3628 if(GetStack(fCountDet) == 2) factor = 12;
3631 // Pointer to the branch
3632 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3633 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3634 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3644 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3645 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3648 // Fill the coefficients found with the fits or other
3649 // methods from the Fit functions
3653 if(GetStack(fCountDet) == 2) factor = 1728;
3656 // Pointer to the branch
3657 for (Int_t k = 0; k < factor; k++) {
3658 fCurrentCoefDetector[k] = fCurrentCoef[0];
3659 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3662 FillFillLinearFitter();
3667 //________________________________________________________________________________
3668 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3671 // DebugStream and fVectorFit
3674 // End of one detector
3675 if ((idect == (fCount-1))) {
3678 for (Int_t k = 0; k < 2304; k++) {
3679 fCurrentCoefDetector[k] = 0.0;
3683 if(fDebugLevel > 1){
3685 if ( !fDebugStreamer ) {
3687 TDirectory *backup = gDirectory;
3688 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3689 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3692 Int_t detector = fCountDet;
3693 Int_t caligroup = idect;
3694 Short_t rowmin = fCalibraMode->GetRowMin(0);
3695 Short_t rowmax = fCalibraMode->GetRowMax(0);
3696 Short_t colmin = fCalibraMode->GetColMin(0);
3697 Short_t colmax = fCalibraMode->GetColMax(0);
3698 Float_t gf = fCurrentCoef[0];
3699 Float_t gfs = fCurrentCoef[1];
3700 Float_t gfE = fCurrentCoefE;
3702 (*fDebugStreamer) << "FillFillCH" <<
3703 "detector=" << detector <<
3704 "caligroup=" << caligroup <<
3705 "rowmin=" << rowmin <<
3706 "rowmax=" << rowmax <<
3707 "colmin=" << colmin <<
3708 "colmax=" << colmax <<
3716 //________________________________________________________________________________
3717 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3720 // DebugStream and fVectorFit and fVectorFit2
3723 // End of one detector
3724 if ((idect == (fCount-1))) {
3728 for (Int_t k = 0; k < 2304; k++) {
3729 fCurrentCoefDetector[k] = 0.0;
3730 fCurrentCoefDetector2[k] = 0.0;
3734 if(fDebugLevel > 1){
3736 if ( !fDebugStreamer ) {
3738 TDirectory *backup = gDirectory;
3739 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3740 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3744 Int_t detector = fCountDet;
3745 Int_t caligroup = idect;
3746 Short_t rowmin = fCalibraMode->GetRowMin(1);
3747 Short_t rowmax = fCalibraMode->GetRowMax(1);
3748 Short_t colmin = fCalibraMode->GetColMin(1);
3749 Short_t colmax = fCalibraMode->GetColMax(1);
3750 Float_t vf = fCurrentCoef[0];
3751 Float_t vs = fCurrentCoef[1];
3752 Float_t vfE = fCurrentCoefE;
3753 Float_t t0f = fCurrentCoef2[0];
3754 Float_t t0s = fCurrentCoef2[1];
3755 Float_t t0E = fCurrentCoefE2;
3759 (* fDebugStreamer) << "FillFillPH"<<
3760 "detector="<<detector<<
3761 "nentries="<<nentries<<
3762 "caligroup="<<caligroup<<
3777 //________________________________________________________________________________
3778 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3781 // DebugStream and fVectorFit
3784 // End of one detector
3785 if ((idect == (fCount-1))) {
3788 for (Int_t k = 0; k < 2304; k++) {
3789 fCurrentCoefDetector[k] = 0.0;
3794 if(fDebugLevel > 1){
3796 if ( !fDebugStreamer ) {
3798 TDirectory *backup = gDirectory;
3799 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3800 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3803 Int_t detector = fCountDet;
3804 Int_t layer = GetLayer(fCountDet);
3805 Int_t caligroup = idect;
3806 Short_t rowmin = fCalibraMode->GetRowMin(2);
3807 Short_t rowmax = fCalibraMode->GetRowMax(2);
3808 Short_t colmin = fCalibraMode->GetColMin(2);
3809 Short_t colmax = fCalibraMode->GetColMax(2);
3810 Float_t widf = fCurrentCoef[0];
3811 Float_t wids = fCurrentCoef[1];
3812 Float_t widfE = fCurrentCoefE;
3814 (* fDebugStreamer) << "FillFillPRF"<<
3815 "detector="<<detector<<
3817 "caligroup="<<caligroup<<
3829 //________________________________________________________________________________
3830 void AliTRDCalibraFit::FillFillLinearFitter()
3833 // DebugStream and fVectorFit
3836 // End of one detector
3842 for (Int_t k = 0; k < 2304; k++) {
3843 fCurrentCoefDetector[k] = 0.0;
3844 fCurrentCoefDetector2[k] = 0.0;
3848 if(fDebugLevel > 1){
3850 if ( !fDebugStreamer ) {
3852 TDirectory *backup = gDirectory;
3853 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3854 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3857 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3858 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3859 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3860 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3861 Float_t tiltangle = padplane->GetTiltingAngle();
3862 Int_t detector = fCountDet;
3863 Int_t stack = GetStack(fCountDet);
3864 Int_t layer = GetLayer(fCountDet);
3865 Float_t vf = fCurrentCoef[0];
3866 Float_t vs = fCurrentCoef[1];
3867 Float_t vfE = fCurrentCoefE;
3868 Float_t lorentzangler = fCurrentCoef2[0];
3869 Float_t elorentzangler = fCurrentCoefE2;
3870 Float_t lorentzangles = fCurrentCoef2[1];
3872 (* fDebugStreamer) << "FillFillLinearFitter"<<
3873 "detector="<<detector<<
3878 "tiltangle="<<tiltangle<<
3882 "lorentzangler="<<lorentzangler<<
3883 "Elorentzangler="<<elorentzangler<<
3884 "lorentzangles="<<lorentzangles<<
3890 //____________Calcul Coef Mean_________________________________________________
3892 //_____________________________________________________________________________
3893 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3896 // For the detector Dect calcul the mean time 0
3897 // for the calibration group idect from the choosen database
3900 fCurrentCoef2[1] = 0.0;
3901 if(fDebugLevel != 1){
3902 if(((fCalibraMode->GetNz(1) > 0) ||
3903 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3905 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3906 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3907 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3911 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3917 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3921 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3922 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3923 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3926 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3934 //_____________________________________________________________________________
3935 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3938 // For the detector Dect calcul the mean gain factor
3939 // for the calibration group idect from the choosen database
3942 fCurrentCoef[1] = 0.0;
3943 if(fDebugLevel != 1){
3944 if (((fCalibraMode->GetNz(0) > 0) ||
3945 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3946 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3947 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3948 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3949 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3952 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3956 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3957 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3962 //_____________________________________________________________________________
3963 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3966 // For the detector Dect calcul the mean sigma of pad response
3967 // function for the calibration group idect from the choosen database
3970 fCurrentCoef[1] = 0.0;
3971 if(fDebugLevel != 1){
3972 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3973 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3974 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3977 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3981 //_____________________________________________________________________________
3982 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3985 // For the detector dect calcul the mean drift velocity for the
3986 // calibration group idect from the choosen database
3989 fCurrentCoef[1] = 0.0;
3990 if(fDebugLevel != 1){
3991 if (((fCalibraMode->GetNz(1) > 0) ||
3992 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3994 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3995 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3996 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4000 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4005 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4010 //_____________________________________________________________________________
4011 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4014 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4017 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
4018 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4022 //_____________________________________________________________________________
4023 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4026 // Default width of the PRF if there is no database as reference
4031 //case 0: return 0.515;
4032 //case 1: return 0.502;
4033 //case 2: return 0.491;
4034 //case 3: return 0.481;
4035 //case 4: return 0.471;
4036 //case 5: return 0.463;
4037 //default: return 0.0;
4040 case 0: return 0.538429;
4041 case 1: return 0.524302;
4042 case 2: return 0.511591;
4043 case 3: return 0.500140;
4044 case 4: return 0.489821;
4045 case 5: return 0.480524;
4046 default: return 0.0;
4049 //________________________________________________________________________________
4050 void AliTRDCalibraFit::SetCalROC(Int_t i)
4053 // Set the calib object for fCountDet
4056 Float_t value = 0.0;
4058 //Get the CalDet object
4060 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4062 AliInfo("Could not get calibDB");
4069 fCalROC->~AliTRDCalROC();
4070 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4071 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4075 fCalROC->~AliTRDCalROC();
4076 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4077 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4079 fCalROC2->~AliTRDCalROC();
4080 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4081 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4085 fCalROC->~AliTRDCalROC();
4086 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4087 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4096 if(fCalROC) delete fCalROC;
4097 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4098 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4099 fCalROC->SetValue(k,1.0);
4103 if(fCalROC) delete fCalROC;
4104 if(fCalROC2) delete fCalROC2;
4105 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4106 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4107 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4108 fCalROC->SetValue(k,1.0);
4109 fCalROC2->SetValue(k,0.0);
4113 if(fCalROC) delete fCalROC;
4114 value = GetPRFDefault(GetLayer(fCountDet));
4115 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4116 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4117 fCalROC->SetValue(k,value);
4125 //____________Fit Methods______________________________________________________
4127 //_____________________________________________________________________________
4128 void AliTRDCalibraFit::FitPente(TH1* projPH)
4131 // Slope methode for the drift velocity
4135 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4142 fCurrentCoefE = 0.0;
4143 fCurrentCoefE2 = 0.0;
4144 fCurrentCoef[0] = 0.0;
4145 fCurrentCoef2[0] = 0.0;
4146 TLine *line = new TLine();
4149 TAxis *xpph = projPH->GetXaxis();
4150 Int_t nbins = xpph->GetNbins();
4151 Double_t lowedge = xpph->GetBinLowEdge(1);
4152 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4153 Double_t widbins = (upedge - lowedge) / nbins;
4154 Double_t limit = upedge + 0.5 * widbins;
4157 // Beginning of the signal
4158 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4159 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4160 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4162 binmax = (Int_t) pentea->GetMaximumBin();
4165 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4167 if (binmax >= nbins) {
4170 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4172 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4173 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4174 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4175 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4176 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4177 if (TMath::Abs(l3P2am) > 0.00000001) {
4178 fPhd[0] = -(l3P1am / (2 * l3P2am));
4181 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4182 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4185 // Amplification region
4188 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4189 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4196 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4198 if (binmax >= nbins) {
4201 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4203 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4204 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4205 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4206 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4207 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4208 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4209 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4211 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4212 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4215 fCurrentCoefE2 = fCurrentCoefE;
4218 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4219 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4220 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4223 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4226 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4228 if (binmin >= nbins) {
4231 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4236 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4237 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4238 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4239 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4240 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4241 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4242 if (TMath::Abs(l3P2dr) > 0.00000001) {
4243 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4245 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4246 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4248 Float_t fPhdt0 = 0.0;
4249 Float_t t0Shift = 0.0;
4252 t0Shift = fT0Shift1;
4256 t0Shift = fT0Shift0;
4259 if ((fPhd[2] > fPhd[0]) &&
4260 (fPhd[2] > fPhd[1]) &&
4261 (fPhd[1] > fPhd[0]) &&
4263 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4264 fNumberFitSuccess++;
4266 if (fPhdt0 >= 0.0) {
4267 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4268 if (fCurrentCoef2[0] < -1.0) {
4269 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4273 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4278 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4279 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4282 if (fDebugLevel == 1) {
4283 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4286 line->SetLineColor(2);
4287 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4288 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4289 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4290 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4291 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4292 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4293 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4294 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4297 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4306 //_____________________________________________________________________________
4307 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4310 // Slope methode but with polynomes de Lagrange
4314 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4317 //Double_t *x = new Double_t[5];
4318 //Double_t *y = new Double_t[5];
4335 fCurrentCoefE = 0.0;
4336 fCurrentCoefE2 = 1.0;
4337 fCurrentCoef[0] = 0.0;
4338 fCurrentCoef2[0] = 0.0;
4339 TLine *line = new TLine();
4340 TF1 * polynome = 0x0;
4341 TF1 * polynomea = 0x0;
4342 TF1 * polynomeb = 0x0;
4350 TAxis *xpph = projPH->GetXaxis();
4351 Int_t nbins = xpph->GetNbins();
4352 Double_t lowedge = xpph->GetBinLowEdge(1);
4353 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4354 Double_t widbins = (upedge - lowedge) / nbins;
4355 Double_t limit = upedge + 0.5 * widbins;
4360 // Beginning of the signal
4361 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4362 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4363 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4366 binmax = (Int_t) pentea->GetMaximumBin();
4368 Double_t minnn = 0.0;
4369 Double_t maxxx = 0.0;
4371 Int_t kase = nbins-binmax;
4379 minnn = pentea->GetBinCenter(binmax-2);
4380 maxxx = pentea->GetBinCenter(binmax);
4381 x[0] = pentea->GetBinCenter(binmax-2);
4382 x[1] = pentea->GetBinCenter(binmax-1);
4383 x[2] = pentea->GetBinCenter(binmax);
4384 y[0] = pentea->GetBinContent(binmax-2);
4385 y[1] = pentea->GetBinContent(binmax-1);
4386 y[2] = pentea->GetBinContent(binmax);
4387 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4388 AliInfo("At the limit for beginning!");
4391 minnn = pentea->GetBinCenter(binmax-2);
4392 maxxx = pentea->GetBinCenter(binmax+1);
4393 x[0] = pentea->GetBinCenter(binmax-2);
4394 x[1] = pentea->GetBinCenter(binmax-1);
4395 x[2] = pentea->GetBinCenter(binmax);
4396 x[3] = pentea->GetBinCenter(binmax+1);
4397 y[0] = pentea->GetBinContent(binmax-2);
4398 y[1] = pentea->GetBinContent(binmax-1);
4399 y[2] = pentea->GetBinContent(binmax);
4400 y[3] = pentea->GetBinContent(binmax+1);
4401 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4409 minnn = pentea->GetBinCenter(binmax);
4410 maxxx = pentea->GetBinCenter(binmax+2);
4411 x[0] = pentea->GetBinCenter(binmax);
4412 x[1] = pentea->GetBinCenter(binmax+1);
4413 x[2] = pentea->GetBinCenter(binmax+2);
4414 y[0] = pentea->GetBinContent(binmax);
4415 y[1] = pentea->GetBinContent(binmax+1);
4416 y[2] = pentea->GetBinContent(binmax+2);
4417 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4420 minnn = pentea->GetBinCenter(binmax-1);
4421 maxxx = pentea->GetBinCenter(binmax+2);
4422 x[0] = pentea->GetBinCenter(binmax-1);
4423 x[1] = pentea->GetBinCenter(binmax);
4424 x[2] = pentea->GetBinCenter(binmax+1);
4425 x[3] = pentea->GetBinCenter(binmax+2);
4426 y[0] = pentea->GetBinContent(binmax-1);
4427 y[1] = pentea->GetBinContent(binmax);
4428 y[2] = pentea->GetBinContent(binmax+1);
4429 y[3] = pentea->GetBinContent(binmax+2);
4430 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4433 minnn = pentea->GetBinCenter(binmax-2);
4434 maxxx = pentea->GetBinCenter(binmax+2);
4435 x[0] = pentea->GetBinCenter(binmax-2);
4436 x[1] = pentea->GetBinCenter(binmax-1);
4437 x[2] = pentea->GetBinCenter(binmax);
4438 x[3] = pentea->GetBinCenter(binmax+1);
4439 x[4] = pentea->GetBinCenter(binmax+2);
4440 y[0] = pentea->GetBinContent(binmax-2);
4441 y[1] = pentea->GetBinContent(binmax-1);
4442 y[2] = pentea->GetBinContent(binmax);
4443 y[3] = pentea->GetBinContent(binmax+1);
4444 y[4] = pentea->GetBinContent(binmax+2);
4445 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4453 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4454 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4456 Double_t step = (maxxx-minnn)/10000;
4458 Double_t maxvalue = 0.0;
4459 Double_t placemaximum = minnn;
4460 for(Int_t o = 0; o < 10000; o++){
4461 if(o == 0) maxvalue = polynomeb->Eval(l);
4462 if(maxvalue < (polynomeb->Eval(l))){
4463 maxvalue = polynomeb->Eval(l);
4468 fPhd[0] = placemaximum;
4471 // Amplification region
4474 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4475 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4481 Double_t minn = 0.0;
4482 Double_t maxx = 0.0;
4494 Int_t kase1 = nbins - binmax;
4496 //Determination of minn and maxx
4497 //case binmax = nbins
4502 minn = projPH->GetBinCenter(binmax-2);
4503 maxx = projPH->GetBinCenter(binmax);
4504 x[0] = projPH->GetBinCenter(binmax-2);
4505 x[1] = projPH->GetBinCenter(binmax-1);
4506 x[2] = projPH->GetBinCenter(binmax);
4507 y[0] = projPH->GetBinContent(binmax-2);
4508 y[1] = projPH->GetBinContent(binmax-1);
4509 y[2] = projPH->GetBinContent(binmax);
4510 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4511 //AliInfo("At the limit for the drift!");
4514 minn = projPH->GetBinCenter(binmax-2);
4515 maxx = projPH->GetBinCenter(binmax+1);
4516 x[0] = projPH->GetBinCenter(binmax-2);
4517 x[1] = projPH->GetBinCenter(binmax-1);
4518 x[2] = projPH->GetBinCenter(binmax);
4519 x[3] = projPH->GetBinCenter(binmax+1);
4520 y[0] = projPH->GetBinContent(binmax-2);
4521 y[1] = projPH->GetBinContent(binmax-1);
4522 y[2] = projPH->GetBinContent(binmax);
4523 y[3] = projPH->GetBinContent(binmax+1);
4524 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4533 minn = projPH->GetBinCenter(binmax);
4534 maxx = projPH->GetBinCenter(binmax+2);
4535 x[0] = projPH->GetBinCenter(binmax);
4536 x[1] = projPH->GetBinCenter(binmax+1);
4537 x[2] = projPH->GetBinCenter(binmax+2);
4538 y[0] = projPH->GetBinContent(binmax);
4539 y[1] = projPH->GetBinContent(binmax+1);
4540 y[2] = projPH->GetBinContent(binmax+2);
4541 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4544 minn = projPH->GetBinCenter(binmax-1);
4545 maxx = projPH->GetBinCenter(binmax+2);
4546 x[0] = projPH->GetBinCenter(binmax-1);
4547 x[1] = projPH->GetBinCenter(binmax);
4548 x[2] = projPH->GetBinCenter(binmax+1);
4549 x[3] = projPH->GetBinCenter(binmax+2);
4550 y[0] = projPH->GetBinContent(binmax-1);
4551 y[1] = projPH->GetBinContent(binmax);
4552 y[2] = projPH->GetBinContent(binmax+1);
4553 y[3] = projPH->GetBinContent(binmax+2);
4554 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4557 minn = projPH->GetBinCenter(binmax-2);
4558 maxx = projPH->GetBinCenter(binmax+2);
4559 x[0] = projPH->GetBinCenter(binmax-2);
4560 x[1] = projPH->GetBinCenter(binmax-1);
4561 x[2] = projPH->GetBinCenter(binmax);
4562 x[3] = projPH->GetBinCenter(binmax+1);
4563 x[4] = projPH->GetBinCenter(binmax+2);
4564 y[0] = projPH->GetBinContent(binmax-2);
4565 y[1] = projPH->GetBinContent(binmax-1);
4566 y[2] = projPH->GetBinContent(binmax);
4567 y[3] = projPH->GetBinContent(binmax+1);
4568 y[4] = projPH->GetBinContent(binmax+2);
4569 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4576 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4577 polynomea->SetParameters(c0,c1,c2,c3,c4);
4579 Double_t step = (maxx-minn)/1000;
4581 Double_t maxvalue = 0.0;
4582 Double_t placemaximum = minn;
4583 for(Int_t o = 0; o < 1000; o++){
4584 if(o == 0) maxvalue = polynomea->Eval(l);
4585 if(maxvalue < (polynomea->Eval(l))){
4586 maxvalue = polynomea->Eval(l);
4591 fPhd[1] = placemaximum;
4595 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4596 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4597 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4600 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4606 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4610 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4611 AliInfo("Too many fluctuations at the end!");
4614 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4615 AliInfo("Too many fluctuations at the end!");
4618 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4619 AliInfo("No entries for the next bin!");
4620 pente->SetBinContent(binmin,0);
4621 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4637 Bool_t case1 = kFALSE;
4638 Bool_t case2 = kFALSE;
4639 Bool_t case4 = kFALSE;
4641 //Determination of min and max
4642 //case binmin <= nbins-3
4644 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4645 min = pente->GetBinCenter(binmin-2);
4646 max = pente->GetBinCenter(binmin+2);
4647 x[0] = pente->GetBinCenter(binmin-2);
4648 x[1] = pente->GetBinCenter(binmin-1);
4649 x[2] = pente->GetBinCenter(binmin);
4650 x[3] = pente->GetBinCenter(binmin+1);
4651 x[4] = pente->GetBinCenter(binmin+2);
4652 y[0] = pente->GetBinContent(binmin-2);
4653 y[1] = pente->GetBinContent(binmin-1);
4654 y[2] = pente->GetBinContent(binmin);
4655 y[3] = pente->GetBinContent(binmin+1);
4656 y[4] = pente->GetBinContent(binmin+2);
4657 //Calcul the polynome de Lagrange
4658 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4660 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4661 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4662 //AliInfo("polynome 4 false 1");
4665 if(((binmin+3) <= (nbins-1)) &&
4666 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4667 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4668 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4669 AliInfo("polynome 4 false 2");
4673 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4674 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4675 //AliInfo("polynome 4 case 1");
4678 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4679 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4680 //AliInfo("polynome 4 case 4");
4685 //case binmin = nbins-2
4687 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4689 min = pente->GetBinCenter(binmin-2);
4690 max = pente->GetBinCenter(binmin+1);
4691 x[0] = pente->GetBinCenter(binmin-2);
4692 x[1] = pente->GetBinCenter(binmin-1);
4693 x[2] = pente->GetBinCenter(binmin);
4694 x[3] = pente->GetBinCenter(binmin+1);
4695 y[0] = pente->GetBinContent(binmin-2);
4696 y[1] = pente->GetBinContent(binmin-1);
4697 y[2] = pente->GetBinContent(binmin);
4698 y[3] = pente->GetBinContent(binmin+1);
4699 //Calcul the polynome de Lagrange
4700 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4701 //richtung +: nothing
4703 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4704 //AliInfo("polynome 3- case 2");
4709 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4711 min = pente->GetBinCenter(binmin-1);
4712 max = pente->GetBinCenter(binmin+2);
4713 x[0] = pente->GetBinCenter(binmin-1);
4714 x[1] = pente->GetBinCenter(binmin);
4715 x[2] = pente->GetBinCenter(binmin+1);
4716 x[3] = pente->GetBinCenter(binmin+2);
4717 y[0] = pente->GetBinContent(binmin-1);
4718 y[1] = pente->GetBinContent(binmin);
4719 y[2] = pente->GetBinContent(binmin+1);
4720 y[3] = pente->GetBinContent(binmin+2);
4721 //Calcul the polynome de Lagrange
4722 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4724 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4725 //AliInfo("polynome 3+ case 2");
4730 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4731 min = pente->GetBinCenter(binmin);
4732 max = pente->GetBinCenter(binmin+2);
4733 x[0] = pente->GetBinCenter(binmin);
4734 x[1] = pente->GetBinCenter(binmin+1);
4735 x[2] = pente->GetBinCenter(binmin+2);
4736 y[0] = pente->GetBinContent(binmin);
4737 y[1] = pente->GetBinContent(binmin+1);
4738 y[2] = pente->GetBinContent(binmin+2);
4739 //Calcul the polynome de Lagrange
4740 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4742 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4743 //AliInfo("polynome 2+ false");
4748 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4750 min = pente->GetBinCenter(binmin-1);
4751 max = pente->GetBinCenter(binmin+1);
4752 x[0] = pente->GetBinCenter(binmin-1);
4753 x[1] = pente->GetBinCenter(binmin);
4754 x[2] = pente->GetBinCenter(binmin+1);
4755 y[0] = pente->GetBinContent(binmin-1);
4756 y[1] = pente->GetBinContent(binmin);
4757 y[2] = pente->GetBinContent(binmin+1);
4758 //Calcul the polynome de Lagrange
4759 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4760 //richtung +: nothing
4761 //richtung -: nothing
4763 //case binmin = nbins-1
4765 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4766 min = pente->GetBinCenter(binmin-2);
4767 max = pente->GetBinCenter(binmin);
4768 x[0] = pente->GetBinCenter(binmin-2);
4769 x[1] = pente->GetBinCenter(binmin-1);
4770 x[2] = pente->GetBinCenter(binmin);
4771 y[0] = pente->GetBinContent(binmin-2);
4772 y[1] = pente->GetBinContent(binmin-1);
4773 y[2] = pente->GetBinContent(binmin);
4774 //Calcul the polynome de Lagrange
4775 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4776 //AliInfo("At the limit for the drift!");
4777 //fluctuation too big!
4778 //richtung +: nothing
4780 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4781 //AliInfo("polynome 2- false ");
4785 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4787 AliInfo("At the limit for the drift and not usable!");
4791 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4793 AliInfo("For the drift...problem!");
4795 //pass but should not happen
4796 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4798 AliInfo("For the drift...problem!");
4802 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4803 polynome->SetParameters(c0,c1,c2,c3,c4);
4804 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4805 Double_t step = (max-min)/1000;
4807 Double_t minvalue = 0.0;
4808 Double_t placeminimum = min;
4809 for(Int_t o = 0; o < 1000; o++){
4810 if(o == 0) minvalue = polynome->Eval(l);
4811 if(minvalue > (polynome->Eval(l))){
4812 minvalue = polynome->Eval(l);
4817 fPhd[2] = placeminimum;
4819 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4820 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4821 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4823 Float_t fPhdt0 = 0.0;
4824 Float_t t0Shift = 0.0;
4827 t0Shift = fT0Shift1;
4831 t0Shift = fT0Shift0;
4834 if ((fPhd[2] > fPhd[0]) &&
4835 (fPhd[2] > fPhd[1]) &&
4836 (fPhd[1] > fPhd[0]) &&
4838 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4839 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4840 else fNumberFitSuccess++;
4841 if (fPhdt0 >= 0.0) {
4842 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4843 if (fCurrentCoef2[0] < -1.0) {
4844 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4848 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4852 //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4853 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4855 if((fPhd[1] > fPhd[0]) &&
4857 if (fPhdt0 >= 0.0) {
4858 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4859 if (fCurrentCoef2[0] < -1.0) {
4860 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4864 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4868 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4869 //printf("Fit failed!\n");
4873 if (fDebugLevel == 1) {
4874 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4877 line->SetLineColor(2);
4878 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4879 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4880 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4881 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4882 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4883 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4884 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4885 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4888 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4895 if(polynome) delete polynome;
4896 if(polynomea) delete polynomea;
4897 if(polynomeb) delete polynomeb;
4898 //if(x) delete [] x;
4899 //if(y) delete [] y;
4900 if(line) delete line;
4905 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4906 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4908 projPH->SetDirectory(0);
4912 //_____________________________________________________________________________
4913 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4916 // Fit methode for the drift velocity
4920 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4923 TAxis *xpph = projPH->GetXaxis();
4924 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4926 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4927 fPH->SetParameter(0,0.469); // Scaling
4928 fPH->SetParameter(1,0.18); // Start
4929 fPH->SetParameter(2,0.0857325); // AR
4930 fPH->SetParameter(3,1.89); // DR
4931 fPH->SetParameter(4,0.08); // QA/QD
4932 fPH->SetParameter(5,0.0); // Baseline
4934 TLine *line = new TLine();
4936 fCurrentCoef[0] = 0.0;
4937 fCurrentCoef2[0] = 0.0;
4938 fCurrentCoefE = 0.0;
4939 fCurrentCoefE2 = 0.0;
4941 if (idect%fFitPHPeriode == 0) {
4943 AliInfo(Form("The detector %d will be fitted",idect));
4944 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4945 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4946 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4947 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4948 fPH->SetParameter(4,0.225); // QA/QD
4949 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4951 if (fDebugLevel != 1) {
4952 projPH->Fit(fPH,"0M","",0.0,upedge);
4955 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4957 projPH->Fit(fPH,"M+","",0.0,upedge);
4959 line->SetLineColor(4);
4960 line->DrawLine(fPH->GetParameter(1)
4962 ,fPH->GetParameter(1)
4963 ,projPH->GetMaximum());
4964 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4966 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4967 ,projPH->GetMaximum());
4968 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4970 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4971 ,projPH->GetMaximum());
4974 if (fPH->GetParameter(3) != 0) {
4975 fNumberFitSuccess++;
4976 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4977 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4978 fCurrentCoef2[0] = fPH->GetParameter(1);
4979 fCurrentCoefE2 = fPH->GetParError(1);
4982 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4983 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4989 // Put the default value
4990 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4991 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4994 if (fDebugLevel != 1) {
4999 //_____________________________________________________________________________
5000 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5003 // Fit methode for the sigma of the pad response function
5008 fCurrentCoef[0] = 0.0;
5009 fCurrentCoefE = 0.0;
5011 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5013 if(TMath::Abs(ret+4) <= 0.000000001){
5014 fCurrentCoef[0] = -fCurrentCoef[1];
5018 fNumberFitSuccess++;
5019 fCurrentCoef[0] = param[2];
5020 fCurrentCoefE = ret;
5024 //_____________________________________________________________________________
5025 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)
5028 // Fit methode for the sigma of the pad response function
5031 //We should have at least 3 points
5032 if(nBins <=3) return -4.0;
5034 TLinearFitter fitter(3,"pol2");
5035 fitter.StoreData(kFALSE);
5036 fitter.ClearPoints();
5038 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5039 Float_t entries = 0;
5040 Int_t nbbinwithentries = 0;
5041 for (Int_t i=0; i<nBins; i++){
5043 if(arraye[i] > 15) nbbinwithentries++;
5044 //printf("entries for i %d: %f\n",i,arraye[i]);
5046 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5047 //printf("entries %f\n",entries);
5048 //printf("nbbinwithentries %d\n",nbbinwithentries);
5051 Float_t errorm = 0.0;
5052 Float_t errorn = 0.0;
5053 Float_t error = 0.0;
5056 for (Int_t ibin=0;ibin<nBins; ibin++){
5057 Float_t entriesI = arraye[ibin];
5058 Float_t valueI = arraym[ibin];
5059 Double_t xcenter = 0.0;
5061 if ((entriesI>15) && (valueI>0.0)){
5062 xcenter = xMin+(ibin+0.5)*binWidth;
5067 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5068 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5069 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5072 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5073 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5074 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5076 if(TMath::Abs(error) < 0.000000001) continue;
5077 val = TMath::Log(Float_t(valueI));
5078 fitter.AddPoint(&xcenter,val,error);
5082 if(fDebugLevel > 1){
5084 if ( !fDebugStreamer ) {
5086 TDirectory *backup = gDirectory;
5087 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5088 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5091 Int_t detector = fCountDet;
5092 Int_t layer = GetLayer(fCountDet);
5095 (* fDebugStreamer) << "FitGausMIFill"<<
5096 "detector="<<detector<<
5100 "entriesI="<<entriesI<<
5103 "xcenter="<<xcenter<<
5113 if(npoints <=3) return -4.0;
5118 fitter.GetParameters(par);
5119 chi2 = fitter.GetChisquare()/Float_t(npoints);
5122 if (!param) param = new TVectorD(3);
5123 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5124 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5125 Double_t deltax = (fitter.GetParError(2))/x;
5126 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5129 (*param)[1] = par[1]/(-2.*par[2]);
5130 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5131 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5132 if ( lnparam0>307 ) return -4;
5133 (*param)[0] = TMath::Exp(lnparam0);
5135 if(fDebugLevel > 1){
5137 if ( !fDebugStreamer ) {
5139 TDirectory *backup = gDirectory;
5140 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5141 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5144 Int_t detector = fCountDet;
5145 Int_t layer = GetLayer(fCountDet);
5148 (* fDebugStreamer) << "FitGausMIFit"<<
5149 "detector="<<detector<<
5152 "errorsigma="<<chi2<<
5153 "mean="<<(*param)[1]<<
5154 "sigma="<<(*param)[2]<<
5155 "constant="<<(*param)[0]<<
5160 if((chi2/(*param)[2]) > 0.1){
5162 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5167 if(fDebugLevel == 1){
5168 TString name("PRF");
5169 name += (Int_t)xMin;
5170 name += (Int_t)xMax;
5171 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5174 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5175 for(Int_t k = 0; k < nBins; k++){
5176 histo->SetBinContent(k+1,arraym[k]);
5177 histo->SetBinError(k+1,arrayme[k]);
5180 name += "functionf";
5181 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5182 f1->SetParameter(0, (*param)[0]);
5183 f1->SetParameter(1, (*param)[1]);
5184 f1->SetParameter(2, (*param)[2]);
5192 //_____________________________________________________________________________
5193 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5196 // Fit methode for the sigma of the pad response function
5199 fCurrentCoef[0] = 0.0;
5200 fCurrentCoefE = 0.0;
5202 if (fDebugLevel != 1) {
5203 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5206 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5208 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5211 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5212 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5213 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5215 fNumberFitSuccess++;
5218 //_____________________________________________________________________________
5219 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5222 // Fit methode for the sigma of the pad response function
5224 fCurrentCoef[0] = 0.0;
5225 fCurrentCoefE = 0.0;
5226 if (fDebugLevel == 1) {
5227 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5231 fCurrentCoef[0] = projPRF->GetRMS();
5232 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5234 fNumberFitSuccess++;
5237 //_____________________________________________________________________________
5238 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5241 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5244 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5247 Int_t nbins = (Int_t)(nybins/(2*nbg));
5248 Float_t lowedge = -3.0*nbg;
5249 Float_t upedge = lowedge + 3.0;
5252 Double_t xvalues = -0.2*nbg+0.1;
5254 Int_t total = 2*nbg;
5257 for(Int_t k = 0; k < total; k++){
5258 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5260 y = fCurrentCoef[0]*fCurrentCoef[0];
5261 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5264 if(fDebugLevel > 1){
5266 if ( !fDebugStreamer ) {
5268 TDirectory *backup = gDirectory;
5269 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5270 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5273 Int_t detector = fCountDet;
5274 Int_t layer = GetLayer(fCountDet);
5275 Int_t nbtotal = total;
5277 Float_t low = lowedge;
5278 Float_t up = upedge;
5279 Float_t tnp = xvalues;
5280 Float_t wid = fCurrentCoef[0];
5281 Float_t widfE = fCurrentCoefE;
5283 (* fDebugStreamer) << "FitTnpRange0"<<
5284 "detector="<<detector<<
5286 "nbtotal="<<nbtotal<<
5304 fCurrentCoefE = 0.0;
5305 fCurrentCoef[0] = 0.0;
5307 //printf("npoints\n",npoints);
5310 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5315 linearfitter.Eval();
5316 linearfitter.GetParameters(pars0);
5317 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5318 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5319 Double_t min0 = 0.0;
5320 Double_t ermin0 = 0.0;
5321 //Double_t prfe0 = 0.0;
5322 Double_t prf0 = 0.0;
5323 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5324 min0 = -pars0[1]/(2*pars0[2]);
5325 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5326 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5329 prfe0 = linearfitter->GetParError(0)*pointError0
5330 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5331 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5332 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5333 fCurrentCoefE = (Float_t) prfe0;
5335 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5338 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5342 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5345 if(fDebugLevel > 1){
5347 if ( !fDebugStreamer ) {
5349 TDirectory *backup = gDirectory;
5350 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5351 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5354 Int_t detector = fCountDet;
5355 Int_t layer = GetLayer(fCountDet);
5356 Int_t nbtotal = total;
5357 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5358 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5360 (* fDebugStreamer) << "FitTnpRange1"<<
5361 "detector="<<detector<<
5363 "nbtotal="<<nbtotal<<
5367 "npoints="<<npoints<<
5370 "sigmaprf="<<fCurrentCoef[0]<<
5371 "sigprf="<<fCurrentCoef[1]<<
5378 //_____________________________________________________________________________
5379 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5382 // Only mean methode for the gain factor
5385 fCurrentCoef[0] = mean;
5386 fCurrentCoefE = 0.0;
5387 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5388 if (fDebugLevel == 1) {
5389 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5393 CalculChargeCoefMean(kTRUE);
5394 fNumberFitSuccess++;
5396 //_____________________________________________________________________________
5397 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5400 // mean w methode for the gain factor
5404 Int_t nybins = projch->GetNbinsX();
5406 //The weight function
5407 Double_t a = 0.00228515;
5408 Double_t b = -0.00231487;
5409 Double_t c = 0.00044298;
5410 Double_t d = -0.00379239;
5411 Double_t e = 0.00338349;
5421 //A arbitrary error for the moment
5422 fCurrentCoefE = 0.0;
5423 fCurrentCoef[0] = 0.0;
5426 Double_t sumw = 0.0;
5428 Float_t sumAll = (Float_t) nentries;
5429 Int_t sumCurrent = 0;
5430 for(Int_t k = 0; k <nybins; k++){
5431 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5432 if (fraction>0.95) break;
5433 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5434 e*fraction*fraction*fraction*fraction;
5435 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5436 sum += weight*projch->GetBinContent(k+1);
5437 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5438 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5440 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5442 if (fDebugLevel == 1) {
5443 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5447 fNumberFitSuccess++;
5448 CalculChargeCoefMean(kTRUE);
5450 //_____________________________________________________________________________
5451 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5454 // mean w methode for the gain factor
5458 Int_t nybins = projch->GetNbinsX();
5460 //The weight function
5461 Double_t a = 0.00228515;
5462 Double_t b = -0.00231487;
5463 Double_t c = 0.00044298;
5464 Double_t d = -0.00379239;
5465 Double_t e = 0.00338349;
5475 //A arbitrary error for the moment
5476 fCurrentCoefE = 0.0;
5477 fCurrentCoef[0] = 0.0;
5480 Double_t sumw = 0.0;
5482 Int_t sumCurrent = 0;
5483 for(Int_t k = 0; k <nybins; k++){
5484 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5485 if (fraction>0.95) break;
5486 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5487 e*fraction*fraction*fraction*fraction;
5488 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5489 sum += weight*projch->GetBinContent(k+1);
5490 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5491 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5493 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5495 if (fDebugLevel == 1) {
5496 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5500 fNumberFitSuccess++;
5502 //_____________________________________________________________________________
5503 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5506 // Fit methode for the gain factor
5509 fCurrentCoef[0] = 0.0;
5510 fCurrentCoefE = 0.0;
5511 Double_t chisqrl = 0.0;
5512 Double_t chisqrg = 0.0;
5513 Double_t chisqr = 0.0;
5514 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5516 projch->Fit("landau","0",""
5517 ,(Double_t) mean/fBeginFitCharge
5518 ,projch->GetBinCenter(projch->GetNbinsX()));
5519 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5520 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5521 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5522 chisqrl = projch->GetFunction("landau")->GetChisquare();
5524 projch->Fit("gaus","0",""
5525 ,(Double_t) mean/fBeginFitCharge
5526 ,projch->GetBinCenter(projch->GetNbinsX()));
5527 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5528 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5529 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5531 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5532 if (fDebugLevel != 1) {
5533 projch->Fit("fLandauGaus","0",""
5534 ,(Double_t) mean/fBeginFitCharge
5535 ,projch->GetBinCenter(projch->GetNbinsX()));
5536 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5539 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5541 projch->Fit("fLandauGaus","+",""
5542 ,(Double_t) mean/fBeginFitCharge
5543 ,projch->GetBinCenter(projch->GetNbinsX()));
5544 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5546 fLandauGaus->Draw("same");
5549 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5550 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5551 fNumberFitSuccess++;
5552 CalculChargeCoefMean(kTRUE);
5553 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5554 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5557 CalculChargeCoefMean(kFALSE);
5558 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5561 if (fDebugLevel != 1) {
5566 //_____________________________________________________________________________
5567 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5570 // Fit methode for the gain factor more time consuming
5574 //Some parameters to initialise
5575 Double_t widthLandau, widthGaus, mPV, integral;
5576 Double_t chisquarel = 0.0;
5577 Double_t chisquareg = 0.0;
5578 projch->Fit("landau","0M+",""
5580 ,projch->GetBinCenter(projch->GetNbinsX()));
5581 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5582 chisquarel = projch->GetFunction("landau")->GetChisquare();
5583 projch->Fit("gaus","0M+",""
5585 ,projch->GetBinCenter(projch->GetNbinsX()));
5586 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5587 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5589 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5590 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5592 // Setting fit range and start values
5594 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5595 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5596 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5597 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5598 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5599 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5600 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5603 fCurrentCoef[0] = 0.0;
5604 fCurrentCoefE = 0.0;
5608 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5613 Double_t projchPeak;
5614 Double_t projchFWHM;
5615 LanGauPro(fp,projchPeak,projchFWHM);
5617 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5618 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5619 fNumberFitSuccess++;
5620 CalculChargeCoefMean(kTRUE);
5621 fCurrentCoef[0] = fp[1];
5622 fCurrentCoefE = fpe[1];
5623 //chargeCoefE2 = chisqr;
5626 CalculChargeCoefMean(kFALSE);
5627 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5629 if (fDebugLevel == 1) {
5630 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5631 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5634 fitsnr->Draw("same");
5640 //_____________________________________________________________________________
5641 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
5644 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5646 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5647 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5648 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5653 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5654 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5658 //_____________________________________________________________________________
5659 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
5662 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5664 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5665 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5666 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5667 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5671 c2 = -(x0*(x[1]+x[2]+x[3])
5672 +x1*(x[0]+x[2]+x[3])
5673 +x2*(x[0]+x[1]+x[3])
5674 +x3*(x[0]+x[1]+x[2]));
5675 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5676 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5677 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5678 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5680 c0 = -(x0*x[1]*x[2]*x[3]
5683 +x3*x[0]*x[1]*x[2]);
5688 //_____________________________________________________________________________
5689 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
5692 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5694 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5695 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5696 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5697 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5698 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5701 c4 = x0+x1+x2+x3+x4;
5702 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5703 +x1*(x[0]+x[2]+x[3]+x[4])
5704 +x2*(x[0]+x[1]+x[3]+x[4])
5705 +x3*(x[0]+x[1]+x[2]+x[4])
5706 +x4*(x[0]+x[1]+x[2]+x[3]));
5707 c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5708 +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5709 +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
5710 +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
5711 +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
5713 c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5714 +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
5715 +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
5716 +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
5717 +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
5719 c0 = (x0*x[1]*x[2]*x[3]*x[4]
5720 +x1*x[0]*x[2]*x[3]*x[4]
5721 +x2*x[0]*x[1]*x[3]*x[4]
5722 +x3*x[0]*x[1]*x[2]*x[4]
5723 +x4*x[0]*x[1]*x[2]*x[3]);
5726 //_____________________________________________________________________________
5727 void AliTRDCalibraFit::NormierungCharge()
5730 // Normalisation of the gain factor resulting for the fits
5733 // Calcul of the mean of choosen method by fFitChargeNDB
5735 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5736 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5738 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5739 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5740 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5741 if (GetStack(detector) == 2) {
5744 if (GetStack(detector) != 2) {
5747 for (Int_t j = 0; j < total; j++) {
5755 fScaleFitFactor = fScaleFitFactor / sum;
5758 fScaleFitFactor = 1.0;
5761 //methode de boeuf mais bon...
5762 Double_t scalefactor = fScaleFitFactor;
5764 if(fDebugLevel > 1){
5766 if ( !fDebugStreamer ) {
5768 TDirectory *backup = gDirectory;
5769 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5770 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5772 (* fDebugStreamer) << "NormierungCharge"<<
5773 "scalefactor="<<scalefactor<<
5777 //_____________________________________________________________________________
5778 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5781 // Rebin of the 1D histo for the gain calibration if needed.
5782 // you have to choose fRebin, divider of fNumberBinCharge
5785 TAxis *xhist = hist->GetXaxis();
5786 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5787 ,xhist->GetBinLowEdge(1)
5788 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5790 AliInfo(Form("fRebin: %d",fRebin));
5792 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5794 for (Int_t ji = i; ji < i+fRebin; ji++) {
5795 sum += hist->GetBinContent(ji);
5798 rehist->SetBinContent(k,sum);
5806 //_____________________________________________________________________________
5807 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5810 // Rebin of the 1D histo for the gain calibration if needed
5811 // you have to choose fRebin divider of fNumberBinCharge
5814 TAxis *xhist = hist->GetXaxis();
5815 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5816 ,xhist->GetBinLowEdge(1)
5817 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5819 AliInfo(Form("fRebin: %d",fRebin));
5821 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5823 for (Int_t ji = i; ji < i+fRebin; ji++) {
5824 sum += hist->GetBinContent(ji);
5827 rehist->SetBinContent(k,sum);
5835 //____________Some basic geometry function_____________________________________
5838 //_____________________________________________________________________________
5839 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5842 // Reconstruct the plane number from the detector number
5845 return ((Int_t) (d % 6));
5849 //_____________________________________________________________________________
5850 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5853 // Reconstruct the stack number from the detector number
5855 const Int_t kNlayer = 6;
5857 return ((Int_t) (d % 30) / kNlayer);
5861 //_____________________________________________________________________________
5862 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5865 // Reconstruct the sector number from the detector number
5869 return ((Int_t) (d / fg));
5874 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5876 //_______________________________________________________________________________
5877 void AliTRDCalibraFit::ResetVectorFit()
5880 // Reset the VectorFits
5883 fVectorFit.SetOwner();
5885 fVectorFit2.SetOwner();
5886 fVectorFit2.Clear();
5890 //____________Private Functions________________________________________________
5893 //_____________________________________________________________________________
5894 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5897 // Function for the fit
5900 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5902 //PARAMETERS FOR FIT PH
5904 //fAsymmGauss->SetParameter(0,0.113755);
5905 //fAsymmGauss->SetParameter(1,0.350706);
5906 //fAsymmGauss->SetParameter(2,0.0604244);
5907 //fAsymmGauss->SetParameter(3,7.65596);
5908 //fAsymmGauss->SetParameter(4,1.00124);
5909 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5917 Double_t dx = 0.005;
5918 Double_t xs = par[1];
5920 Double_t paras[2] = { 0.0, 0.0 };
5923 if ((xs >= par[1]) &&
5924 (xs < (par[1]+par[2]))) {
5925 //fAsymmGauss->SetParameter(0,par[0]);
5926 //fAsymmGauss->SetParameter(1,xs);
5927 //ss += fAsymmGauss->Eval(xx);
5930 ss += AsymmGauss(&xx,paras);
5932 if ((xs >= (par[1]+par[2])) &&
5933 (xs < (par[1]+par[2]+par[3]))) {
5934 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5935 //fAsymmGauss->SetParameter(1,xs);
5936 //ss += fAsymmGauss->Eval(xx);
5937 paras[0] = par[0]*par[4];
5939 ss += AsymmGauss(&xx,paras);
5948 //_____________________________________________________________________________
5949 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5952 // Function for the fit
5955 //par[0] = normalization
5963 Double_t par1save = par[1];
5964 //Double_t par2save = par[2];
5965 Double_t par2save = 0.0604244;
5966 //Double_t par3save = par[3];
5967 Double_t par3save = 7.65596;
5968 //Double_t par5save = par[5];
5969 Double_t par5save = 0.870597;
5970 Double_t dx = x[0] - par1save;
5972 Double_t sigma2 = par2save*par2save;
5973 Double_t sqrt2 = TMath::Sqrt(2.0);
5974 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5975 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5976 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5977 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5979 //return par[0]*(exp1+par[4]*exp2);
5980 return par[0] * (exp1 + 1.00124 * exp2);
5984 //_____________________________________________________________________________
5985 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5988 // Sum Landau + Gaus with identical mean
5991 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5992 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5993 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5994 Double_t val = valLandau + valGaus;
6000 //_____________________________________________________________________________
6001 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6004 // Function for the fit
6007 // par[0]=Width (scale) parameter of Landau density
6008 // par[1]=Most Probable (MP, location) parameter of Landau density
6009 // par[2]=Total area (integral -inf to inf, normalization constant)
6010 // par[3]=Width (sigma) of convoluted Gaussian function
6012 // In the Landau distribution (represented by the CERNLIB approximation),
6013 // the maximum is located at x=-0.22278298 with the location parameter=0.
6014 // This shift is corrected within this function, so that the actual
6015 // maximum is identical to the MP parameter.
6018 // Numeric constants
6019 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6020 Double_t mpshift = -0.22278298; // Landau maximum location
6022 // Control constants
6023 Double_t np = 100.0; // Number of convolution steps
6024 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6036 // MP shift correction
6037 mpc = par[1] - mpshift * par[0];
6039 // Range of convolution integral
6040 xlow = x[0] - sc * par[3];
6041 xupp = x[0] + sc * par[3];
6043 step = (xupp - xlow) / np;
6045 // Convolution integral of Landau and Gaussian by sum
6046 for (i = 1.0; i <= np/2; i++) {
6048 xx = xlow + (i-.5) * step;
6049 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6050 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6052 xx = xupp - (i-.5) * step;
6053 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6054 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6058 return (par[2] * step * sum * invsq2pi / par[3]);
6061 //_____________________________________________________________________________
6062 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6063 , const Double_t *parlimitslo, const Double_t *parlimitshi
6064 , Double_t *fitparams, Double_t *fiterrors
6065 , Double_t *chiSqr, Int_t *ndf) const
6068 // Function for the fit
6072 Char_t funname[100];
6074 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6079 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6080 ffit->SetParameters(startvalues);
6081 ffit->SetParNames("Width","MP","Area","GSigma");
6083 for (i = 0; i < 4; i++) {
6084 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6087 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6089 ffit->GetParameters(fitparams); // Obtain fit parameters
6090 for (i = 0; i < 4; i++) {
6091 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6093 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6094 ndf[0] = ffit->GetNDF(); // Obtain ndf
6096 return (ffit); // Return fit function
6100 //_____________________________________________________________________________
6101 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6104 // Function for the fit
6117 Int_t maxcalls = 10000;
6119 // Search for maximum
6120 p = params[1] - 0.1 * params[0];
6121 step = 0.05 * params[0];
6125 while ((l != lold) && (i < maxcalls)) {
6129 l = LanGauFun(&x,params);
6131 step = -step / 10.0;
6136 if (i == maxcalls) {
6142 // Search for right x location of fy
6143 p = maxx + params[0];
6149 while ( (l != lold) && (i < maxcalls) ) {
6154 l = TMath::Abs(LanGauFun(&x,params) - fy);
6168 // Search for left x location of fy
6170 p = maxx - 0.5 * params[0];
6176 while ((l != lold) && (i < maxcalls)) {
6180 l = TMath::Abs(LanGauFun(&x,params) - fy);
6182 step = -step / 10.0;
6187 if (i == maxcalls) {
6196 //_____________________________________________________________________________
6197 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6200 // Gaus with identical mean
6203 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];