1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
63 #include <TTreeStream.h>
64 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
71 #include "AliTRDCalibraFit.h"
72 #include "AliTRDCalibraMode.h"
73 #include "AliTRDCalibraVector.h"
74 #include "AliTRDCalibraVdriftLinearFit.h"
75 #include "AliTRDcalibDB.h"
76 #include "AliTRDgeometry.h"
77 #include "AliTRDpadPlane.h"
78 #include "AliTRDgeometry.h"
79 #include "AliTRDCommonParam.h"
80 #include "./Cal/AliTRDCalROC.h"
81 #include "./Cal/AliTRDCalPad.h"
82 #include "./Cal/AliTRDCalDet.h"
85 ClassImp(AliTRDCalibraFit)
87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
90 //_____________singleton implementation_________________________________________________
91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
94 // Singleton implementation
97 if (fgTerminated != kFALSE) {
101 if (fgInstance == 0) {
102 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
124 //______________________________________________________________________________________
125 AliTRDCalibraFit::AliTRDCalibraFit()
128 ,fNumberOfBinsExpected(0)
130 ,fBeginFitCharge(3.5)
132 ,fTakeTheMaxPH(kTRUE)
140 ,fNumberFitSuccess(0)
147 ,fCalibraMode(new AliTRDCalibraMode())
152 ,fScaleFitFactor(0.0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
167 // Default constructor
170 fGeo = new AliTRDgeometry();
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
177 for (Int_t i = 0; i < 3; i++) {
183 //______________________________________________________________________________________
184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
189 ,fBeginFitCharge(c.fBeginFitCharge)
190 ,fFitPHPeriode(c.fFitPHPeriode)
191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
192 ,fT0Shift0(c.fT0Shift0)
193 ,fT0Shift1(c.fT0Shift1)
194 ,fRangeFitPRF(c.fRangeFitPRF)
196 ,fMinEntries(c.fMinEntries)
198 ,fNumberFit(c.fNumberFit)
199 ,fNumberFitSuccess(c.fNumberFitSuccess)
200 ,fNumberEnt(c.fNumberEnt)
201 ,fStatisticMean(c.fStatisticMean)
203 ,fDebugLevel(c.fDebugLevel)
204 ,fFitVoir(c.fFitVoir)
205 ,fMagneticField(c.fMagneticField)
207 ,fCurrentCoefE(c.fCurrentCoefE)
208 ,fCurrentCoefE2(c.fCurrentCoefE2)
211 ,fScaleFitFactor(c.fScaleFitFactor)
212 ,fEntriesCurrent(c.fEntriesCurrent)
213 ,fCountDet(c.fCountDet)
220 ,fCurrentCoefDetector(0x0)
221 ,fCurrentCoefDetector2(0x0)
229 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
231 //Current variables initialised
232 for (Int_t k = 0; k < 2; k++) {
233 fCurrentCoef[k] = 0.0;
234 fCurrentCoef2[k] = 0.0;
236 for (Int_t i = 0; i < 3; i++) {
240 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
241 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
243 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
244 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
246 fVectorFit.SetName(c.fVectorFit.GetName());
247 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
248 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
249 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
251 if (GetStack(detector) == 2) {
257 Float_t *coef = new Float_t[ntotal];
258 for (Int_t i = 0; i < ntotal; i++) {
259 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
261 fitInfo->SetCoef(coef);
262 fitInfo->SetDetector(detector);
263 fVectorFit.Add((TObject *) fitInfo);
265 fVectorFit.SetName(c.fVectorFit.GetName());
266 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
267 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
268 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
270 if (GetStack(detector) == 2) {
276 Float_t *coef = new Float_t[ntotal];
277 for (Int_t i = 0; i < ntotal; i++) {
278 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
280 fitInfo->SetCoef(coef);
281 fitInfo->SetDetector(detector);
282 fVectorFit2.Add((TObject *) fitInfo);
287 fGeo = new AliTRDgeometry();
290 //____________________________________________________________________________________
291 AliTRDCalibraFit::~AliTRDCalibraFit()
294 // AliTRDCalibraFit destructor
296 if ( fDebugStreamer ) delete fDebugStreamer;
297 if ( fCalDet ) delete fCalDet;
298 if ( fCalDet2 ) delete fCalDet2;
299 if ( fCalROC ) delete fCalROC;
300 if ( fCalROC2 ) delete fCalROC2;
301 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
302 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
304 fVectorFit2.Delete();
310 //_____________________________________________________________________________
311 void AliTRDCalibraFit::Destroy()
323 //_____________________________________________________________________________
324 void AliTRDCalibraFit::DestroyDebugStreamer()
327 // Delete DebugStreamer
330 if ( fDebugStreamer ) delete fDebugStreamer;
331 fDebugStreamer = 0x0;
334 //__________________________________________________________________________________
335 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
338 // From the drift velocity and t0
339 // return the position of the peak and maximum negative slope
342 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
343 Double_t widbins = 0.1; // 0.1 mus
345 //peak and maxnegslope in mus
346 Double_t begind = t0*widbins + fT0Shift0;
347 Double_t peakd = t0*widbins + fT0Shift1;
348 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
350 // peak and maxnegslope in timebin
351 begin = TMath::Nint(begind*widbins);
352 peak = TMath::Nint(peakd*widbins);
353 end = TMath::Nint(maxnegslope*widbins);
356 //____________Functions fit Online CH2d________________________________________
357 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
360 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
361 // calibration group normalized the resulted coefficients (to 1 normally)
364 // Set the calibration mode
365 const char *name = ch->GetTitle();
366 if(!SetModeCalibration(name,0)) return kFALSE;
368 // Number of Ybins (detectors or groups of pads)
369 Int_t nbins = ch->GetNbinsX();// charge
370 Int_t nybins = ch->GetNbinsY();// groups number
371 if (!InitFit(nybins,0)) {
377 fStatisticMean = 0.0;
379 fNumberFitSuccess = 0;
381 // Init fCountDet and fCount
382 InitfCountDetAndfCount(0);
383 // Beginning of the loop betwwen dect1 and dect2
384 for (Int_t idect = fDect1; idect < fDect2; idect++) {
385 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
386 UpdatefCountDetAndfCount(idect,0);
387 ReconstructFitRowMinRowMax(idect, 0);
389 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
390 projch->SetDirectory(0);
391 // Number of entries for this calibration group
392 Double_t nentries = 0.0;
394 for (Int_t k = 0; k < nbins; k++) {
395 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
396 nentries += ch->GetBinContent(binnb);
397 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
398 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
400 projch->SetEntries(nentries);
401 //printf("The number of entries for the group %d is %f\n",idect,nentries);
406 // Rebin and statistic stuff
408 projch = ReBin((TH1I *) projch);
410 // This detector has not enough statistics or was off
411 if (nentries <= fMinEntries) {
412 NotEnoughStatisticCH(idect);
413 if (fDebugLevel != 1) {
418 // Statistics of the group fitted
419 fStatisticMean += nentries;
424 case 0: FitMeanW((TH1 *) projch, nentries); break;
425 case 1: FitMean((TH1 *) projch, nentries, mean); break;
426 case 2: FitCH((TH1 *) projch, mean); break;
427 case 3: FitBisCH((TH1 *) projch, mean); break;
428 default: return kFALSE;
431 FillInfosFitCH(idect);
433 if (fDebugLevel != 1) {
438 if (fDebugLevel != 1) {
442 if (fNumberFit > 0) {
443 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
444 fStatisticMean = fStatisticMean / fNumberFit;
447 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
449 delete fDebugStreamer;
450 fDebugStreamer = 0x0;
454 //____________Functions fit Online CH2d________________________________________
455 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
458 // Reconstruct a 1D histo from the vectorCH for each calibration group,
459 // fit the histo, normalized the resulted coefficients (to 1 normally)
462 // Set the calibraMode
463 const char *name = calvect->GetNameCH();
464 if(!SetModeCalibration(name,0)) return kFALSE;
466 // Number of Xbins (detectors or groups of pads)
467 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
473 fStatisticMean = 0.0;
475 fNumberFitSuccess = 0;
477 // Init fCountDet and fCount
478 InitfCountDetAndfCount(0);
479 // Beginning of the loop between dect1 and dect2
480 for (Int_t idect = fDect1; idect < fDect2; idect++) {
481 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
482 UpdatefCountDetAndfCount(idect,0);
483 ReconstructFitRowMinRowMax(idect,0);
485 Double_t nentries = 0.0;
487 if(!calvect->GetCHEntries(fCountDet)) {
488 NotEnoughStatisticCH(idect);
494 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
495 projch->SetDirectory(0);
496 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
497 nentries += projch->GetBinContent(k+1);
498 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
504 //printf("The number of entries for the group %d is %f\n",idect,nentries);
507 projch = ReBin((TH1F *) projch);
509 // This detector has not enough statistics or was not found in VectorCH
510 if (nentries <= fMinEntries) {
511 NotEnoughStatisticCH(idect);
514 // Statistic of the histos fitted
515 fStatisticMean += nentries;
520 case 0: FitMeanW((TH1 *) projch, nentries); break;
521 case 1: FitMean((TH1 *) projch, nentries, mean); break;
522 case 2: FitCH((TH1 *) projch, mean); break;
523 case 3: FitBisCH((TH1 *) projch, mean); break;
524 default: return kFALSE;
527 FillInfosFitCH(idect);
530 if (fDebugLevel != 1) {
534 if (fNumberFit > 0) {
535 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));
536 fStatisticMean = fStatisticMean / fNumberFit;
539 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
541 delete fDebugStreamer;
542 fDebugStreamer = 0x0;
545 //________________functions fit Online PH2d____________________________________
546 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
549 // Take the 1D profiles (average pulse height), projections of the 2D PH
550 // on the Xaxis, for each calibration group
551 // Reconstruct a drift velocity
552 // A first calibration of T0 is also made using the same method
555 // Set the calibration mode
556 const char *name = ph->GetTitle();
557 if(!SetModeCalibration(name,1)) return kFALSE;
559 // Number of Xbins (detectors or groups of pads)
560 Int_t nbins = ph->GetNbinsX();// time
561 Int_t nybins = ph->GetNbinsY();// calibration group
562 if (!InitFit(nybins,1)) {
568 fStatisticMean = 0.0;
570 fNumberFitSuccess = 0;
572 // Init fCountDet and fCount
573 InitfCountDetAndfCount(1);
574 // Beginning of the loop
575 for (Int_t idect = fDect1; idect < fDect2; idect++) {
576 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
577 UpdatefCountDetAndfCount(idect,1);
578 ReconstructFitRowMinRowMax(idect,1);
580 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
581 projph->SetDirectory(0);
582 // Number of entries for this calibration group
583 Double_t nentries = 0;
584 for (Int_t k = 0; k < nbins; k++) {
585 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
586 nentries += ph->GetBinEntries(binnb);
591 //printf("The number of entries for the group %d is %f\n",idect,nentries);
592 // This detector has not enough statistics or was off
593 if (nentries <= fMinEntries) {
594 //printf("Not enough statistic!\n");
595 NotEnoughStatisticPH(idect,nentries);
596 if (fDebugLevel != 1) {
601 // Statistics of the histos fitted
603 fStatisticMean += nentries;
604 // Calcul of "real" coef
605 CalculVdriftCoefMean();
610 case 0: FitLagrangePoly((TH1 *) projph); break;
611 case 1: FitPente((TH1 *) projph); break;
612 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
613 default: return kFALSE;
615 // Fill the tree if end of a detector or only the pointer to the branch!!!
616 FillInfosFitPH(idect,nentries);
618 if (fDebugLevel != 1) {
623 if (fNumberFit > 0) {
624 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));
625 fStatisticMean = fStatisticMean / fNumberFit;
628 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
630 delete fDebugStreamer;
631 fDebugStreamer = 0x0;
634 //____________Functions fit Online PH2d________________________________________
635 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
638 // Reconstruct the average pulse height from the vectorPH for each
640 // Reconstruct a drift velocity
641 // A first calibration of T0 is also made using the same method (slope method)
644 // Set the calibration mode
645 const char *name = calvect->GetNamePH();
646 if(!SetModeCalibration(name,1)) return kFALSE;
648 // Number of Xbins (detectors or groups of pads)
649 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
655 fStatisticMean = 0.0;
657 fNumberFitSuccess = 0;
659 // Init fCountDet and fCount
660 InitfCountDetAndfCount(1);
661 // Beginning of the loop
662 for (Int_t idect = fDect1; idect < fDect2; idect++) {
663 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
664 UpdatefCountDetAndfCount(idect,1);
665 ReconstructFitRowMinRowMax(idect,1);
668 if(!calvect->GetPHEntries(fCountDet)) {
669 NotEnoughStatisticPH(idect,fEntriesCurrent);
674 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
675 projph->SetDirectory(0);
676 if(fEntriesCurrent > 0) fNumberEnt++;
677 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
678 // This detector has not enough statistics or was off
679 if (fEntriesCurrent <= fMinEntries) {
680 //printf("Not enough stat!\n");
681 NotEnoughStatisticPH(idect,fEntriesCurrent);
684 // Statistic of the histos fitted
686 fStatisticMean += fEntriesCurrent;
687 // Calcul of "real" coef
688 CalculVdriftCoefMean();
693 case 0: FitLagrangePoly((TH1 *) projph); break;
694 case 1: FitPente((TH1 *) projph); break;
695 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
696 default: return kFALSE;
698 // Fill the tree if end of a detector or only the pointer to the branch!!!
699 FillInfosFitPH(idect,fEntriesCurrent);
703 if (fNumberFit > 0) {
704 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));
705 fStatisticMean = fStatisticMean / fNumberFit;
708 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
710 delete fDebugStreamer;
711 fDebugStreamer = 0x0;
714 //____________Functions fit Online PRF2d_______________________________________
715 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
718 // Take the 1D profiles (pad response function), projections of the 2D PRF
719 // on the Xaxis, for each calibration group
720 // Fit with a gaussian to reconstruct the sigma of the pad response function
723 // Set the calibration mode
724 const char *name = prf->GetTitle();
725 if(!SetModeCalibration(name,2)) return kFALSE;
727 // Number of Ybins (detectors or groups of pads)
728 Int_t nybins = prf->GetNbinsY();// calibration groups
729 Int_t nbins = prf->GetNbinsX();// bins
730 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
731 if((nbg > 0) || (nbg == -1)) return kFALSE;
732 if (!InitFit(nybins,2)) {
738 fStatisticMean = 0.0;
740 fNumberFitSuccess = 0;
742 // Init fCountDet and fCount
743 InitfCountDetAndfCount(2);
744 // Beginning of the loop
745 for (Int_t idect = fDect1; idect < fDect2; idect++) {
746 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
747 UpdatefCountDetAndfCount(idect,2);
748 ReconstructFitRowMinRowMax(idect,2);
750 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
751 projprf->SetDirectory(0);
752 // Number of entries for this calibration group
753 Double_t nentries = 0;
754 for (Int_t k = 0; k < nbins; k++) {
755 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
756 nentries += prf->GetBinEntries(binnb);
758 if(nentries > 0) fNumberEnt++;
759 // This detector has not enough statistics or was off
760 if (nentries <= fMinEntries) {
761 NotEnoughStatisticPRF(idect);
762 if (fDebugLevel != 1) {
767 // Statistics of the histos fitted
769 fStatisticMean += nentries;
770 // Calcul of "real" coef
775 case 0: FitPRF((TH1 *) projprf); break;
776 case 1: RmsPRF((TH1 *) projprf); break;
777 default: return kFALSE;
779 // Fill the tree if end of a detector or only the pointer to the branch!!!
780 FillInfosFitPRF(idect);
782 if (fDebugLevel != 1) {
787 if (fNumberFit > 0) {
788 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
789 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
790 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
791 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
792 fStatisticMean = fStatisticMean / fNumberFit;
795 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
797 delete fDebugStreamer;
798 fDebugStreamer = 0x0;
801 //____________Functions fit Online PRF2d_______________________________________
802 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
805 // Take the 1D profiles (pad response function), projections of the 2D PRF
806 // on the Xaxis, for each calibration group
807 // Fit with a gaussian to reconstruct the sigma of the pad response function
810 // Set the calibration mode
811 const char *name = prf->GetTitle();
812 if(!SetModeCalibration(name,2)) return kFALSE;
814 // Number of Ybins (detectors or groups of pads)
815 TAxis *xprf = prf->GetXaxis();
816 TAxis *yprf = prf->GetYaxis();
817 Int_t nybins = yprf->GetNbins();// calibration groups
818 Int_t nbins = xprf->GetNbins();// bins
819 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
820 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
821 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
822 if(nbg == -1) return kFALSE;
823 if(nbg > 0) fMethod = 1;
825 if (!InitFit(nybins,2)) {
831 fStatisticMean = 0.0;
833 fNumberFitSuccess = 0;
835 // Init fCountDet and fCount
836 InitfCountDetAndfCount(2);
837 // Beginning of the loop
838 for (Int_t idect = fDect1; idect < fDect2; idect++) {
839 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
840 UpdatefCountDetAndfCount(idect,2);
841 ReconstructFitRowMinRowMax(idect,2);
842 // Build the array of entries and sum
843 TArrayD arraye = TArrayD(nbins);
844 TArrayD arraym = TArrayD(nbins);
845 TArrayD arrayme = TArrayD(nbins);
846 Double_t nentries = 0;
847 //printf("nbins %d\n",nbins);
848 for (Int_t k = 0; k < nbins; k++) {
849 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
850 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
851 Double_t mean = (Double_t)prf->GetBinContent(binnb);
852 Double_t error = (Double_t)prf->GetBinError(binnb);
853 //printf("for %d we have %f\n",k,entries);
855 arraye.AddAt(entries,k);
856 arraym.AddAt(mean,k);
857 arrayme.AddAt(error,k);
859 if(nentries > 0) fNumberEnt++;
860 //printf("The number of entries for the group %d is %f\n",idect,nentries);
861 // This detector has not enough statistics or was off
862 if (nentries <= fMinEntries) {
863 NotEnoughStatisticPRF(idect);
866 // Statistics of the histos fitted
868 fStatisticMean += nentries;
869 // Calcul of "real" coef
874 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
875 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
876 default: return kFALSE;
878 // Fill the tree if end of a detector or only the pointer to the branch!!!
879 FillInfosFitPRF(idect);
882 if (fNumberFit > 0) {
883 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
884 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
885 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
886 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
887 fStatisticMean = fStatisticMean / fNumberFit;
890 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
892 delete fDebugStreamer;
893 fDebugStreamer = 0x0;
896 //____________Functions fit Online PRF2d_______________________________________
897 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
900 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
901 // each calibration group
902 // Fit with a gaussian to reconstruct the sigma of the pad response function
905 // Set the calibra mode
906 const char *name = calvect->GetNamePRF();
907 if(!SetModeCalibration(name,2)) return kFALSE;
908 //printf("test0 %s\n",name);
910 // Number of Xbins (detectors or groups of pads)
911 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
916 ///printf("test2\n");
919 fStatisticMean = 0.0;
921 fNumberFitSuccess = 0;
923 // Init fCountDet and fCount
924 InitfCountDetAndfCount(2);
925 // Beginning of the loop
926 for (Int_t idect = fDect1; idect < fDect2; idect++) {
927 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
928 UpdatefCountDetAndfCount(idect,2);
929 ReconstructFitRowMinRowMax(idect,2);
932 if(!calvect->GetPRFEntries(fCountDet)) {
933 NotEnoughStatisticPRF(idect);
936 TString tname("PRF");
938 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
939 projprf->SetDirectory(0);
940 if(fEntriesCurrent > 0) fNumberEnt++;
941 // This detector has not enough statistics or was off
942 if (fEntriesCurrent <= fMinEntries) {
943 NotEnoughStatisticPRF(idect);
946 // Statistic of the histos fitted
948 fStatisticMean += fEntriesCurrent;
949 // Calcul of "real" coef
954 case 1: FitPRF((TH1 *) projprf); break;
955 case 2: RmsPRF((TH1 *) projprf); break;
956 default: return kFALSE;
958 // Fill the tree if end of a detector or only the pointer to the branch!!!
959 FillInfosFitPRF(idect);
962 if (fNumberFit > 0) {
963 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));
966 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
968 delete fDebugStreamer;
969 fDebugStreamer = 0x0;
972 //____________Functions fit Online PRF2d_______________________________________
973 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
976 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
977 // each calibration group
978 // Fit with a gaussian to reconstruct the sigma of the pad response function
981 // Set the calibra mode
982 const char *name = calvect->GetNamePRF();
983 if(!SetModeCalibration(name,2)) return kFALSE;
984 //printf("test0 %s\n",name);
985 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
986 //printf("test1 %d\n",nbg);
987 if(nbg == -1) return kFALSE;
988 if(nbg > 0) fMethod = 1;
990 // Number of Xbins (detectors or groups of pads)
991 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
999 fStatisticMean = 0.0;
1001 fNumberFitSuccess = 0;
1005 Double_t *arrayx = 0;
1006 Double_t *arraye = 0;
1007 Double_t *arraym = 0;
1008 Double_t *arrayme = 0;
1009 Float_t lowedge = 0.0;
1010 Float_t upedge = 0.0;
1011 // Init fCountDet and fCount
1012 InitfCountDetAndfCount(2);
1013 // Beginning of the loop
1014 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1015 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1016 UpdatefCountDetAndfCount(idect,2);
1017 ReconstructFitRowMinRowMax(idect,2);
1019 fEntriesCurrent = 0;
1020 if(!calvect->GetPRFEntries(fCountDet)) {
1021 NotEnoughStatisticPRF(idect);
1024 TString tname("PRF");
1026 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1027 nbins = projprftree->GetN();
1028 arrayx = (Double_t *)projprftree->GetX();
1029 arraye = (Double_t *)projprftree->GetEX();
1030 arraym = (Double_t *)projprftree->GetY();
1031 arrayme = (Double_t *)projprftree->GetEY();
1032 Float_t step = arrayx[1]-arrayx[0];
1033 lowedge = arrayx[0] - step/2.0;
1034 upedge = arrayx[(nbins-1)] + step/2.0;
1035 //printf("nbins est %d\n",nbins);
1036 for(Int_t k = 0; k < nbins; k++){
1037 fEntriesCurrent += (Int_t)arraye[k];
1038 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1039 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1041 if(fEntriesCurrent > 0) fNumberEnt++;
1042 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1043 // This detector has not enough statistics or was off
1044 if (fEntriesCurrent <= fMinEntries) {
1045 NotEnoughStatisticPRF(idect);
1048 // Statistic of the histos fitted
1050 fStatisticMean += fEntriesCurrent;
1051 // Calcul of "real" coef
1052 CalculPRFCoefMean();
1056 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1057 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1058 default: return kFALSE;
1060 // Fill the tree if end of a detector or only the pointer to the branch!!!
1061 FillInfosFitPRF(idect);
1064 if (fNumberFit > 0) {
1065 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));
1068 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1070 delete fDebugStreamer;
1071 fDebugStreamer = 0x0;
1074 //____________Functions fit Online CH2d________________________________________
1075 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1078 // The linear method
1081 fStatisticMean = 0.0;
1083 fNumberFitSuccess = 0;
1085 if(!InitFitLinearFitter()) return kFALSE;
1088 for(Int_t idet = 0; idet < 540; idet++){
1091 //printf("detector number %d\n",idet);
1096 fEntriesCurrent = 0;
1098 Bool_t here = calivdli->GetParam(idet,¶m);
1099 Bool_t heree = calivdli->GetError(idet,&error);
1100 //printf("here %d and heree %d\n",here, heree);
1102 fEntriesCurrent = (Int_t) error[2];
1105 //printf("Number of entries %d\n",fEntriesCurrent);
1106 // Nothing found or not enough statistic
1107 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1108 NotEnoughStatisticLinearFitter();
1115 fStatisticMean += fEntriesCurrent;
1118 if((-(param[1])) <= 0.0) {
1119 NotEnoughStatisticLinearFitter();
1123 // CalculDatabaseVdriftandTan
1124 CalculVdriftLorentzCoef();
1127 fNumberFitSuccess ++;
1129 // Put the fCurrentCoef
1130 fCurrentCoef[0] = -param[1];
1131 // here the database must be the one of the reconstruction for the lorentz angle....
1132 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1133 fCurrentCoefE = error[1];
1134 fCurrentCoefE2 = error[0];
1135 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1136 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1140 FillInfosFitLinearFitter();
1145 if (fNumberFit > 0) {
1146 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));
1149 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1151 delete fDebugStreamer;
1152 fDebugStreamer = 0x0;
1156 //____________Functions for seeing if the pad is really okey___________________
1157 //_____________________________________________________________________________
1158 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1161 // Get numberofgroupsprf
1165 const Char_t *pattern0 = "Ngp0";
1166 const Char_t *pattern1 = "Ngp1";
1167 const Char_t *pattern2 = "Ngp2";
1168 const Char_t *pattern3 = "Ngp3";
1169 const Char_t *pattern4 = "Ngp4";
1170 const Char_t *pattern5 = "Ngp5";
1171 const Char_t *pattern6 = "Ngp6";
1174 if (strstr(nametitle,pattern0)) {
1177 if (strstr(nametitle,pattern1)) {
1180 if (strstr(nametitle,pattern2)) {
1183 if (strstr(nametitle,pattern3)) {
1186 if (strstr(nametitle,pattern4)) {
1189 if (strstr(nametitle,pattern5)) {
1192 if (strstr(nametitle,pattern6)){
1199 //_____________________________________________________________________________
1200 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1203 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1204 // corresponding to the given name
1207 if(!SetNzFromTObject(name,i)) return kFALSE;
1208 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1213 //_____________________________________________________________________________
1214 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1217 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1218 // corresponding to the given TObject
1222 const Char_t *patternrphi0 = "Nrphi0";
1223 const Char_t *patternrphi1 = "Nrphi1";
1224 const Char_t *patternrphi2 = "Nrphi2";
1225 const Char_t *patternrphi3 = "Nrphi3";
1226 const Char_t *patternrphi4 = "Nrphi4";
1227 const Char_t *patternrphi5 = "Nrphi5";
1228 const Char_t *patternrphi6 = "Nrphi6";
1231 const Char_t *patternrphi10 = "Nrphi10";
1232 const Char_t *patternrphi100 = "Nrphi100";
1233 const Char_t *patternz10 = "Nz10";
1234 const Char_t *patternz100 = "Nz100";
1237 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1238 fCalibraMode->SetAllTogether(i);
1240 if (fDebugLevel > 1) {
1241 AliInfo(Form("fNbDet %d and 100",fNbDet));
1245 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1246 fCalibraMode->SetPerSuperModule(i);
1248 if (fDebugLevel > 1) {
1249 AliInfo(Form("fNDet %d and 100",fNbDet));
1254 if (strstr(name,patternrphi0)) {
1255 fCalibraMode->SetNrphi(i ,0);
1256 if (fDebugLevel > 1) {
1257 AliInfo(Form("fNbDet %d and 0",fNbDet));
1261 if (strstr(name,patternrphi1)) {
1262 fCalibraMode->SetNrphi(i, 1);
1263 if (fDebugLevel > 1) {
1264 AliInfo(Form("fNbDet %d and 1",fNbDet));
1268 if (strstr(name,patternrphi2)) {
1269 fCalibraMode->SetNrphi(i, 2);
1270 if (fDebugLevel > 1) {
1271 AliInfo(Form("fNbDet %d and 2",fNbDet));
1275 if (strstr(name,patternrphi3)) {
1276 fCalibraMode->SetNrphi(i, 3);
1277 if (fDebugLevel > 1) {
1278 AliInfo(Form("fNbDet %d and 3",fNbDet));
1282 if (strstr(name,patternrphi4)) {
1283 fCalibraMode->SetNrphi(i, 4);
1284 if (fDebugLevel > 1) {
1285 AliInfo(Form("fNbDet %d and 4",fNbDet));
1289 if (strstr(name,patternrphi5)) {
1290 fCalibraMode->SetNrphi(i, 5);
1291 if (fDebugLevel > 1) {
1292 AliInfo(Form("fNbDet %d and 5",fNbDet));
1296 if (strstr(name,patternrphi6)) {
1297 fCalibraMode->SetNrphi(i, 6);
1298 if (fDebugLevel > 1) {
1299 AliInfo(Form("fNbDet %d and 6",fNbDet));
1304 if (fDebugLevel > 1) {
1305 AliInfo(Form("fNbDet %d and rest",fNbDet));
1307 fCalibraMode->SetNrphi(i ,0);
1311 //_____________________________________________________________________________
1312 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1315 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1316 // corresponding to the given TObject
1320 const Char_t *patternz0 = "Nz0";
1321 const Char_t *patternz1 = "Nz1";
1322 const Char_t *patternz2 = "Nz2";
1323 const Char_t *patternz3 = "Nz3";
1324 const Char_t *patternz4 = "Nz4";
1326 const Char_t *patternrphi10 = "Nrphi10";
1327 const Char_t *patternrphi100 = "Nrphi100";
1328 const Char_t *patternz10 = "Nz10";
1329 const Char_t *patternz100 = "Nz100";
1331 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1332 fCalibraMode->SetAllTogether(i);
1334 if (fDebugLevel > 1) {
1335 AliInfo(Form("fNbDet %d and 100",fNbDet));
1339 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1340 fCalibraMode->SetPerSuperModule(i);
1342 if (fDebugLevel > 1) {
1343 AliInfo(Form("fNbDet %d and 10",fNbDet));
1347 if (strstr(name,patternz0)) {
1348 fCalibraMode->SetNz(i, 0);
1349 if (fDebugLevel > 1) {
1350 AliInfo(Form("fNbDet %d and 0",fNbDet));
1354 if (strstr(name,patternz1)) {
1355 fCalibraMode->SetNz(i ,1);
1356 if (fDebugLevel > 1) {
1357 AliInfo(Form("fNbDet %d and 1",fNbDet));
1361 if (strstr(name,patternz2)) {
1362 fCalibraMode->SetNz(i ,2);
1363 if (fDebugLevel > 1) {
1364 AliInfo(Form("fNbDet %d and 2",fNbDet));
1368 if (strstr(name,patternz3)) {
1369 fCalibraMode->SetNz(i ,3);
1370 if (fDebugLevel > 1) {
1371 AliInfo(Form("fNbDet %d and 3",fNbDet));
1375 if (strstr(name,patternz4)) {
1376 fCalibraMode->SetNz(i ,4);
1377 if (fDebugLevel > 1) {
1378 AliInfo(Form("fNbDet %d and 4",fNbDet));
1383 if (fDebugLevel > 1) {
1384 AliInfo(Form("fNbDet %d and rest",fNbDet));
1386 fCalibraMode->SetNz(i ,0);
1389 //______________________________________________________________________
1390 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1392 // ofwhat is equaled to 0: mean value of all passing detectors
1393 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1396 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1398 AliInfo("The Vector Fit is not complete!");
1401 Int_t detector = -1;
1403 Float_t value = 0.0;
1405 /////////////////////////////////
1406 // Calculate the mean values
1407 ////////////////////////////////
1409 ////////////////////////
1410 Double_t meanAll = 0.0;
1411 Double_t meanSupermodule[18];
1412 Double_t meanDetector[540];
1414 Int_t countSupermodule[18];
1415 Int_t countDetector[540];
1416 for(Int_t sm = 0; sm < 18; sm++){
1417 meanSupermodule[sm] = 0.0;
1418 countSupermodule[sm] = 0;
1420 for(Int_t det = 0; det < 540; det++){
1421 meanDetector[det] = 0.0;
1422 countDetector[det] = 0;
1426 for (Int_t k = 0; k < loop; k++) {
1427 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1428 sector = GetSector(detector);
1430 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1432 meanDetector[detector] += value;
1433 countDetector[detector]++;
1434 meanSupermodule[sector] += value;
1435 countSupermodule[sector]++;
1441 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1442 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1443 for (Int_t row = 0; row < rowMax; row++) {
1444 for (Int_t col = 0; col < colMax; col++) {
1445 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1447 meanDetector[detector] += value;
1448 countDetector[detector]++;
1449 meanSupermodule[sector] += value;
1450 countSupermodule[sector]++;
1459 if(countAll > 0) meanAll = meanAll/countAll;
1460 for(Int_t sm = 0; sm < 18; sm++){
1461 if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1463 for(Int_t det = 0; det < 540; det++){
1464 if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1466 // Put the mean value for the no-fitted
1467 /////////////////////////////////////////////
1468 for (Int_t k = 0; k < loop; k++) {
1469 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1470 sector = GetSector(detector);
1471 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1472 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1473 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1475 for (Int_t row = 0; row < rowMax; row++) {
1476 for (Int_t col = 0; col < colMax; col++) {
1477 value = coef[(Int_t)(col*rowMax+row)];
1479 if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1481 if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1482 else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1483 else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1487 if(fDebugLevel > 1){
1489 if ( !fDebugStreamer ) {
1491 TDirectory *backup = gDirectory;
1492 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1493 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1496 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1498 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1499 "detector="<<detector<<
1512 //______________________________________________________________________
1513 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1515 // ofwhat is equaled to 0: mean value of all passing detectors
1516 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1519 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1521 AliInfo("The Vector Fit is not complete!");
1524 Int_t detector = -1;
1526 Float_t value = 0.0;
1528 /////////////////////////////////
1529 // Calculate the mean values
1530 ////////////////////////////////
1532 ////////////////////////
1533 Double_t meanAll = 0.0;
1534 Double_t meanSupermodule[18];
1535 Double_t meanDetector[540];
1537 Int_t countSupermodule[18];
1538 Int_t countDetector[540];
1539 for(Int_t sm = 0; sm < 18; sm++){
1540 meanSupermodule[sm] = 0.0;
1541 countSupermodule[sm] = 0;
1543 for(Int_t det = 0; det < 540; det++){
1544 meanDetector[det] = 0.0;
1545 countDetector[det] = 0;
1549 for (Int_t k = 0; k < loop; k++) {
1550 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1551 sector = GetSector(detector);
1553 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1555 meanDetector[detector] += value;
1556 countDetector[detector]++;
1557 meanSupermodule[sector] += value;
1558 countSupermodule[sector]++;
1564 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1565 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1566 for (Int_t row = 0; row < rowMax; row++) {
1567 for (Int_t col = 0; col < colMax; col++) {
1568 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1570 meanDetector[detector] += value;
1571 countDetector[detector]++;
1572 meanSupermodule[sector] += value;
1573 countSupermodule[sector]++;
1582 if(countAll > 0) meanAll = meanAll/countAll;
1583 for(Int_t sm = 0; sm < 18; sm++){
1584 if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1586 for(Int_t det = 0; det < 540; det++){
1587 if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1589 // Put the mean value for the no-fitted
1590 /////////////////////////////////////////////
1591 for (Int_t k = 0; k < loop; k++) {
1592 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1593 sector = GetSector(detector);
1594 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1595 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1596 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1598 for (Int_t row = 0; row < rowMax; row++) {
1599 for (Int_t col = 0; col < colMax; col++) {
1600 value = coef[(Int_t)(col*rowMax+row)];
1602 if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1604 if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1605 else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1606 else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1610 if(fDebugLevel > 1){
1612 if ( !fDebugStreamer ) {
1614 TDirectory *backup = gDirectory;
1615 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1616 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1619 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1621 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1622 "detector="<<detector<<
1635 //_____________________________________________________________________________
1636 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1639 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1640 // It takes the mean value of the coefficients per detector
1641 // This object has to be written in the database
1644 // Create the DetObject
1645 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1647 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1648 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1649 Int_t detector = -1;
1650 Float_t value = 0.0;
1653 for (Int_t k = 0; k < loop; k++) {
1654 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1657 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1661 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1662 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1663 for (Int_t row = 0; row < rowMax; row++) {
1664 for (Int_t col = 0; col < colMax; col++) {
1665 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1666 mean += TMath::Abs(value);
1670 if(count > 0) mean = mean/count;
1672 object->SetValue(detector,mean);
1677 //_____________________________________________________________________________
1678 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1681 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1682 // It takes the mean value of the coefficients per detector
1683 // This object has to be written in the database
1686 // Create the DetObject
1687 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1690 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1691 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1692 Int_t detector = -1;
1693 Float_t value = 0.0;
1695 for (Int_t k = 0; k < loop; k++) {
1696 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1699 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1700 if(!meanOtherBefore){
1701 if(value > 0) value = value*scaleFitFactor;
1703 else value = value*scaleFitFactor;
1704 mean = TMath::Abs(value);
1708 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1709 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1710 for (Int_t row = 0; row < rowMax; row++) {
1711 for (Int_t col = 0; col < colMax; col++) {
1712 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1713 if(!meanOtherBefore) {
1714 if(value > 0) value = value*scaleFitFactor;
1716 else value = value*scaleFitFactor;
1717 mean += TMath::Abs(value);
1721 if(count > 0) mean = mean/count;
1723 object->SetValue(detector,mean);
1728 //_____________________________________________________________________________
1729 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1732 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1733 // It takes the min value of the coefficients per detector
1734 // This object has to be written in the database
1737 // Create the DetObject
1738 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1740 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1741 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1742 Int_t detector = -1;
1743 Float_t value = 0.0;
1745 for (Int_t k = 0; k < loop; k++) {
1746 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1747 Float_t min = 100.0;
1749 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1751 if(value > 70.0) value = value-100.0;
1756 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1757 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1758 for (Int_t row = 0; row < rowMax; row++) {
1759 for (Int_t col = 0; col < colMax; col++) {
1760 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1762 if(value > 70.0) value = value-100.0;
1764 if(min > value) min = value;
1768 object->SetValue(detector,min);
1774 //_____________________________________________________________________________
1775 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
1778 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1779 // It takes the min value of the coefficients per detector
1780 // This object has to be written in the database
1783 // Create the DetObject
1784 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1787 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1788 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1789 Int_t detector = -1;
1790 Float_t value = 0.0;
1792 for (Int_t k = 0; k < loop; k++) {
1793 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1795 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1796 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1797 Float_t min = 100.0;
1798 for (Int_t row = 0; row < rowMax; row++) {
1799 for (Int_t col = 0; col < colMax; col++) {
1800 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1801 mean += -TMath::Abs(value);
1805 if(count > 0) mean = mean/count;
1807 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1808 object->SetValue(detector,-TMath::Abs(value));
1814 //_____________________________________________________________________________
1815 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
1818 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1819 // You need first to create the object for the detectors,
1820 // where the mean value is put.
1821 // This object has to be written in the database
1824 // Create the DetObject
1825 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1828 for(Int_t k = 0; k < 540; k++){
1829 AliTRDCalROC *calROC = object->GetCalROC(k);
1830 Int_t nchannels = calROC->GetNchannels();
1831 for(Int_t ch = 0; ch < nchannels; ch++){
1832 calROC->SetValue(ch,1.0);
1838 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1839 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1840 Int_t detector = -1;
1841 Float_t value = 0.0;
1843 for (Int_t k = 0; k < loop; k++) {
1844 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1845 AliTRDCalROC *calROC = object->GetCalROC(detector);
1846 Float_t mean = detobject->GetValue(detector);
1847 if(mean == 0) continue;
1848 Int_t rowMax = calROC->GetNrows();
1849 Int_t colMax = calROC->GetNcols();
1850 for (Int_t row = 0; row < rowMax; row++) {
1851 for (Int_t col = 0; col < colMax; col++) {
1852 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1853 if(value > 0) value = value*scaleFitFactor;
1854 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1862 //_____________________________________________________________________________
1863 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1866 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1867 // You need first to create the object for the detectors,
1868 // where the mean value is put.
1869 // This object has to be written in the database
1872 // Create the DetObject
1873 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1876 for(Int_t k = 0; k < 540; k++){
1877 AliTRDCalROC *calROC = object->GetCalROC(k);
1878 Int_t nchannels = calROC->GetNchannels();
1879 for(Int_t ch = 0; ch < nchannels; ch++){
1880 calROC->SetValue(ch,1.0);
1886 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1887 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1888 Int_t detector = -1;
1889 Float_t value = 0.0;
1891 for (Int_t k = 0; k < loop; k++) {
1892 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1893 AliTRDCalROC *calROC = object->GetCalROC(detector);
1894 Float_t mean = detobject->GetValue(detector);
1895 if(mean == 0) continue;
1896 Int_t rowMax = calROC->GetNrows();
1897 Int_t colMax = calROC->GetNcols();
1898 for (Int_t row = 0; row < rowMax; row++) {
1899 for (Int_t col = 0; col < colMax; col++) {
1900 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1901 calROC->SetValue(col,row,TMath::Abs(value)/mean);
1909 //_____________________________________________________________________________
1910 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1913 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1914 // You need first to create the object for the detectors,
1915 // where the mean value is put.
1916 // This object has to be written in the database
1919 // Create the DetObject
1920 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1923 for(Int_t k = 0; k < 540; k++){
1924 AliTRDCalROC *calROC = object->GetCalROC(k);
1925 Int_t nchannels = calROC->GetNchannels();
1926 for(Int_t ch = 0; ch < nchannels; ch++){
1927 calROC->SetValue(ch,0.0);
1933 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1934 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1935 Int_t detector = -1;
1936 Float_t value = 0.0;
1938 for (Int_t k = 0; k < loop; k++) {
1939 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1940 AliTRDCalROC *calROC = object->GetCalROC(detector);
1941 Float_t min = detobject->GetValue(detector);
1942 Int_t rowMax = calROC->GetNrows();
1943 Int_t colMax = calROC->GetNcols();
1944 for (Int_t row = 0; row < rowMax; row++) {
1945 for (Int_t col = 0; col < colMax; col++) {
1946 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1948 if(value > 70.0) value = value - 100.0;
1950 calROC->SetValue(col,row,value-min);
1958 //_____________________________________________________________________________
1959 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
1962 // It Creates the AliTRDCalPad object from AliTRDFitInfo
1963 // This object has to be written in the database
1966 // Create the DetObject
1967 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1969 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1970 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1971 Int_t detector = -1;
1972 Float_t value = 0.0;
1974 for (Int_t k = 0; k < loop; k++) {
1975 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1976 AliTRDCalROC *calROC = object->GetCalROC(detector);
1977 Int_t rowMax = calROC->GetNrows();
1978 Int_t colMax = calROC->GetNcols();
1979 for (Int_t row = 0; row < rowMax; row++) {
1980 for (Int_t col = 0; col < colMax; col++) {
1981 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1982 calROC->SetValue(col,row,TMath::Abs(value));
1990 //_____________________________________________________________________________
1991 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
1994 // It Creates the AliTRDCalDet object from AliTRDFitInfo
1995 // 0 successful fit 1 not successful fit
1996 // mean is the mean value over the successful fit
1997 // do not use it for t0: no meaning
2000 // Create the CalObject
2001 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2005 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2007 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2008 for(Int_t k = 0; k < 540; k++){
2009 object->SetValue(k,1.0);
2012 Int_t detector = -1;
2013 Float_t value = 0.0;
2015 for (Int_t k = 0; k < loop; k++) {
2016 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2017 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2018 if(value <= 0) object->SetValue(detector,1.0);
2020 object->SetValue(detector,0.0);
2025 if(count > 0) mean /= count;
2028 //_____________________________________________________________________________
2029 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2032 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2033 // 0 not successful fit 1 successful fit
2034 // mean mean value over the successful fit
2037 // Create the CalObject
2038 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2042 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2044 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2045 for(Int_t k = 0; k < 540; k++){
2046 AliTRDCalROC *calROC = object->GetCalROC(k);
2047 Int_t nchannels = calROC->GetNchannels();
2048 for(Int_t ch = 0; ch < nchannels; ch++){
2049 calROC->SetValue(ch,1.0);
2053 Int_t detector = -1;
2054 Float_t value = 0.0;
2056 for (Int_t k = 0; k < loop; k++) {
2057 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2058 AliTRDCalROC *calROC = object->GetCalROC(detector);
2059 Int_t nchannels = calROC->GetNchannels();
2060 for (Int_t ch = 0; ch < nchannels; ch++) {
2061 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2062 if(value <= 0) calROC->SetValue(ch,1.0);
2064 calROC->SetValue(ch,0.0);
2070 if(count > 0) mean /= count;
2073 //_____________________________________________________________________________
2074 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2077 // Set FitPH if 1 then each detector will be fitted
2080 if (periodeFitPH > 0) {
2081 fFitPHPeriode = periodeFitPH;
2084 AliInfo("periodeFitPH must be higher than 0!");
2088 //_____________________________________________________________________________
2089 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2092 // The fit of the deposited charge distribution begins at
2093 // histo->Mean()/beginFitCharge
2094 // You can here set beginFitCharge
2097 if (beginFitCharge > 0) {
2098 fBeginFitCharge = beginFitCharge;
2101 AliInfo("beginFitCharge must be strict positif!");
2106 //_____________________________________________________________________________
2107 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2110 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2111 // You can here set t0Shift0
2115 fT0Shift0 = t0Shift;
2118 AliInfo("t0Shift0 must be strict positif!");
2123 //_____________________________________________________________________________
2124 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2127 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2128 // You can here set t0Shift1
2132 fT0Shift1 = t0Shift;
2135 AliInfo("t0Shift must be strict positif!");
2140 //_____________________________________________________________________________
2141 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2144 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2145 // You can here set rangeFitPRF
2148 if ((rangeFitPRF > 0) &&
2149 (rangeFitPRF <= 1.5)) {
2150 fRangeFitPRF = rangeFitPRF;
2153 AliInfo("rangeFitPRF must be between 0 and 1.0");
2158 //_____________________________________________________________________________
2159 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2162 // Minimum entries for fitting
2165 if (minEntries > 0) {
2166 fMinEntries = minEntries;
2169 AliInfo("fMinEntries must be >= 0.");
2174 //_____________________________________________________________________________
2175 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2178 // Rebin with rebin time less bins the Ch histo
2179 // You can set here rebin that should divide the number of bins of CH histo
2184 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2187 AliInfo("You have to choose a positiv value!");
2191 //_____________________________________________________________________________
2192 Bool_t AliTRDCalibraFit::FillVectorFit()
2195 // For the Fit functions fill the vector Fit
2198 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2201 if (GetStack(fCountDet) == 2) {
2208 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2209 Float_t *coef = new Float_t[ntotal];
2210 for (Int_t i = 0; i < ntotal; i++) {
2211 coef[i] = fCurrentCoefDetector[i];
2214 Int_t detector = fCountDet;
2216 fitInfo->SetCoef(coef);
2217 fitInfo->SetDetector(detector);
2218 fVectorFit.Add((TObject *) fitInfo);
2223 //_____________________________________________________________________________
2224 Bool_t AliTRDCalibraFit::FillVectorFit2()
2227 // For the Fit functions fill the vector Fit
2230 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2233 if (GetStack(fCountDet) == 2) {
2240 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2241 Float_t *coef = new Float_t[ntotal];
2242 for (Int_t i = 0; i < ntotal; i++) {
2243 coef[i] = fCurrentCoefDetector2[i];
2246 Int_t detector = fCountDet;
2248 fitInfo->SetCoef(coef);
2249 fitInfo->SetDetector(detector);
2250 fVectorFit2.Add((TObject *) fitInfo);
2255 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2256 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2259 // Init the number of expected bins and fDect1[i] fDect2[i]
2262 gStyle->SetPalette(1);
2263 gStyle->SetOptStat(1111);
2264 gStyle->SetPadBorderMode(0);
2265 gStyle->SetCanvasColor(10);
2266 gStyle->SetPadLeftMargin(0.13);
2267 gStyle->SetPadRightMargin(0.01);
2269 // Mode groups of pads: the total number of bins!
2270 CalculNumberOfBinsExpected(i);
2272 // Quick verification that we have the good pad calibration mode!
2273 if (fNumberOfBinsExpected != nbins) {
2274 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2278 // Security for fDebug 3 and 4
2279 if ((fDebugLevel >= 3) &&
2283 AliInfo("This detector doesn't exit!");
2287 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2288 CalculDect1Dect2(i);
2293 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2294 Bool_t AliTRDCalibraFit::InitFitCH()
2297 // Init the fVectorFitCH for normalisation
2298 // Init the histo for debugging
2303 fScaleFitFactor = 0.0;
2304 fCurrentCoefDetector = new Float_t[2304];
2305 for (Int_t k = 0; k < 2304; k++) {
2306 fCurrentCoefDetector[k] = 0.0;
2308 fVectorFit.SetName("gainfactorscoefficients");
2310 // fDebug == 0 nothing
2311 // fDebug == 1 and fFitVoir no histo
2312 if (fDebugLevel == 1) {
2313 if(!CheckFitVoir()) return kFALSE;
2315 //Get the CalDet object
2317 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2319 AliInfo("Could not get calibDB");
2322 if(fCalDet) delete fCalDet;
2323 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2326 Float_t devalue = 1.0;
2327 if(fCalDet) delete fCalDet;
2328 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2329 for(Int_t k = 0; k < 540; k++){
2330 fCalDet->SetValue(k,devalue);
2336 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2337 Bool_t AliTRDCalibraFit::InitFitPH()
2340 // Init the arrays of results
2341 // Init the histos for debugging
2345 fVectorFit.SetName("driftvelocitycoefficients");
2346 fVectorFit2.SetName("t0coefficients");
2348 fCurrentCoefDetector = new Float_t[2304];
2349 for (Int_t k = 0; k < 2304; k++) {
2350 fCurrentCoefDetector[k] = 0.0;
2353 fCurrentCoefDetector2 = new Float_t[2304];
2354 for (Int_t k = 0; k < 2304; k++) {
2355 fCurrentCoefDetector2[k] = 0.0;
2358 //fDebug == 0 nothing
2359 // fDebug == 1 and fFitVoir no histo
2360 if (fDebugLevel == 1) {
2361 if(!CheckFitVoir()) return kFALSE;
2363 //Get the CalDet object
2365 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2367 AliInfo("Could not get calibDB");
2370 if(fCalDet) delete fCalDet;
2371 if(fCalDet2) delete fCalDet2;
2372 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2373 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2376 Float_t devalue = 1.5;
2377 Float_t devalue2 = 0.0;
2378 if(fCalDet) delete fCalDet;
2379 if(fCalDet2) delete fCalDet2;
2380 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2381 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2382 for(Int_t k = 0; k < 540; k++){
2383 fCalDet->SetValue(k,devalue);
2384 fCalDet2->SetValue(k,devalue2);
2389 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2390 Bool_t AliTRDCalibraFit::InitFitPRF()
2393 // Init the calibration mode (Nz, Nrphi), the histograms for
2394 // debugging the fit methods if fDebug > 0,
2398 fVectorFit.SetName("prfwidthcoefficients");
2400 fCurrentCoefDetector = new Float_t[2304];
2401 for (Int_t k = 0; k < 2304; k++) {
2402 fCurrentCoefDetector[k] = 0.0;
2405 // fDebug == 0 nothing
2406 // fDebug == 1 and fFitVoir no histo
2407 if (fDebugLevel == 1) {
2408 if(!CheckFitVoir()) return kFALSE;
2412 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2413 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2416 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2421 fCurrentCoefDetector = new Float_t[2304];
2422 fCurrentCoefDetector2 = new Float_t[2304];
2423 for (Int_t k = 0; k < 2304; k++) {
2424 fCurrentCoefDetector[k] = 0.0;
2425 fCurrentCoefDetector2[k] = 0.0;
2428 //printf("test0\n");
2430 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2432 AliInfo("Could not get calibDB");
2436 //Get the CalDet object
2438 if(fCalDet) delete fCalDet;
2439 if(fCalDet2) delete fCalDet2;
2440 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2441 //printf("test1\n");
2442 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2443 //printf("test2\n");
2444 for(Int_t k = 0; k < 540; k++){
2445 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2447 //printf("test3\n");
2450 Float_t devalue = 1.5;
2451 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2452 if(fCalDet) delete fCalDet;
2453 if(fCalDet2) delete fCalDet2;
2454 //printf("test1\n");
2455 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2456 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2457 //printf("test2\n");
2458 for(Int_t k = 0; k < 540; k++){
2459 fCalDet->SetValue(k,devalue);
2460 fCalDet2->SetValue(k,devalue2);
2462 //printf("test3\n");
2467 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2468 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2471 // Init the current detector where we are fCountDet and the
2472 // next fCount for the functions Fit...
2475 // Loop on the Xbins of ch!!
2476 fCountDet = -1; // Current detector
2477 fCount = 0; // To find the next detector
2480 if (fDebugLevel >= 3) {
2481 // Set countdet to the detector
2482 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2483 // Set counter to write at the end of the detector
2485 // Get the right calib objects
2488 if(fDebugLevel == 1) {
2490 fCalibraMode->CalculXBins(fCountDet,i);
2491 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2493 fCalibraMode->CalculXBins(fCountDet,i);
2495 fCount = fCalibraMode->GetXbins(i);
2497 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2498 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2499 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2500 ,(Int_t) GetStack(fCountDet)
2501 ,(Int_t) GetSector(fCountDet),i);
2504 //_______________________________________________________________________________
2505 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2508 // Calculate the number of bins expected (calibration groups)
2511 fNumberOfBinsExpected = 0;
2513 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2514 fNumberOfBinsExpected = 1;
2518 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2519 fNumberOfBinsExpected = 18;
2523 fCalibraMode->ModePadCalibration(2,i);
2524 fCalibraMode->ModePadFragmentation(0,2,0,i);
2525 fCalibraMode->SetDetChamb2(i);
2526 if (fDebugLevel > 1) {
2527 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2529 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2530 fCalibraMode->ModePadCalibration(0,i);
2531 fCalibraMode->ModePadFragmentation(0,0,0,i);
2532 fCalibraMode->SetDetChamb0(i);
2533 if (fDebugLevel > 1) {
2534 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2536 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2539 //_______________________________________________________________________________
2540 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2543 // Calculate the range of fits
2548 if (fDebugLevel == 1) {
2552 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2554 fDect2 = fNumberOfBinsExpected;
2556 if (fDebugLevel >= 3) {
2557 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2558 fCalibraMode->CalculXBins(fCountDet,i);
2559 fDect1 = fCalibraMode->GetXbins(i);
2560 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2561 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2562 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2563 ,(Int_t) GetStack(fCountDet)
2564 ,(Int_t) GetSector(fCountDet),i);
2565 // Set for the next detector
2566 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2569 //_______________________________________________________________________________
2570 Bool_t AliTRDCalibraFit::CheckFitVoir()
2573 // Check if fFitVoir is in the range
2576 if (fFitVoir < fNumberOfBinsExpected) {
2577 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2580 AliInfo("fFitVoir is out of range of the histo!");
2585 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2586 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2589 // See if we are in a new detector and update the
2590 // variables fNfragZ and fNfragRphi if yes
2591 // Will never happen for only one detector (3 and 4)
2592 // Doesn't matter for 2
2594 if (fCount == idect) {
2595 // On en est au detector (or first detector in the group)
2597 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2598 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2599 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2600 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2601 ,(Int_t) GetStack(fCountDet)
2602 ,(Int_t) GetSector(fCountDet),i);
2603 // Set for the next detector
2604 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2609 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2610 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2613 // Reconstruct the min pad row, max pad row, min pad col and
2614 // max pad col of the calibration group for the Fit functions
2615 // idect is the calibration group inside the detector
2617 if (fDebugLevel != 1) {
2618 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2620 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2621 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2623 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2624 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2627 // For the case where there are not enough entries in the histograms
2628 // of the calibration group, the value present in the choosen database
2629 // will be put. A negativ sign enables to know that a fit was not possible.
2632 if (fDebugLevel == 1) {
2633 AliInfo("The element has not enough statistic to be fitted");
2635 else if (fNbDet > 0){
2636 Int_t firstdetector = fCountDet;
2637 Int_t lastdetector = fCountDet+fNbDet;
2638 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2639 ,idect,firstdetector,lastdetector));
2640 // loop over detectors
2641 for(Int_t det = firstdetector; det < lastdetector; det++){
2643 //Set the calibration object again
2647 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2649 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2650 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2651 ,(Int_t) GetStack(fCountDet)
2652 ,(Int_t) GetSector(fCountDet),0);
2653 // Reconstruct row min row max
2654 ReconstructFitRowMinRowMax(idect,0);
2656 // Calcul the coef from the database choosen for the detector
2657 CalculChargeCoefMean(kFALSE);
2659 //stack 2, not stack 2
2661 if(GetStack(fCountDet) == 2) factor = 12;
2664 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2665 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2666 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2667 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2671 //Put default value negative
2672 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2673 fCurrentCoefE = 0.0;
2678 if(fDebugLevel > 1){
2680 if ( !fDebugStreamer ) {
2682 TDirectory *backup = gDirectory;
2683 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2684 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2687 Int_t detector = fCountDet;
2688 Int_t caligroup = idect;
2689 Short_t rowmin = fCalibraMode->GetRowMin(0);
2690 Short_t rowmax = fCalibraMode->GetRowMax(0);
2691 Short_t colmin = fCalibraMode->GetColMin(0);
2692 Short_t colmax = fCalibraMode->GetColMax(0);
2693 Float_t gf = fCurrentCoef[0];
2694 Float_t gfs = fCurrentCoef[1];
2695 Float_t gfE = fCurrentCoefE;
2697 (*fDebugStreamer) << "FillFillCH" <<
2698 "detector=" << detector <<
2699 "caligroup=" << caligroup <<
2700 "rowmin=" << rowmin <<
2701 "rowmax=" << rowmax <<
2702 "colmin=" << colmin <<
2703 "colmax=" << colmax <<
2711 for (Int_t k = 0; k < 2304; k++) {
2712 fCurrentCoefDetector[k] = 0.0;
2716 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2720 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2721 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2723 // Calcul the coef from the database choosen
2724 CalculChargeCoefMean(kFALSE);
2726 //stack 2, not stack 2
2728 if(GetStack(fCountDet) == 2) factor = 12;
2731 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2732 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2733 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2734 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2738 //Put default value negative
2739 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2740 fCurrentCoefE = 0.0;
2749 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2750 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2753 // For the case where there are not enough entries in the histograms
2754 // of the calibration group, the value present in the choosen database
2755 // will be put. A negativ sign enables to know that a fit was not possible.
2757 if (fDebugLevel == 1) {
2758 AliInfo("The element has not enough statistic to be fitted");
2760 else if (fNbDet > 0) {
2762 Int_t firstdetector = fCountDet;
2763 Int_t lastdetector = fCountDet+fNbDet;
2764 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2765 ,idect,firstdetector,lastdetector));
2766 // loop over detectors
2767 for(Int_t det = firstdetector; det < lastdetector; det++){
2769 //Set the calibration object again
2773 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2775 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
2776 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2777 ,(Int_t) GetStack(fCountDet)
2778 ,(Int_t) GetSector(fCountDet),1);
2779 // Reconstruct row min row max
2780 ReconstructFitRowMinRowMax(idect,1);
2782 // Calcul the coef from the database choosen for the detector
2783 CalculVdriftCoefMean();
2786 //stack 2, not stack 2
2788 if(GetStack(fCountDet) == 2) factor = 12;
2791 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2792 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2793 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2794 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2795 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2799 //Put default value negative
2800 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2801 fCurrentCoefE = 0.0;
2802 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2803 fCurrentCoefE2 = 0.0;
2809 if(fDebugLevel > 1){
2811 if ( !fDebugStreamer ) {
2813 TDirectory *backup = gDirectory;
2814 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
2815 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2819 Int_t detector = fCountDet;
2820 Int_t caligroup = idect;
2821 Short_t rowmin = fCalibraMode->GetRowMin(1);
2822 Short_t rowmax = fCalibraMode->GetRowMax(1);
2823 Short_t colmin = fCalibraMode->GetColMin(1);
2824 Short_t colmax = fCalibraMode->GetColMax(1);
2825 Float_t vf = fCurrentCoef[0];
2826 Float_t vs = fCurrentCoef[1];
2827 Float_t vfE = fCurrentCoefE;
2828 Float_t t0f = fCurrentCoef2[0];
2829 Float_t t0s = fCurrentCoef2[1];
2830 Float_t t0E = fCurrentCoefE2;
2834 (* fDebugStreamer) << "FillFillPH"<<
2835 "detector="<<detector<<
2836 "nentries="<<nentries<<
2837 "caligroup="<<caligroup<<
2851 for (Int_t k = 0; k < 2304; k++) {
2852 fCurrentCoefDetector[k] = 0.0;
2853 fCurrentCoefDetector2[k] = 0.0;
2857 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2861 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2862 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2864 CalculVdriftCoefMean();
2867 //stack 2 and not stack 2
2869 if(GetStack(fCountDet) == 2) factor = 12;
2873 // Fill the fCurrentCoefDetector 2
2874 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2875 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2876 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2877 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2881 // Put the default value
2882 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2883 fCurrentCoefE = 0.0;
2884 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2885 fCurrentCoefE2 = 0.0;
2887 FillFillPH(idect,nentries);
2896 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2897 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2900 // For the case where there are not enough entries in the histograms
2901 // of the calibration group, the value present in the choosen database
2902 // will be put. A negativ sign enables to know that a fit was not possible.
2905 if (fDebugLevel == 1) {
2906 AliInfo("The element has not enough statistic to be fitted");
2908 else if (fNbDet > 0){
2910 Int_t firstdetector = fCountDet;
2911 Int_t lastdetector = fCountDet+fNbDet;
2912 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2913 ,idect,firstdetector,lastdetector));
2915 // loop over detectors
2916 for(Int_t det = firstdetector; det < lastdetector; det++){
2918 //Set the calibration object again
2922 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2924 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
2925 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2926 ,(Int_t) GetStack(fCountDet)
2927 ,(Int_t) GetSector(fCountDet),2);
2928 // Reconstruct row min row max
2929 ReconstructFitRowMinRowMax(idect,2);
2931 // Calcul the coef from the database choosen for the detector
2932 CalculPRFCoefMean();
2934 //stack 2, not stack 2
2936 if(GetStack(fCountDet) == 2) factor = 12;
2939 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2940 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2941 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2942 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2946 //Put default value negative
2947 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2948 fCurrentCoefE = 0.0;
2953 if(fDebugLevel > 1){
2955 if ( !fDebugStreamer ) {
2957 TDirectory *backup = gDirectory;
2958 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
2959 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2962 Int_t detector = fCountDet;
2963 Int_t layer = GetLayer(fCountDet);
2964 Int_t caligroup = idect;
2965 Short_t rowmin = fCalibraMode->GetRowMin(2);
2966 Short_t rowmax = fCalibraMode->GetRowMax(2);
2967 Short_t colmin = fCalibraMode->GetColMin(2);
2968 Short_t colmax = fCalibraMode->GetColMax(2);
2969 Float_t widf = fCurrentCoef[0];
2970 Float_t wids = fCurrentCoef[1];
2971 Float_t widfE = fCurrentCoefE;
2973 (* fDebugStreamer) << "FillFillPRF"<<
2974 "detector="<<detector<<
2976 "caligroup="<<caligroup<<
2987 for (Int_t k = 0; k < 2304; k++) {
2988 fCurrentCoefDetector[k] = 0.0;
2992 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2996 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2997 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2999 CalculPRFCoefMean();
3001 // stack 2 and not stack 2
3003 if(GetStack(fCountDet) == 2) factor = 12;
3007 // Fill the fCurrentCoefDetector
3008 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3009 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3010 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3014 // Put the default value
3015 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3016 fCurrentCoefE = 0.0;
3024 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3025 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3028 // For the case where there are not enough entries in the histograms
3029 // of the calibration group, the value present in the choosen database
3030 // will be put. A negativ sign enables to know that a fit was not possible.
3033 // Calcul the coef from the database choosen
3034 CalculVdriftLorentzCoef();
3037 if(GetStack(fCountDet) == 2) factor = 1728;
3041 // Fill the fCurrentCoefDetector
3042 for (Int_t k = 0; k < factor; k++) {
3043 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3044 // should be negative
3045 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3049 //Put default opposite sign
3050 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3051 fCurrentCoefE = 0.0;
3052 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3053 fCurrentCoefE2 = 0.0;
3055 FillFillLinearFitter();
3060 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3061 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3064 // Fill the coefficients found with the fits or other
3065 // methods from the Fit functions
3068 if (fDebugLevel != 1) {
3070 Int_t firstdetector = fCountDet;
3071 Int_t lastdetector = fCountDet+fNbDet;
3072 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3073 ,idect,firstdetector,lastdetector));
3074 // loop over detectors
3075 for(Int_t det = firstdetector; det < lastdetector; det++){
3077 //Set the calibration object again
3081 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3083 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3084 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3085 ,(Int_t) GetStack(fCountDet)
3086 ,(Int_t) GetSector(fCountDet),0);
3087 // Reconstruct row min row max
3088 ReconstructFitRowMinRowMax(idect,0);
3090 // Calcul the coef from the database choosen for the detector
3091 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3092 else CalculChargeCoefMean(kTRUE);
3094 //stack 2, not stack 2
3096 if(GetStack(fCountDet) == 2) factor = 12;
3099 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3100 Double_t coeftoput = 1.0;
3101 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3102 else coeftoput = fCurrentCoef[0];
3103 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3104 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3105 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3112 if(fDebugLevel > 1){
3114 if ( !fDebugStreamer ) {
3116 TDirectory *backup = gDirectory;
3117 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3118 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3121 Int_t detector = fCountDet;
3122 Int_t caligroup = idect;
3123 Short_t rowmin = fCalibraMode->GetRowMin(0);
3124 Short_t rowmax = fCalibraMode->GetRowMax(0);
3125 Short_t colmin = fCalibraMode->GetColMin(0);
3126 Short_t colmax = fCalibraMode->GetColMax(0);
3127 Float_t gf = fCurrentCoef[0];
3128 Float_t gfs = fCurrentCoef[1];
3129 Float_t gfE = fCurrentCoefE;
3131 (*fDebugStreamer) << "FillFillCH" <<
3132 "detector=" << detector <<
3133 "caligroup=" << caligroup <<
3134 "rowmin=" << rowmin <<
3135 "rowmax=" << rowmax <<
3136 "colmin=" << colmin <<
3137 "colmax=" << colmax <<
3145 for (Int_t k = 0; k < 2304; k++) {
3146 fCurrentCoefDetector[k] = 0.0;
3150 //printf("Check the count now: fCountDet %d\n",fCountDet);
3155 if(GetStack(fCountDet) == 2) factor = 12;
3158 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3159 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3160 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3171 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3172 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3175 // Fill the coefficients found with the fits or other
3176 // methods from the Fit functions
3179 if (fDebugLevel != 1) {
3182 Int_t firstdetector = fCountDet;
3183 Int_t lastdetector = fCountDet+fNbDet;
3184 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3185 ,idect,firstdetector,lastdetector));
3187 // loop over detectors
3188 for(Int_t det = firstdetector; det < lastdetector; det++){
3190 //Set the calibration object again
3194 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3196 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3197 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3198 ,(Int_t) GetStack(fCountDet)
3199 ,(Int_t) GetSector(fCountDet),1);
3200 // Reconstruct row min row max
3201 ReconstructFitRowMinRowMax(idect,1);
3203 // Calcul the coef from the database choosen for the detector
3204 CalculVdriftCoefMean();
3207 //stack 2, not stack 2
3209 if(GetStack(fCountDet) == 2) factor = 12;
3212 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3213 Double_t coeftoput = 1.5;
3214 Double_t coeftoput2 = 0.0;
3216 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3217 else coeftoput = fCurrentCoef[0];
3219 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3220 else coeftoput2 = fCurrentCoef2[0];
3222 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3223 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3224 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3225 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3233 if(fDebugLevel > 1){
3235 if ( !fDebugStreamer ) {
3237 TDirectory *backup = gDirectory;
3238 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3239 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3243 Int_t detector = fCountDet;
3244 Int_t caligroup = idect;
3245 Short_t rowmin = fCalibraMode->GetRowMin(1);
3246 Short_t rowmax = fCalibraMode->GetRowMax(1);
3247 Short_t colmin = fCalibraMode->GetColMin(1);
3248 Short_t colmax = fCalibraMode->GetColMax(1);
3249 Float_t vf = fCurrentCoef[0];
3250 Float_t vs = fCurrentCoef[1];
3251 Float_t vfE = fCurrentCoefE;
3252 Float_t t0f = fCurrentCoef2[0];
3253 Float_t t0s = fCurrentCoef2[1];
3254 Float_t t0E = fCurrentCoefE2;
3258 (* fDebugStreamer) << "FillFillPH"<<
3259 "detector="<<detector<<
3260 "nentries="<<nentries<<
3261 "caligroup="<<caligroup<<
3275 for (Int_t k = 0; k < 2304; k++) {
3276 fCurrentCoefDetector[k] = 0.0;
3277 fCurrentCoefDetector2[k] = 0.0;
3281 //printf("Check the count now: fCountDet %d\n",fCountDet);
3286 if(GetStack(fCountDet) == 2) factor = 12;
3289 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3290 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3291 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3292 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3296 FillFillPH(idect,nentries);
3301 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3302 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3305 // Fill the coefficients found with the fits or other
3306 // methods from the Fit functions
3309 if (fDebugLevel != 1) {
3312 Int_t firstdetector = fCountDet;
3313 Int_t lastdetector = fCountDet+fNbDet;
3314 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3315 ,idect,firstdetector,lastdetector));
3317 // loop over detectors
3318 for(Int_t det = firstdetector; det < lastdetector; det++){
3320 //Set the calibration object again
3324 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3326 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3327 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3328 ,(Int_t) GetStack(fCountDet)
3329 ,(Int_t) GetSector(fCountDet),2);
3330 // Reconstruct row min row max
3331 ReconstructFitRowMinRowMax(idect,2);
3333 // Calcul the coef from the database choosen for the detector
3334 CalculPRFCoefMean();
3336 //stack 2, not stack 2
3338 if(GetStack(fCountDet) == 2) factor = 12;
3341 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3342 Double_t coeftoput = 1.0;
3343 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3344 else coeftoput = fCurrentCoef[0];
3345 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3346 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3347 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3354 if(fDebugLevel > 1){
3356 if ( !fDebugStreamer ) {
3358 TDirectory *backup = gDirectory;
3359 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3360 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3363 Int_t detector = fCountDet;
3364 Int_t layer = GetLayer(fCountDet);
3365 Int_t caligroup = idect;
3366 Short_t rowmin = fCalibraMode->GetRowMin(2);
3367 Short_t rowmax = fCalibraMode->GetRowMax(2);
3368 Short_t colmin = fCalibraMode->GetColMin(2);
3369 Short_t colmax = fCalibraMode->GetColMax(2);
3370 Float_t widf = fCurrentCoef[0];
3371 Float_t wids = fCurrentCoef[1];
3372 Float_t widfE = fCurrentCoefE;
3374 (* fDebugStreamer) << "FillFillPRF"<<
3375 "detector="<<detector<<
3377 "caligroup="<<caligroup<<
3388 for (Int_t k = 0; k < 2304; k++) {
3389 fCurrentCoefDetector[k] = 0.0;
3393 //printf("Check the count now: fCountDet %d\n",fCountDet);
3398 if(GetStack(fCountDet) == 2) factor = 12;
3401 // Pointer to the branch
3402 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3403 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3404 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3414 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3415 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3418 // Fill the coefficients found with the fits or other
3419 // methods from the Fit functions
3423 if(GetStack(fCountDet) == 2) factor = 1728;
3426 // Pointer to the branch
3427 for (Int_t k = 0; k < factor; k++) {
3428 fCurrentCoefDetector[k] = fCurrentCoef[0];
3429 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3432 FillFillLinearFitter();
3437 //________________________________________________________________________________
3438 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3441 // DebugStream and fVectorFit
3444 // End of one detector
3445 if ((idect == (fCount-1))) {
3448 for (Int_t k = 0; k < 2304; k++) {
3449 fCurrentCoefDetector[k] = 0.0;
3453 if(fDebugLevel > 1){
3455 if ( !fDebugStreamer ) {
3457 TDirectory *backup = gDirectory;
3458 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3459 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3462 Int_t detector = fCountDet;
3463 Int_t caligroup = idect;
3464 Short_t rowmin = fCalibraMode->GetRowMin(0);
3465 Short_t rowmax = fCalibraMode->GetRowMax(0);
3466 Short_t colmin = fCalibraMode->GetColMin(0);
3467 Short_t colmax = fCalibraMode->GetColMax(0);
3468 Float_t gf = fCurrentCoef[0];
3469 Float_t gfs = fCurrentCoef[1];
3470 Float_t gfE = fCurrentCoefE;
3472 (*fDebugStreamer) << "FillFillCH" <<
3473 "detector=" << detector <<
3474 "caligroup=" << caligroup <<
3475 "rowmin=" << rowmin <<
3476 "rowmax=" << rowmax <<
3477 "colmin=" << colmin <<
3478 "colmax=" << colmax <<
3486 //________________________________________________________________________________
3487 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3490 // DebugStream and fVectorFit and fVectorFit2
3493 // End of one detector
3494 if ((idect == (fCount-1))) {
3498 for (Int_t k = 0; k < 2304; k++) {
3499 fCurrentCoefDetector[k] = 0.0;
3500 fCurrentCoefDetector2[k] = 0.0;
3504 if(fDebugLevel > 1){
3506 if ( !fDebugStreamer ) {
3508 TDirectory *backup = gDirectory;
3509 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3510 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3514 Int_t detector = fCountDet;
3515 Int_t caligroup = idect;
3516 Short_t rowmin = fCalibraMode->GetRowMin(1);
3517 Short_t rowmax = fCalibraMode->GetRowMax(1);
3518 Short_t colmin = fCalibraMode->GetColMin(1);
3519 Short_t colmax = fCalibraMode->GetColMax(1);
3520 Float_t vf = fCurrentCoef[0];
3521 Float_t vs = fCurrentCoef[1];
3522 Float_t vfE = fCurrentCoefE;
3523 Float_t t0f = fCurrentCoef2[0];
3524 Float_t t0s = fCurrentCoef2[1];
3525 Float_t t0E = fCurrentCoefE2;
3529 (* fDebugStreamer) << "FillFillPH"<<
3530 "detector="<<detector<<
3531 "nentries="<<nentries<<
3532 "caligroup="<<caligroup<<
3547 //________________________________________________________________________________
3548 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3551 // DebugStream and fVectorFit
3554 // End of one detector
3555 if ((idect == (fCount-1))) {
3558 for (Int_t k = 0; k < 2304; k++) {
3559 fCurrentCoefDetector[k] = 0.0;
3564 if(fDebugLevel > 1){
3566 if ( !fDebugStreamer ) {
3568 TDirectory *backup = gDirectory;
3569 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3570 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3573 Int_t detector = fCountDet;
3574 Int_t layer = GetLayer(fCountDet);
3575 Int_t caligroup = idect;
3576 Short_t rowmin = fCalibraMode->GetRowMin(2);
3577 Short_t rowmax = fCalibraMode->GetRowMax(2);
3578 Short_t colmin = fCalibraMode->GetColMin(2);
3579 Short_t colmax = fCalibraMode->GetColMax(2);
3580 Float_t widf = fCurrentCoef[0];
3581 Float_t wids = fCurrentCoef[1];
3582 Float_t widfE = fCurrentCoefE;
3584 (* fDebugStreamer) << "FillFillPRF"<<
3585 "detector="<<detector<<
3587 "caligroup="<<caligroup<<
3599 //________________________________________________________________________________
3600 void AliTRDCalibraFit::FillFillLinearFitter()
3603 // DebugStream and fVectorFit
3606 // End of one detector
3612 for (Int_t k = 0; k < 2304; k++) {
3613 fCurrentCoefDetector[k] = 0.0;
3614 fCurrentCoefDetector2[k] = 0.0;
3618 if(fDebugLevel > 1){
3620 if ( !fDebugStreamer ) {
3622 TDirectory *backup = gDirectory;
3623 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3624 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3627 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3628 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3629 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3630 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3631 Float_t tiltangle = padplane->GetTiltingAngle();
3632 Int_t detector = fCountDet;
3633 Int_t stack = GetStack(fCountDet);
3634 Int_t layer = GetLayer(fCountDet);
3635 Float_t vf = fCurrentCoef[0];
3636 Float_t vs = fCurrentCoef[1];
3637 Float_t vfE = fCurrentCoefE;
3638 Float_t lorentzangler = fCurrentCoef2[0];
3639 Float_t elorentzangler = fCurrentCoefE2;
3640 Float_t lorentzangles = fCurrentCoef2[1];
3642 (* fDebugStreamer) << "FillFillLinearFitter"<<
3643 "detector="<<detector<<
3648 "tiltangle="<<tiltangle<<
3652 "lorentzangler="<<lorentzangler<<
3653 "Elorentzangler="<<elorentzangler<<
3654 "lorentzangles="<<lorentzangles<<
3660 //____________Calcul Coef Mean_________________________________________________
3662 //_____________________________________________________________________________
3663 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3666 // For the detector Dect calcul the mean time 0
3667 // for the calibration group idect from the choosen database
3670 fCurrentCoef2[1] = 0.0;
3671 if(fDebugLevel != 1){
3672 if(((fCalibraMode->GetNz(1) > 0) ||
3673 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3675 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3676 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3677 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3681 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3687 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3691 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3692 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3693 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3696 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3704 //_____________________________________________________________________________
3705 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3708 // For the detector Dect calcul the mean gain factor
3709 // for the calibration group idect from the choosen database
3712 fCurrentCoef[1] = 0.0;
3713 if(fDebugLevel != 1){
3714 if (((fCalibraMode->GetNz(0) > 0) ||
3715 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3716 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3717 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3718 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3719 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3722 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3726 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3727 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3732 //_____________________________________________________________________________
3733 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3736 // For the detector Dect calcul the mean sigma of pad response
3737 // function for the calibration group idect from the choosen database
3740 fCurrentCoef[1] = 0.0;
3741 if(fDebugLevel != 1){
3742 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3743 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3744 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3747 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3751 //_____________________________________________________________________________
3752 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3755 // For the detector dect calcul the mean drift velocity for the
3756 // calibration group idect from the choosen database
3759 fCurrentCoef[1] = 0.0;
3760 if(fDebugLevel != 1){
3761 if (((fCalibraMode->GetNz(1) > 0) ||
3762 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3764 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3765 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3766 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3770 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3775 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3780 //_____________________________________________________________________________
3781 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
3784 // For the detector fCountDet, mean drift velocity and tan lorentzangle
3787 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
3788 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3792 //_____________________________________________________________________________
3793 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
3796 // Default width of the PRF if there is no database as reference
3801 //case 0: return 0.515;
3802 //case 1: return 0.502;
3803 //case 2: return 0.491;
3804 //case 3: return 0.481;
3805 //case 4: return 0.471;
3806 //case 5: return 0.463;
3807 //default: return 0.0;
3810 case 0: return 0.538429;
3811 case 1: return 0.524302;
3812 case 2: return 0.511591;
3813 case 3: return 0.500140;
3814 case 4: return 0.489821;
3815 case 5: return 0.480524;
3816 default: return 0.0;
3819 //________________________________________________________________________________
3820 void AliTRDCalibraFit::SetCalROC(Int_t i)
3823 // Set the calib object for fCountDet
3826 Float_t value = 0.0;
3828 //Get the CalDet object
3830 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
3832 AliInfo("Could not get calibDB");
3839 fCalROC->~AliTRDCalROC();
3840 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3841 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3845 fCalROC->~AliTRDCalROC();
3846 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3847 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3849 fCalROC2->~AliTRDCalROC();
3850 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3851 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3855 fCalROC->~AliTRDCalROC();
3856 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3857 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3866 if(fCalROC) delete fCalROC;
3867 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3868 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3869 fCalROC->SetValue(k,1.0);
3873 if(fCalROC) delete fCalROC;
3874 if(fCalROC2) delete fCalROC2;
3875 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3876 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3877 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3878 fCalROC->SetValue(k,1.0);
3879 fCalROC2->SetValue(k,0.0);
3883 if(fCalROC) delete fCalROC;
3884 value = GetPRFDefault(GetLayer(fCountDet));
3885 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3886 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3887 fCalROC->SetValue(k,value);
3895 //____________Fit Methods______________________________________________________
3897 //_____________________________________________________________________________
3898 void AliTRDCalibraFit::FitPente(TH1* projPH)
3901 // Slope methode for the drift velocity
3905 const Float_t kDrWidth = AliTRDgeometry::DrThick();
3912 fCurrentCoefE = 0.0;
3913 fCurrentCoefE2 = 0.0;
3914 fCurrentCoef[0] = 0.0;
3915 fCurrentCoef2[0] = 0.0;
3916 TLine *line = new TLine();
3919 TAxis *xpph = projPH->GetXaxis();
3920 Int_t nbins = xpph->GetNbins();
3921 Double_t lowedge = xpph->GetBinLowEdge(1);
3922 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3923 Double_t widbins = (upedge - lowedge) / nbins;
3924 Double_t limit = upedge + 0.5 * widbins;
3927 // Beginning of the signal
3928 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3929 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
3930 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3932 binmax = (Int_t) pentea->GetMaximumBin();
3935 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3937 if (binmax >= nbins) {
3940 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3942 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3943 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3944 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3945 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3946 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3948 fPhd[0] = -(l3P1am / (2 * l3P2am));
3951 if((l3P1am != 0.0) && (l3P2am != 0.0)){
3952 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3955 // Amplification region
3958 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3959 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3966 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3968 if (binmax >= nbins) {
3971 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3973 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3974 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3975 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3976 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3977 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3979 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3981 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3982 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3985 fCurrentCoefE2 = fCurrentCoefE;
3988 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3989 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
3990 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3993 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3996 AliInfo("Put the binmax from 1 to 2 to enable the fit");
3998 if (binmin >= nbins) {
4001 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4006 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4007 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4008 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4009 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4010 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4011 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4013 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4015 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
4016 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4018 Float_t fPhdt0 = 0.0;
4019 Float_t t0Shift = 0.0;
4022 t0Shift = fT0Shift1;
4026 t0Shift = fT0Shift0;
4029 if ((fPhd[2] > fPhd[0]) &&
4030 (fPhd[2] > fPhd[1]) &&
4031 (fPhd[1] > fPhd[0]) &&
4033 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4034 fNumberFitSuccess++;
4036 if (fPhdt0 >= 0.0) {
4037 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4038 if (fCurrentCoef2[0] < -1.0) {
4039 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4043 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4048 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4049 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4052 if (fDebugLevel == 1) {
4053 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4056 line->SetLineColor(2);
4057 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4058 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4059 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4060 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4061 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4062 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4063 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4064 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4067 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4076 //_____________________________________________________________________________
4077 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4080 // Slope methode but with polynomes de Lagrange
4084 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4087 Double_t *x = new Double_t[5];
4088 Double_t *y = new Double_t[5];
4103 fCurrentCoefE = 0.0;
4104 fCurrentCoefE2 = 1.0;
4105 fCurrentCoef[0] = 0.0;
4106 fCurrentCoef2[0] = 0.0;
4107 TLine *line = new TLine();
4108 TF1 * polynome = 0x0;
4109 TF1 * polynomea = 0x0;
4110 TF1 * polynomeb = 0x0;
4114 TAxis *xpph = projPH->GetXaxis();
4115 Int_t nbins = xpph->GetNbins();
4116 Double_t lowedge = xpph->GetBinLowEdge(1);
4117 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4118 Double_t widbins = (upedge - lowedge) / nbins;
4119 Double_t limit = upedge + 0.5 * widbins;
4124 // Beginning of the signal
4125 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4126 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4127 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4130 binmax = (Int_t) pentea->GetMaximumBin();
4132 Double_t minnn = 0.0;
4133 Double_t maxxx = 0.0;
4135 Int_t kase = nbins-binmax;
4143 minnn = pentea->GetBinCenter(binmax-2);
4144 maxxx = pentea->GetBinCenter(binmax);
4145 x[0] = pentea->GetBinCenter(binmax-2);
4146 x[1] = pentea->GetBinCenter(binmax-1);
4147 x[2] = pentea->GetBinCenter(binmax);
4148 y[0] = pentea->GetBinContent(binmax-2);
4149 y[1] = pentea->GetBinContent(binmax-1);
4150 y[2] = pentea->GetBinContent(binmax);
4151 c = CalculPolynomeLagrange2(x,y);
4152 AliInfo("At the limit for beginning!");
4155 minnn = pentea->GetBinCenter(binmax-2);
4156 maxxx = pentea->GetBinCenter(binmax+1);
4157 x[0] = pentea->GetBinCenter(binmax-2);
4158 x[1] = pentea->GetBinCenter(binmax-1);
4159 x[2] = pentea->GetBinCenter(binmax);
4160 x[3] = pentea->GetBinCenter(binmax+1);
4161 y[0] = pentea->GetBinContent(binmax-2);
4162 y[1] = pentea->GetBinContent(binmax-1);
4163 y[2] = pentea->GetBinContent(binmax);
4164 y[3] = pentea->GetBinContent(binmax+1);
4165 c = CalculPolynomeLagrange3(x,y);
4173 minnn = pentea->GetBinCenter(binmax);
4174 maxxx = pentea->GetBinCenter(binmax+2);
4175 x[0] = pentea->GetBinCenter(binmax);
4176 x[1] = pentea->GetBinCenter(binmax+1);
4177 x[2] = pentea->GetBinCenter(binmax+2);
4178 y[0] = pentea->GetBinContent(binmax);
4179 y[1] = pentea->GetBinContent(binmax+1);
4180 y[2] = pentea->GetBinContent(binmax+2);
4181 c = CalculPolynomeLagrange2(x,y);
4184 minnn = pentea->GetBinCenter(binmax-1);
4185 maxxx = pentea->GetBinCenter(binmax+2);
4186 x[0] = pentea->GetBinCenter(binmax-1);
4187 x[1] = pentea->GetBinCenter(binmax);
4188 x[2] = pentea->GetBinCenter(binmax+1);
4189 x[3] = pentea->GetBinCenter(binmax+2);
4190 y[0] = pentea->GetBinContent(binmax-1);
4191 y[1] = pentea->GetBinContent(binmax);
4192 y[2] = pentea->GetBinContent(binmax+1);
4193 y[3] = pentea->GetBinContent(binmax+2);
4194 c = CalculPolynomeLagrange3(x,y);
4197 minnn = pentea->GetBinCenter(binmax-2);
4198 maxxx = pentea->GetBinCenter(binmax+2);
4199 x[0] = pentea->GetBinCenter(binmax-2);
4200 x[1] = pentea->GetBinCenter(binmax-1);
4201 x[2] = pentea->GetBinCenter(binmax);
4202 x[3] = pentea->GetBinCenter(binmax+1);
4203 x[4] = pentea->GetBinCenter(binmax+2);
4204 y[0] = pentea->GetBinContent(binmax-2);
4205 y[1] = pentea->GetBinContent(binmax-1);
4206 y[2] = pentea->GetBinContent(binmax);
4207 y[3] = pentea->GetBinContent(binmax+1);
4208 y[4] = pentea->GetBinContent(binmax+2);
4209 c = CalculPolynomeLagrange4(x,y);
4217 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4218 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4220 Double_t step = (maxxx-minnn)/10000;
4222 Double_t maxvalue = 0.0;
4223 Double_t placemaximum = minnn;
4224 for(Int_t o = 0; o < 10000; o++){
4225 if(o == 0) maxvalue = polynomeb->Eval(l);
4226 if(maxvalue < (polynomeb->Eval(l))){
4227 maxvalue = polynomeb->Eval(l);
4232 fPhd[0] = placemaximum;
4235 // Amplification region
4238 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4239 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4245 Double_t minn = 0.0;
4246 Double_t maxx = 0.0;
4258 Int_t kase1 = nbins - binmax;
4260 //Determination of minn and maxx
4261 //case binmax = nbins
4266 minn = projPH->GetBinCenter(binmax-2);
4267 maxx = projPH->GetBinCenter(binmax);
4268 x[0] = projPH->GetBinCenter(binmax-2);
4269 x[1] = projPH->GetBinCenter(binmax-1);
4270 x[2] = projPH->GetBinCenter(binmax);
4271 y[0] = projPH->GetBinContent(binmax-2);
4272 y[1] = projPH->GetBinContent(binmax-1);
4273 y[2] = projPH->GetBinContent(binmax);
4274 c = CalculPolynomeLagrange2(x,y);
4275 //AliInfo("At the limit for the drift!");
4278 minn = projPH->GetBinCenter(binmax-2);
4279 maxx = projPH->GetBinCenter(binmax+1);
4280 x[0] = projPH->GetBinCenter(binmax-2);
4281 x[1] = projPH->GetBinCenter(binmax-1);
4282 x[2] = projPH->GetBinCenter(binmax);
4283 x[3] = projPH->GetBinCenter(binmax+1);
4284 y[0] = projPH->GetBinContent(binmax-2);
4285 y[1] = projPH->GetBinContent(binmax-1);
4286 y[2] = projPH->GetBinContent(binmax);
4287 y[3] = projPH->GetBinContent(binmax+1);
4288 c = CalculPolynomeLagrange3(x,y);
4297 minn = projPH->GetBinCenter(binmax);
4298 maxx = projPH->GetBinCenter(binmax+2);
4299 x[0] = projPH->GetBinCenter(binmax);
4300 x[1] = projPH->GetBinCenter(binmax+1);
4301 x[2] = projPH->GetBinCenter(binmax+2);
4302 y[0] = projPH->GetBinContent(binmax);
4303 y[1] = projPH->GetBinContent(binmax+1);
4304 y[2] = projPH->GetBinContent(binmax+2);
4305 c = CalculPolynomeLagrange2(x,y);
4308 minn = projPH->GetBinCenter(binmax-1);
4309 maxx = projPH->GetBinCenter(binmax+2);
4310 x[0] = projPH->GetBinCenter(binmax-1);
4311 x[1] = projPH->GetBinCenter(binmax);
4312 x[2] = projPH->GetBinCenter(binmax+1);
4313 x[3] = projPH->GetBinCenter(binmax+2);
4314 y[0] = projPH->GetBinContent(binmax-1);
4315 y[1] = projPH->GetBinContent(binmax);
4316 y[2] = projPH->GetBinContent(binmax+1);
4317 y[3] = projPH->GetBinContent(binmax+2);
4318 c = CalculPolynomeLagrange3(x,y);
4321 minn = projPH->GetBinCenter(binmax-2);
4322 maxx = projPH->GetBinCenter(binmax+2);
4323 x[0] = projPH->GetBinCenter(binmax-2);
4324 x[1] = projPH->GetBinCenter(binmax-1);
4325 x[2] = projPH->GetBinCenter(binmax);
4326 x[3] = projPH->GetBinCenter(binmax+1);
4327 x[4] = projPH->GetBinCenter(binmax+2);
4328 y[0] = projPH->GetBinContent(binmax-2);
4329 y[1] = projPH->GetBinContent(binmax-1);
4330 y[2] = projPH->GetBinContent(binmax);
4331 y[3] = projPH->GetBinContent(binmax+1);
4332 y[4] = projPH->GetBinContent(binmax+2);
4333 c = CalculPolynomeLagrange4(x,y);
4340 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4341 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4343 Double_t step = (maxx-minn)/1000;
4345 Double_t maxvalue = 0.0;
4346 Double_t placemaximum = minn;
4347 for(Int_t o = 0; o < 1000; o++){
4348 if(o == 0) maxvalue = polynomea->Eval(l);
4349 if(maxvalue < (polynomea->Eval(l))){
4350 maxvalue = polynomea->Eval(l);
4355 fPhd[1] = placemaximum;
4359 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4360 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4361 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4364 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4370 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4374 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4375 AliInfo("Too many fluctuations at the end!");
4378 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4379 AliInfo("Too many fluctuations at the end!");
4382 if(pente->GetBinContent(binmin+1)==0){
4383 AliInfo("No entries for the next bin!");
4384 pente->SetBinContent(binmin,0);
4385 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4401 Bool_t case1 = kFALSE;
4402 Bool_t case2 = kFALSE;
4403 Bool_t case4 = kFALSE;
4405 //Determination of min and max
4406 //case binmin <= nbins-3
4408 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4409 min = pente->GetBinCenter(binmin-2);
4410 max = pente->GetBinCenter(binmin+2);
4411 x[0] = pente->GetBinCenter(binmin-2);
4412 x[1] = pente->GetBinCenter(binmin-1);
4413 x[2] = pente->GetBinCenter(binmin);
4414 x[3] = pente->GetBinCenter(binmin+1);
4415 x[4] = pente->GetBinCenter(binmin+2);
4416 y[0] = pente->GetBinContent(binmin-2);
4417 y[1] = pente->GetBinContent(binmin-1);
4418 y[2] = pente->GetBinContent(binmin);
4419 y[3] = pente->GetBinContent(binmin+1);
4420 y[4] = pente->GetBinContent(binmin+2);
4421 //Calcul the polynome de Lagrange
4422 c = CalculPolynomeLagrange4(x,y);
4424 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4425 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4426 //AliInfo("polynome 4 false 1");
4429 if(((binmin+3) <= (nbins-1)) &&
4430 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4431 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4432 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4433 AliInfo("polynome 4 false 2");
4437 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4438 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4439 //AliInfo("polynome 4 case 1");
4442 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4443 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4444 //AliInfo("polynome 4 case 4");
4449 //case binmin = nbins-2
4451 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4453 min = pente->GetBinCenter(binmin-2);
4454 max = pente->GetBinCenter(binmin+1);
4455 x[0] = pente->GetBinCenter(binmin-2);
4456 x[1] = pente->GetBinCenter(binmin-1);
4457 x[2] = pente->GetBinCenter(binmin);
4458 x[3] = pente->GetBinCenter(binmin+1);
4459 y[0] = pente->GetBinContent(binmin-2);
4460 y[1] = pente->GetBinContent(binmin-1);
4461 y[2] = pente->GetBinContent(binmin);
4462 y[3] = pente->GetBinContent(binmin+1);
4463 //Calcul the polynome de Lagrange
4464 c = CalculPolynomeLagrange3(x,y);
4465 //richtung +: nothing
4467 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4468 //AliInfo("polynome 3- case 2");
4473 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4475 min = pente->GetBinCenter(binmin-1);
4476 max = pente->GetBinCenter(binmin+2);
4477 x[0] = pente->GetBinCenter(binmin-1);
4478 x[1] = pente->GetBinCenter(binmin);
4479 x[2] = pente->GetBinCenter(binmin+1);
4480 x[3] = pente->GetBinCenter(binmin+2);
4481 y[0] = pente->GetBinContent(binmin-1);
4482 y[1] = pente->GetBinContent(binmin);
4483 y[2] = pente->GetBinContent(binmin+1);
4484 y[3] = pente->GetBinContent(binmin+2);
4485 //Calcul the polynome de Lagrange
4486 c = CalculPolynomeLagrange3(x,y);
4488 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4489 //AliInfo("polynome 3+ case 2");
4494 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4495 min = pente->GetBinCenter(binmin);
4496 max = pente->GetBinCenter(binmin+2);
4497 x[0] = pente->GetBinCenter(binmin);
4498 x[1] = pente->GetBinCenter(binmin+1);
4499 x[2] = pente->GetBinCenter(binmin+2);
4500 y[0] = pente->GetBinContent(binmin);
4501 y[1] = pente->GetBinContent(binmin+1);
4502 y[2] = pente->GetBinContent(binmin+2);
4503 //Calcul the polynome de Lagrange
4504 c = CalculPolynomeLagrange2(x,y);
4506 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4507 //AliInfo("polynome 2+ false");
4512 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4514 min = pente->GetBinCenter(binmin-1);
4515 max = pente->GetBinCenter(binmin+1);
4516 x[0] = pente->GetBinCenter(binmin-1);
4517 x[1] = pente->GetBinCenter(binmin);
4518 x[2] = pente->GetBinCenter(binmin+1);
4519 y[0] = pente->GetBinContent(binmin-1);
4520 y[1] = pente->GetBinContent(binmin);
4521 y[2] = pente->GetBinContent(binmin+1);
4522 //Calcul the polynome de Lagrange
4523 c = CalculPolynomeLagrange2(x,y);
4524 //richtung +: nothing
4525 //richtung -: nothing
4527 //case binmin = nbins-1
4529 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4530 min = pente->GetBinCenter(binmin-2);
4531 max = pente->GetBinCenter(binmin);
4532 x[0] = pente->GetBinCenter(binmin-2);
4533 x[1] = pente->GetBinCenter(binmin-1);
4534 x[2] = pente->GetBinCenter(binmin);
4535 y[0] = pente->GetBinContent(binmin-2);
4536 y[1] = pente->GetBinContent(binmin-1);
4537 y[2] = pente->GetBinContent(binmin);
4538 //Calcul the polynome de Lagrange
4539 c = CalculPolynomeLagrange2(x,y);
4540 //AliInfo("At the limit for the drift!");
4541 //fluctuation too big!
4542 //richtung +: nothing
4544 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4545 //AliInfo("polynome 2- false ");
4549 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4551 AliInfo("At the limit for the drift and not usable!");
4555 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4557 AliInfo("For the drift...problem!");
4559 //pass but should not happen
4560 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4562 AliInfo("For the drift...problem!");
4566 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4567 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4568 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4569 Double_t step = (max-min)/1000;
4571 Double_t minvalue = 0.0;
4572 Double_t placeminimum = min;
4573 for(Int_t o = 0; o < 1000; o++){
4574 if(o == 0) minvalue = polynome->Eval(l);
4575 if(minvalue > (polynome->Eval(l))){
4576 minvalue = polynome->Eval(l);
4581 fPhd[2] = placeminimum;
4583 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4584 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4585 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4587 Float_t fPhdt0 = 0.0;
4588 Float_t t0Shift = 0.0;
4591 t0Shift = fT0Shift1;
4595 t0Shift = fT0Shift0;
4598 if ((fPhd[2] > fPhd[0]) &&
4599 (fPhd[2] > fPhd[1]) &&
4600 (fPhd[1] > fPhd[0]) &&
4602 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4603 fNumberFitSuccess++;
4604 if (fPhdt0 >= 0.0) {
4605 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4606 if (fCurrentCoef2[0] < -1.0) {
4607 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4611 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4615 if((fPhd[1] > fPhd[0]) &&
4617 if (fPhdt0 >= 0.0) {
4618 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4619 if (fCurrentCoef2[0] < -1.0) {
4620 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4624 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4628 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4629 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4630 //printf("Fit failed!\n");
4634 if (fDebugLevel == 1) {
4635 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4638 line->SetLineColor(2);
4639 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4640 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4641 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4642 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4643 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4644 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4645 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4646 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4649 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4654 if(pentea) delete pentea;
4655 if(pente) delete pente;
4656 if(polynome) delete polynome;
4657 if(polynomea) delete polynomea;
4658 if(polynomeb) delete polynomeb;
4662 if(line) delete line;
4667 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4668 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4670 projPH->SetDirectory(0);
4674 //_____________________________________________________________________________
4675 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4678 // Fit methode for the drift velocity
4682 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4685 TAxis *xpph = projPH->GetXaxis();
4686 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4688 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4689 fPH->SetParameter(0,0.469); // Scaling
4690 fPH->SetParameter(1,0.18); // Start
4691 fPH->SetParameter(2,0.0857325); // AR
4692 fPH->SetParameter(3,1.89); // DR
4693 fPH->SetParameter(4,0.08); // QA/QD
4694 fPH->SetParameter(5,0.0); // Baseline
4696 TLine *line = new TLine();
4698 fCurrentCoef[0] = 0.0;
4699 fCurrentCoef2[0] = 0.0;
4700 fCurrentCoefE = 0.0;
4701 fCurrentCoefE2 = 0.0;
4703 if (idect%fFitPHPeriode == 0) {
4705 AliInfo(Form("The detector %d will be fitted",idect));
4706 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4707 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4708 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4709 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4710 fPH->SetParameter(4,0.225); // QA/QD
4711 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4713 if (fDebugLevel != 1) {
4714 projPH->Fit(fPH,"0M","",0.0,upedge);
4717 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4719 projPH->Fit(fPH,"M+","",0.0,upedge);
4721 line->SetLineColor(4);
4722 line->DrawLine(fPH->GetParameter(1)
4724 ,fPH->GetParameter(1)
4725 ,projPH->GetMaximum());
4726 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4728 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4729 ,projPH->GetMaximum());
4730 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4732 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4733 ,projPH->GetMaximum());
4736 if (fPH->GetParameter(3) != 0) {
4737 fNumberFitSuccess++;
4738 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4739 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4740 fCurrentCoef2[0] = fPH->GetParameter(1);
4741 fCurrentCoefE2 = fPH->GetParError(1);
4744 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4745 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4751 // Put the default value
4752 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4753 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4756 if (fDebugLevel != 1) {
4761 //_____________________________________________________________________________
4762 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
4765 // Fit methode for the sigma of the pad response function
4770 fCurrentCoef[0] = 0.0;
4771 fCurrentCoefE = 0.0;
4773 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
4776 fCurrentCoef[0] = -fCurrentCoef[1];
4780 fNumberFitSuccess++;
4781 fCurrentCoef[0] = param[2];
4782 fCurrentCoefE = ret;
4786 //_____________________________________________________________________________
4787 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)
4790 // Fit methode for the sigma of the pad response function
4793 //We should have at least 3 points
4794 if(nBins <=3) return -4.0;
4796 TLinearFitter fitter(3,"pol2");
4797 fitter.StoreData(kFALSE);
4798 fitter.ClearPoints();
4800 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
4801 Float_t entries = 0;
4802 Int_t nbbinwithentries = 0;
4803 for (Int_t i=0; i<nBins; i++){
4805 if(arraye[i] > 15) nbbinwithentries++;
4806 //printf("entries for i %d: %f\n",i,arraye[i]);
4808 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
4809 //printf("entries %f\n",entries);
4810 //printf("nbbinwithentries %d\n",nbbinwithentries);
4813 Float_t errorm = 0.0;
4814 Float_t errorn = 0.0;
4815 Float_t error = 0.0;
4818 for (Int_t ibin=0;ibin<nBins; ibin++){
4819 Float_t entriesI = arraye[ibin];
4820 Float_t valueI = arraym[ibin];
4821 Double_t xcenter = 0.0;
4823 if ((entriesI>15) && (valueI>0.0)){
4824 xcenter = xMin+(ibin+0.5)*binWidth;
4829 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
4830 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
4831 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4834 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
4835 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
4836 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4838 if(error == 0.0) continue;
4839 val = TMath::Log(Float_t(valueI));
4840 fitter.AddPoint(&xcenter,val,error);
4844 if(fDebugLevel > 1){
4846 if ( !fDebugStreamer ) {
4848 TDirectory *backup = gDirectory;
4849 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4850 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4853 Int_t detector = fCountDet;
4854 Int_t layer = GetLayer(fCountDet);
4857 (* fDebugStreamer) << "FitGausMIFill"<<
4858 "detector="<<detector<<
4862 "entriesI="<<entriesI<<
4865 "xcenter="<<xcenter<<
4875 if(npoints <=3) return -4.0;
4880 fitter.GetParameters(par);
4881 chi2 = fitter.GetChisquare()/Float_t(npoints);
4884 if (!param) param = new TVectorD(3);
4885 if(par[2] == 0.0) return -4.0;
4886 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
4887 Double_t deltax = (fitter.GetParError(2))/x;
4888 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
4891 (*param)[1] = par[1]/(-2.*par[2]);
4892 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
4893 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
4894 if ( lnparam0>307 ) return -4;
4895 (*param)[0] = TMath::Exp(lnparam0);
4897 if(fDebugLevel > 1){
4899 if ( !fDebugStreamer ) {
4901 TDirectory *backup = gDirectory;
4902 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4903 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4906 Int_t detector = fCountDet;
4907 Int_t layer = GetLayer(fCountDet);
4910 (* fDebugStreamer) << "FitGausMIFit"<<
4911 "detector="<<detector<<
4914 "errorsigma="<<chi2<<
4915 "mean="<<(*param)[1]<<
4916 "sigma="<<(*param)[2]<<
4917 "constant="<<(*param)[0]<<
4922 if((chi2/(*param)[2]) > 0.1){
4924 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
4929 if(fDebugLevel == 1){
4930 TString name("PRF");
4931 name += (Int_t)xMin;
4932 name += (Int_t)xMax;
4933 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
4936 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
4937 for(Int_t k = 0; k < nBins; k++){
4938 histo->SetBinContent(k+1,arraym[k]);
4939 histo->SetBinError(k+1,arrayme[k]);
4942 name += "functionf";
4943 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
4944 f1->SetParameter(0, (*param)[0]);
4945 f1->SetParameter(1, (*param)[1]);
4946 f1->SetParameter(2, (*param)[2]);
4954 //_____________________________________________________________________________
4955 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
4958 // Fit methode for the sigma of the pad response function
4961 fCurrentCoef[0] = 0.0;
4962 fCurrentCoefE = 0.0;
4964 if (fDebugLevel != 1) {
4965 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
4968 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4970 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
4973 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
4974 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
4975 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4977 fNumberFitSuccess++;
4980 //_____________________________________________________________________________
4981 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
4984 // Fit methode for the sigma of the pad response function
4986 fCurrentCoef[0] = 0.0;
4987 fCurrentCoefE = 0.0;
4988 if (fDebugLevel == 1) {
4989 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4993 fCurrentCoef[0] = projPRF->GetRMS();
4994 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4996 fNumberFitSuccess++;
4999 //_____________________________________________________________________________
5000 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5003 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5006 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5009 Int_t nbins = (Int_t)(nybins/(2*nbg));
5010 Float_t lowedge = -3.0*nbg;
5011 Float_t upedge = lowedge + 3.0;
5014 Double_t xvalues = -0.2*nbg+0.1;
5016 Int_t total = 2*nbg;
5019 for(Int_t k = 0; k < total; k++){
5020 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5022 y = fCurrentCoef[0]*fCurrentCoef[0];
5023 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5026 if(fDebugLevel > 1){
5028 if ( !fDebugStreamer ) {
5030 TDirectory *backup = gDirectory;
5031 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5032 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5035 Int_t detector = fCountDet;
5036 Int_t layer = GetLayer(fCountDet);
5037 Int_t nbtotal = total;
5039 Float_t low = lowedge;
5040 Float_t up = upedge;
5041 Float_t tnp = xvalues;
5042 Float_t wid = fCurrentCoef[0];
5043 Float_t widfE = fCurrentCoefE;
5045 (* fDebugStreamer) << "FitTnpRange0"<<
5046 "detector="<<detector<<
5048 "nbtotal="<<nbtotal<<
5066 fCurrentCoefE = 0.0;
5067 fCurrentCoef[0] = 0.0;
5069 //printf("npoints\n",npoints);
5072 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5077 linearfitter.Eval();
5078 linearfitter.GetParameters(pars0);
5079 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5080 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5081 Double_t min0 = 0.0;
5082 Double_t ermin0 = 0.0;
5083 //Double_t prfe0 = 0.0;
5084 Double_t prf0 = 0.0;
5085 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
5086 min0 = -pars0[1]/(2*pars0[2]);
5087 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5088 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5091 prfe0 = linearfitter->GetParError(0)*pointError0
5092 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5093 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5094 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5095 fCurrentCoefE = (Float_t) prfe0;
5097 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5100 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5104 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5107 if(fDebugLevel > 1){
5109 if ( !fDebugStreamer ) {
5111 TDirectory *backup = gDirectory;
5112 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5113 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5116 Int_t detector = fCountDet;
5117 Int_t layer = GetLayer(fCountDet);
5118 Int_t nbtotal = total;
5119 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5120 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5122 (* fDebugStreamer) << "FitTnpRange1"<<
5123 "detector="<<detector<<
5125 "nbtotal="<<nbtotal<<
5129 "npoints="<<npoints<<
5132 "sigmaprf="<<fCurrentCoef[0]<<
5133 "sigprf="<<fCurrentCoef[1]<<
5140 //_____________________________________________________________________________
5141 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5144 // Only mean methode for the gain factor
5147 fCurrentCoef[0] = mean;
5148 fCurrentCoefE = 0.0;
5149 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5150 if (fDebugLevel == 1) {
5151 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5155 CalculChargeCoefMean(kTRUE);
5156 fNumberFitSuccess++;
5158 //_____________________________________________________________________________
5159 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5162 // mean w methode for the gain factor
5166 Int_t nybins = projch->GetNbinsX();
5168 //The weight function
5169 Double_t a = 0.00228515;
5170 Double_t b = -0.00231487;
5171 Double_t c = 0.00044298;
5172 Double_t d = -0.00379239;
5173 Double_t e = 0.00338349;
5183 //A arbitrary error for the moment
5184 fCurrentCoefE = 0.0;
5185 fCurrentCoef[0] = 0.0;
5188 Double_t sumw = 0.0;
5190 Float_t sumAll = (Float_t) nentries;
5191 Int_t sumCurrent = 0;
5192 for(Int_t k = 0; k <nybins; k++){
5193 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5194 if (fraction>0.95) break;
5195 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5196 e*fraction*fraction*fraction*fraction;
5197 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5198 sum += weight*projch->GetBinContent(k+1);
5199 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5200 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5202 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5204 if (fDebugLevel == 1) {
5205 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5209 fNumberFitSuccess++;
5210 CalculChargeCoefMean(kTRUE);
5212 //_____________________________________________________________________________
5213 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5216 // mean w methode for the gain factor
5220 Int_t nybins = projch->GetNbinsX();
5222 //The weight function
5223 Double_t a = 0.00228515;
5224 Double_t b = -0.00231487;
5225 Double_t c = 0.00044298;
5226 Double_t d = -0.00379239;
5227 Double_t e = 0.00338349;
5237 //A arbitrary error for the moment
5238 fCurrentCoefE = 0.0;
5239 fCurrentCoef[0] = 0.0;
5242 Double_t sumw = 0.0;
5244 Int_t sumCurrent = 0;
5245 for(Int_t k = 0; k <nybins; k++){
5246 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5247 if (fraction>0.95) break;
5248 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5249 e*fraction*fraction*fraction*fraction;
5250 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5251 sum += weight*projch->GetBinContent(k+1);
5252 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5253 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5255 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5257 if (fDebugLevel == 1) {
5258 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5262 fNumberFitSuccess++;
5264 //_____________________________________________________________________________
5265 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5268 // Fit methode for the gain factor
5271 fCurrentCoef[0] = 0.0;
5272 fCurrentCoefE = 0.0;
5273 Double_t chisqrl = 0.0;
5274 Double_t chisqrg = 0.0;
5275 Double_t chisqr = 0.0;
5276 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5278 projch->Fit("landau","0",""
5279 ,(Double_t) mean/fBeginFitCharge
5280 ,projch->GetBinCenter(projch->GetNbinsX()));
5281 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5282 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5283 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5284 chisqrl = projch->GetFunction("landau")->GetChisquare();
5286 projch->Fit("gaus","0",""
5287 ,(Double_t) mean/fBeginFitCharge
5288 ,projch->GetBinCenter(projch->GetNbinsX()));
5289 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5290 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5291 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5293 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5294 if (fDebugLevel != 1) {
5295 projch->Fit("fLandauGaus","0",""
5296 ,(Double_t) mean/fBeginFitCharge
5297 ,projch->GetBinCenter(projch->GetNbinsX()));
5298 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5301 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5303 projch->Fit("fLandauGaus","+",""
5304 ,(Double_t) mean/fBeginFitCharge
5305 ,projch->GetBinCenter(projch->GetNbinsX()));
5306 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5308 fLandauGaus->Draw("same");
5311 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5312 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5313 fNumberFitSuccess++;
5314 CalculChargeCoefMean(kTRUE);
5315 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5316 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5319 CalculChargeCoefMean(kFALSE);
5320 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5323 if (fDebugLevel != 1) {
5328 //_____________________________________________________________________________
5329 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5332 // Fit methode for the gain factor more time consuming
5336 //Some parameters to initialise
5337 Double_t widthLandau, widthGaus, mPV, integral;
5338 Double_t chisquarel = 0.0;
5339 Double_t chisquareg = 0.0;
5340 projch->Fit("landau","0M+",""
5342 ,projch->GetBinCenter(projch->GetNbinsX()));
5343 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5344 chisquarel = projch->GetFunction("landau")->GetChisquare();
5345 projch->Fit("gaus","0M+",""
5347 ,projch->GetBinCenter(projch->GetNbinsX()));
5348 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5349 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5351 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5352 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5354 // Setting fit range and start values
5356 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5357 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5358 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5359 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5360 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5361 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5362 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5365 fCurrentCoef[0] = 0.0;
5366 fCurrentCoefE = 0.0;
5370 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5375 Double_t projchPeak;
5376 Double_t projchFWHM;
5377 LanGauPro(fp,projchPeak,projchFWHM);
5379 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5380 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5381 fNumberFitSuccess++;
5382 CalculChargeCoefMean(kTRUE);
5383 fCurrentCoef[0] = fp[1];
5384 fCurrentCoefE = fpe[1];
5385 //chargeCoefE2 = chisqr;
5388 CalculChargeCoefMean(kFALSE);
5389 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5391 if (fDebugLevel == 1) {
5392 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5393 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5396 fitsnr->Draw("same");
5402 //_____________________________________________________________________________
5403 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
5406 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5408 Double_t *c = new Double_t[5];
5409 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5410 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5411 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5416 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5417 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5424 //_____________________________________________________________________________
5425 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
5428 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5430 Double_t *c = new Double_t[5];
5431 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5432 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5433 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5434 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5438 c[2] = -(x0*(x[1]+x[2]+x[3])
5439 +x1*(x[0]+x[2]+x[3])
5440 +x2*(x[0]+x[1]+x[3])
5441 +x3*(x[0]+x[1]+x[2]));
5442 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5443 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5444 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5445 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5447 c[0] = -(x0*x[1]*x[2]*x[3]
5450 +x3*x[0]*x[1]*x[2]);
5458 //_____________________________________________________________________________
5459 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
5462 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5464 Double_t *c = new Double_t[5];
5465 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5466 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5467 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5468 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5469 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5472 c[4] = x0+x1+x2+x3+x4;
5473 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
5474 +x1*(x[0]+x[2]+x[3]+x[4])
5475 +x2*(x[0]+x[1]+x[3]+x[4])
5476 +x3*(x[0]+x[1]+x[2]+x[4])
5477 +x4*(x[0]+x[1]+x[2]+x[3]));
5478 c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5479 +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])
5480 +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])
5481 +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])
5482 +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]));
5484 c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5485 +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])
5486 +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])
5487 +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])
5488 +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]));
5490 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
5491 +x1*x[0]*x[2]*x[3]*x[4]
5492 +x2*x[0]*x[1]*x[3]*x[4]
5493 +x3*x[0]*x[1]*x[2]*x[4]
5494 +x4*x[0]*x[1]*x[2]*x[3]);
5500 //_____________________________________________________________________________
5501 void AliTRDCalibraFit::NormierungCharge()
5504 // Normalisation of the gain factor resulting for the fits
5507 // Calcul of the mean of choosen method by fFitChargeNDB
5509 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5510 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5512 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5513 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5514 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5515 if (GetStack(detector) == 2) {
5518 if (GetStack(detector) != 2) {
5521 for (Int_t j = 0; j < total; j++) {
5529 fScaleFitFactor = fScaleFitFactor / sum;
5532 fScaleFitFactor = 1.0;
5535 //methode de boeuf mais bon...
5536 Double_t scalefactor = fScaleFitFactor;
5538 if(fDebugLevel > 1){
5540 if ( !fDebugStreamer ) {
5542 TDirectory *backup = gDirectory;
5543 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5544 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5546 (* fDebugStreamer) << "NormierungCharge"<<
5547 "scalefactor="<<scalefactor<<
5551 //_____________________________________________________________________________
5552 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5555 // Rebin of the 1D histo for the gain calibration if needed.
5556 // you have to choose fRebin, divider of fNumberBinCharge
5559 TAxis *xhist = hist->GetXaxis();
5560 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5561 ,xhist->GetBinLowEdge(1)
5562 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5564 AliInfo(Form("fRebin: %d",fRebin));
5566 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5568 for (Int_t ji = i; ji < i+fRebin; ji++) {
5569 sum += hist->GetBinContent(ji);
5572 rehist->SetBinContent(k,sum);
5580 //_____________________________________________________________________________
5581 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5584 // Rebin of the 1D histo for the gain calibration if needed
5585 // you have to choose fRebin divider of fNumberBinCharge
5588 TAxis *xhist = hist->GetXaxis();
5589 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5590 ,xhist->GetBinLowEdge(1)
5591 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5593 AliInfo(Form("fRebin: %d",fRebin));
5595 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5597 for (Int_t ji = i; ji < i+fRebin; ji++) {
5598 sum += hist->GetBinContent(ji);
5601 rehist->SetBinContent(k,sum);
5609 //____________Some basic geometry function_____________________________________
5612 //_____________________________________________________________________________
5613 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5616 // Reconstruct the plane number from the detector number
5619 return ((Int_t) (d % 6));
5623 //_____________________________________________________________________________
5624 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5627 // Reconstruct the stack number from the detector number
5629 const Int_t kNlayer = 6;
5631 return ((Int_t) (d % 30) / kNlayer);
5635 //_____________________________________________________________________________
5636 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5639 // Reconstruct the sector number from the detector number
5643 return ((Int_t) (d / fg));
5648 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5650 //_______________________________________________________________________________
5651 void AliTRDCalibraFit::ResetVectorFit()
5654 // Reset the VectorFits
5657 fVectorFit.SetOwner();
5659 fVectorFit2.SetOwner();
5660 fVectorFit2.Clear();
5664 //____________Private Functions________________________________________________
5667 //_____________________________________________________________________________
5668 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5671 // Function for the fit
5674 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5676 //PARAMETERS FOR FIT PH
5678 //fAsymmGauss->SetParameter(0,0.113755);
5679 //fAsymmGauss->SetParameter(1,0.350706);
5680 //fAsymmGauss->SetParameter(2,0.0604244);
5681 //fAsymmGauss->SetParameter(3,7.65596);
5682 //fAsymmGauss->SetParameter(4,1.00124);
5683 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5691 Double_t dx = 0.005;
5692 Double_t xs = par[1];
5694 Double_t paras[2] = { 0.0, 0.0 };
5697 if ((xs >= par[1]) &&
5698 (xs < (par[1]+par[2]))) {
5699 //fAsymmGauss->SetParameter(0,par[0]);
5700 //fAsymmGauss->SetParameter(1,xs);
5701 //ss += fAsymmGauss->Eval(xx);
5704 ss += AsymmGauss(&xx,paras);
5706 if ((xs >= (par[1]+par[2])) &&
5707 (xs < (par[1]+par[2]+par[3]))) {
5708 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5709 //fAsymmGauss->SetParameter(1,xs);
5710 //ss += fAsymmGauss->Eval(xx);
5711 paras[0] = par[0]*par[4];
5713 ss += AsymmGauss(&xx,paras);
5722 //_____________________________________________________________________________
5723 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5726 // Function for the fit
5729 //par[0] = normalization
5737 Double_t par1save = par[1];
5738 //Double_t par2save = par[2];
5739 Double_t par2save = 0.0604244;
5740 //Double_t par3save = par[3];
5741 Double_t par3save = 7.65596;
5742 //Double_t par5save = par[5];
5743 Double_t par5save = 0.870597;
5744 Double_t dx = x[0] - par1save;
5746 Double_t sigma2 = par2save*par2save;
5747 Double_t sqrt2 = TMath::Sqrt(2.0);
5748 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5749 * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5750 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5751 * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5753 //return par[0]*(exp1+par[4]*exp2);
5754 return par[0] * (exp1 + 1.00124 * exp2);
5758 //_____________________________________________________________________________
5759 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5762 // Sum Landau + Gaus with identical mean
5765 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5766 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5767 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5768 Double_t val = valLandau + valGaus;
5774 //_____________________________________________________________________________
5775 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
5778 // Function for the fit
5781 // par[0]=Width (scale) parameter of Landau density
5782 // par[1]=Most Probable (MP, location) parameter of Landau density
5783 // par[2]=Total area (integral -inf to inf, normalization constant)
5784 // par[3]=Width (sigma) of convoluted Gaussian function
5786 // In the Landau distribution (represented by the CERNLIB approximation),
5787 // the maximum is located at x=-0.22278298 with the location parameter=0.
5788 // This shift is corrected within this function, so that the actual
5789 // maximum is identical to the MP parameter.
5792 // Numeric constants
5793 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
5794 Double_t mpshift = -0.22278298; // Landau maximum location
5796 // Control constants
5797 Double_t np = 100.0; // Number of convolution steps
5798 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
5810 // MP shift correction
5811 mpc = par[1] - mpshift * par[0];
5813 // Range of convolution integral
5814 xlow = x[0] - sc * par[3];
5815 xupp = x[0] + sc * par[3];
5817 step = (xupp - xlow) / np;
5819 // Convolution integral of Landau and Gaussian by sum
5820 for (i = 1.0; i <= np/2; i++) {
5822 xx = xlow + (i-.5) * step;
5823 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5824 sum += fland * TMath::Gaus(x[0],xx,par[3]);
5826 xx = xupp - (i-.5) * step;
5827 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5828 sum += fland * TMath::Gaus(x[0],xx,par[3]);
5832 return (par[2] * step * sum * invsq2pi / par[3]);
5835 //_____________________________________________________________________________
5836 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
5837 , const Double_t *parlimitslo, const Double_t *parlimitshi
5838 , Double_t *fitparams, Double_t *fiterrors
5839 , Double_t *chiSqr, Int_t *ndf) const
5842 // Function for the fit
5846 Char_t funname[100];
5848 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
5853 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
5854 ffit->SetParameters(startvalues);
5855 ffit->SetParNames("Width","MP","Area","GSigma");
5857 for (i = 0; i < 4; i++) {
5858 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
5861 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
5863 ffit->GetParameters(fitparams); // Obtain fit parameters
5864 for (i = 0; i < 4; i++) {
5865 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
5867 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
5868 ndf[0] = ffit->GetNDF(); // Obtain ndf
5870 return (ffit); // Return fit function
5874 //_____________________________________________________________________________
5875 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
5878 // Function for the fit
5891 Int_t maxcalls = 10000;
5893 // Search for maximum
5894 p = params[1] - 0.1 * params[0];
5895 step = 0.05 * params[0];
5899 while ((l != lold) && (i < maxcalls)) {
5903 l = LanGauFun(&x,params);
5905 step = -step / 10.0;
5910 if (i == maxcalls) {
5916 // Search for right x location of fy
5917 p = maxx + params[0];
5923 while ( (l != lold) && (i < maxcalls) ) {
5928 l = TMath::Abs(LanGauFun(&x,params) - fy);
5942 // Search for left x location of fy
5944 p = maxx - 0.5 * params[0];
5950 while ((l != lold) && (i < maxcalls)) {
5954 l = TMath::Abs(LanGauFun(&x,params) - fy);
5956 step = -step / 10.0;
5961 if (i == maxcalls) {
5970 //_____________________________________________________________________________
5971 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
5974 // Gaus with identical mean
5977 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];