1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
62 #include <TTreeStream.h>
63 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
71 #include "AliTRDCalibraFit.h"
72 #include "AliTRDCalibraMode.h"
73 #include "AliTRDCalibraVector.h"
74 #include "AliTRDCalibraVdriftLinearFit.h"
75 #include "AliTRDCalibraExbAltFit.h"
76 #include "AliTRDcalibDB.h"
77 #include "AliTRDgeometry.h"
78 #include "AliTRDpadPlane.h"
79 #include "AliTRDgeometry.h"
80 #include "AliTRDCommonParam.h"
81 #include "./Cal/AliTRDCalROC.h"
82 #include "./Cal/AliTRDCalPad.h"
83 #include "./Cal/AliTRDCalDet.h"
86 ClassImp(AliTRDCalibraFit)
88 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
89 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
91 //_____________singleton implementation_________________________________________________
92 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
95 // Singleton implementation
98 if (fgTerminated != kFALSE) {
102 if (fgInstance == 0) {
103 fgInstance = new AliTRDCalibraFit();
109 //______________________________________________________________________________________
110 void AliTRDCalibraFit::Terminate()
113 // Singleton implementation
114 // Deletes the instance of this class
117 fgTerminated = kTRUE;
119 if (fgInstance != 0) {
125 //______________________________________________________________________________________
126 AliTRDCalibraFit::AliTRDCalibraFit()
129 ,fNumberOfBinsExpected(0)
131 ,fBeginFitCharge(3.5)
133 ,fTakeTheMaxPH(kTRUE)
142 ,fNumberFitSuccess(0)
149 ,fCalibraMode(new AliTRDCalibraMode())
154 ,fScaleFitFactor(0.0)
163 ,fCalDetVdriftUsed(0x0)
165 ,fCurrentCoefDetector(0x0)
166 ,fCurrentCoefDetector2(0x0)
171 // Default constructor
174 fGeo = new AliTRDgeometry();
176 // Current variables initialised
177 for (Int_t k = 0; k < 2; k++) {
178 fCurrentCoef[k] = 0.0;
179 fCurrentCoef2[k] = 0.0;
181 for (Int_t i = 0; i < 3; i++) {
187 //______________________________________________________________________________________
188 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
191 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
193 ,fBeginFitCharge(c.fBeginFitCharge)
194 ,fFitPHPeriode(c.fFitPHPeriode)
195 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
196 ,fT0Shift0(c.fT0Shift0)
197 ,fT0Shift1(c.fT0Shift1)
198 ,fRangeFitPRF(c.fRangeFitPRF)
200 ,fMinEntries(c.fMinEntries)
202 ,fScaleGain(c.fScaleGain)
203 ,fNumberFit(c.fNumberFit)
204 ,fNumberFitSuccess(c.fNumberFitSuccess)
205 ,fNumberEnt(c.fNumberEnt)
206 ,fStatisticMean(c.fStatisticMean)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fFitVoir(c.fFitVoir)
210 ,fMagneticField(c.fMagneticField)
212 ,fCurrentCoefE(c.fCurrentCoefE)
213 ,fCurrentCoefE2(c.fCurrentCoefE2)
216 ,fScaleFitFactor(c.fScaleFitFactor)
217 ,fEntriesCurrent(c.fEntriesCurrent)
218 ,fCountDet(c.fCountDet)
225 ,fCalDetVdriftUsed(0x0)
227 ,fCurrentCoefDetector(0x0)
228 ,fCurrentCoefDetector2(0x0)
236 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
238 //Current variables initialised
239 for (Int_t k = 0; k < 2; k++) {
240 fCurrentCoef[k] = 0.0;
241 fCurrentCoef2[k] = 0.0;
243 for (Int_t i = 0; i < 3; i++) {
247 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
248 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
250 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
251 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
253 if(c.fCalDetVdriftUsed) fCalDetVdriftUsed = new AliTRDCalDet(*c.fCalDetVdriftUsed);
254 if(c.fCalDetExBUsed) fCalDetExBUsed = new AliTRDCalDet(*c.fCalDetExBUsed);
256 fVectorFit.SetName(c.fVectorFit.GetName());
257 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
258 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
259 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
261 if (GetStack(detector) == 2) {
267 Float_t *coef = new Float_t[ntotal];
268 for (Int_t i = 0; i < ntotal; i++) {
269 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
271 fitInfo->SetCoef(coef);
272 fitInfo->SetDetector(detector);
273 fVectorFit.Add((TObject *) fitInfo);
275 fVectorFit.SetName(c.fVectorFit.GetName());
276 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
277 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
278 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
280 if (GetStack(detector) == 2) {
286 Float_t *coef = new Float_t[ntotal];
287 for (Int_t i = 0; i < ntotal; i++) {
288 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
290 fitInfo->SetCoef(coef);
291 fitInfo->SetDetector(detector);
292 fVectorFit2.Add((TObject *) fitInfo);
297 fGeo = new AliTRDgeometry();
300 //____________________________________________________________________________________
301 AliTRDCalibraFit::~AliTRDCalibraFit()
304 // AliTRDCalibraFit destructor
306 if ( fDebugStreamer ) delete fDebugStreamer;
307 if ( fCalDet ) delete fCalDet;
308 if ( fCalDet2 ) delete fCalDet2;
309 if ( fCalROC ) delete fCalROC;
310 if ( fCalROC2 ) delete fCalROC2;
311 if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed;
312 if ( fCalDetExBUsed) delete fCalDetExBUsed;
313 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
314 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
316 fVectorFit2.Delete();
322 //_____________________________________________________________________________
323 void AliTRDCalibraFit::Destroy()
335 //_____________________________________________________________________________
336 void AliTRDCalibraFit::DestroyDebugStreamer()
339 // Delete DebugStreamer
342 if ( fDebugStreamer ) delete fDebugStreamer;
343 fDebugStreamer = 0x0;
346 //__________________________________________________________________________________
347 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
350 // From the drift velocity and t0
351 // return the position of the peak and maximum negative slope
354 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
355 Double_t widbins = 0.1; // 0.1 mus
357 //peak and maxnegslope in mus
358 Double_t begind = t0*widbins + fT0Shift0;
359 Double_t peakd = t0*widbins + fT0Shift1;
360 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
362 // peak and maxnegslope in timebin
363 begin = TMath::Nint(begind*widbins);
364 peak = TMath::Nint(peakd*widbins);
365 end = TMath::Nint(maxnegslope*widbins);
368 //____________Functions fit Online CH2d________________________________________
369 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
372 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
373 // calibration group normalized the resulted coefficients (to 1 normally)
376 // Set the calibration mode
377 //const char *name = ch->GetTitle();
378 TString name = ch->GetTitle();
379 if(!SetModeCalibration(name,0)) return kFALSE;
381 // Number of Ybins (detectors or groups of pads)
382 Int_t nbins = ch->GetNbinsX();// charge
383 Int_t nybins = ch->GetNbinsY();// groups number
384 if (!InitFit(nybins,0)) {
390 fStatisticMean = 0.0;
392 fNumberFitSuccess = 0;
394 // Init fCountDet and fCount
395 InitfCountDetAndfCount(0);
396 // Beginning of the loop betwwen dect1 and dect2
397 for (Int_t idect = fDect1; idect < fDect2; idect++) {
398 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
399 UpdatefCountDetAndfCount(idect,0);
400 ReconstructFitRowMinRowMax(idect, 0);
402 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
403 projch->SetDirectory(0);
404 // Number of entries for this calibration group
405 Double_t nentries = 0.0;
407 for (Int_t k = 0; k < nbins; k++) {
408 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
409 nentries += ch->GetBinContent(binnb);
410 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
411 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
413 projch->SetEntries(nentries);
414 //printf("The number of entries for the group %d is %f\n",idect,nentries);
419 // Rebin and statistic stuff
421 projch = ReBin((TH1I *) projch);
423 // This detector has not enough statistics or was off
424 if (nentries <= fMinEntries) {
425 NotEnoughStatisticCH(idect);
426 if (fDebugLevel != 1) {
431 // Statistics of the group fitted
432 fStatisticMean += nentries;
437 case 0: FitMeanW((TH1 *) projch, nentries); break;
438 case 1: FitMean((TH1 *) projch, nentries, mean); break;
439 case 2: FitCH((TH1 *) projch, mean); break;
440 case 3: FitBisCH((TH1 *) projch, mean); break;
441 default: return kFALSE;
444 FillInfosFitCH(idect);
446 if (fDebugLevel != 1) {
451 if (fDebugLevel != 1) {
455 if (fNumberFit > 0) {
456 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));
457 fStatisticMean = fStatisticMean / fNumberFit;
460 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
462 delete fDebugStreamer;
463 fDebugStreamer = 0x0;
467 //____________Functions fit Online CH2d________________________________________
468 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
471 // Reconstruct a 1D histo from the vectorCH for each calibration group,
472 // fit the histo, normalized the resulted coefficients (to 1 normally)
475 // Set the calibraMode
476 //const char *name = calvect->GetNameCH();
477 TString name = calvect->GetNameCH();
478 if(!SetModeCalibration(name,0)) return kFALSE;
480 // Number of Xbins (detectors or groups of pads)
481 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
487 fStatisticMean = 0.0;
489 fNumberFitSuccess = 0;
491 // Init fCountDet and fCount
492 InitfCountDetAndfCount(0);
493 // Beginning of the loop between dect1 and dect2
494 for (Int_t idect = fDect1; idect < fDect2; idect++) {
495 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
496 UpdatefCountDetAndfCount(idect,0);
497 ReconstructFitRowMinRowMax(idect,0);
499 Double_t nentries = 0.0;
501 if(!calvect->GetCHEntries(fCountDet)) {
502 NotEnoughStatisticCH(idect);
508 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
509 projch->SetDirectory(0);
510 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
511 nentries += projch->GetBinContent(k+1);
512 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
518 //printf("The number of entries for the group %d is %f\n",idect,nentries);
521 projch = ReBin((TH1F *) projch);
523 // This detector has not enough statistics or was not found in VectorCH
524 if (nentries <= fMinEntries) {
525 NotEnoughStatisticCH(idect);
528 // Statistic of the histos fitted
529 fStatisticMean += nentries;
534 case 0: FitMeanW((TH1 *) projch, nentries); break;
535 case 1: FitMean((TH1 *) projch, nentries, mean); break;
536 case 2: FitCH((TH1 *) projch, mean); break;
537 case 3: FitBisCH((TH1 *) projch, mean); break;
538 default: return kFALSE;
541 FillInfosFitCH(idect);
544 if (fDebugLevel != 1) {
548 if (fNumberFit > 0) {
549 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));
550 fStatisticMean = fStatisticMean / fNumberFit;
553 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
555 delete fDebugStreamer;
556 fDebugStreamer = 0x0;
559 //____________Functions fit Online CH2d________________________________________
560 Double_t AliTRDCalibraFit::AnalyseCHAllTogether(const TH2I *ch)
563 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
564 // calibration group normalized the resulted coefficients (to 1 normally)
566 Int_t nbins = ch->GetNbinsX();// charge
567 Int_t nybins = ch->GetNbinsY();// groups number
569 TH1I *projch = (TH1I *) ch->ProjectionX("projch",1,nybins+1,(Option_t *)"e");
570 projch->SetDirectory(0);
571 // Number of entries for this calibration group
572 Double_t nentries = 0.0;
574 for (Int_t k = 0; k < nbins; k++) {
575 nentries += projch->GetBinContent(k+1);
576 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
577 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
579 projch->SetEntries(nentries);
580 //printf("The number of entries for the group %d is %f\n",idect,nentries);
584 // This detector has not enough statistics or was off
585 if (nentries <= fMinEntries) {
587 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
593 case 0: FitMeanW((TH1 *) projch, nentries); break;
594 case 1: FitMean((TH1 *) projch, nentries, mean); break;
595 case 2: FitCH((TH1 *) projch, mean); break;
596 case 3: FitBisCH((TH1 *) projch, mean); break;
597 default: return -100.0;
599 delete fDebugStreamer;
600 fDebugStreamer = 0x0;
602 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
606 //________________functions fit Online PH2d____________________________________
607 Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph)
610 // Take the 1D profiles (average pulse height), projections of the 2D PH
611 // on the Xaxis, for each calibration group
612 // Reconstruct a drift velocity
613 // A first calibration of T0 is also made using the same method
616 // Number of Xbins (detectors or groups of pads)
617 Int_t nbins = ph->GetNbinsX();// time
618 Int_t nybins = ph->GetNbinsY();// calibration group
621 TH1D *projph = (TH1D *) ph->ProjectionX("projph",1,nybins+1,(Option_t *) "e");
622 projph->SetDirectory(0);
623 // Number of entries for this calibration group
624 Double_t nentries = 0;
625 for(Int_t idect = 0; idect < nybins; idect++){
626 for (Int_t k = 0; k < nbins; k++) {
627 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
628 nentries += ph->GetBinEntries(binnb);
631 //printf("AnalysePHAllTogether:: the number of entries is %f\n",nentries);
632 // This detector has not enough statistics or was off
633 if (nentries <= fMinEntries) {
634 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
635 if (fDebugLevel != 1) {
641 //printf("Method\n");
644 case 0: FitLagrangePoly((TH1 *) projph); break;
645 case 1: FitPente((TH1 *) projph); break;
646 case 2: FitPH((TH1 *) projph,0); break;
647 default: return -100.0;
650 if (fDebugLevel != 1) {
653 delete fDebugStreamer;
654 fDebugStreamer = 0x0;
656 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
660 //____________Functions fit Online PH2d________________________________________
661 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
664 // Reconstruct the average pulse height from the vectorPH for each
666 // Reconstruct a drift velocity
667 // A first calibration of T0 is also made using the same method (slope method)
670 // Set the calibration mode
671 //const char *name = calvect->GetNamePH();
672 TString name = calvect->GetNamePH();
673 if(!SetModeCalibration(name,1)) return kFALSE;
675 // Number of Xbins (detectors or groups of pads)
676 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
682 fStatisticMean = 0.0;
684 fNumberFitSuccess = 0;
686 // Init fCountDet and fCount
687 InitfCountDetAndfCount(1);
688 // Beginning of the loop
689 for (Int_t idect = fDect1; idect < fDect2; idect++) {
690 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
691 UpdatefCountDetAndfCount(idect,1);
692 ReconstructFitRowMinRowMax(idect,1);
695 if(!calvect->GetPHEntries(fCountDet)) {
696 NotEnoughStatisticPH(idect,fEntriesCurrent);
701 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
702 projph->SetDirectory(0);
703 if(fEntriesCurrent > 0) fNumberEnt++;
704 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
705 // This detector has not enough statistics or was off
706 if (fEntriesCurrent <= fMinEntries) {
707 //printf("Not enough stat!\n");
708 NotEnoughStatisticPH(idect,fEntriesCurrent);
711 // Statistic of the histos fitted
713 fStatisticMean += fEntriesCurrent;
714 // Calcul of "real" coef
715 CalculVdriftCoefMean();
720 case 0: FitLagrangePoly((TH1 *) projph); break;
721 case 1: FitPente((TH1 *) projph); break;
722 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
723 default: return kFALSE;
725 // Fill the tree if end of a detector or only the pointer to the branch!!!
726 FillInfosFitPH(idect,fEntriesCurrent);
730 if (fNumberFit > 0) {
731 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));
732 fStatisticMean = fStatisticMean / fNumberFit;
735 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
737 delete fDebugStreamer;
738 fDebugStreamer = 0x0;
741 //________________functions fit Online PH2d____________________________________
742 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
745 // Take the 1D profiles (average pulse height), projections of the 2D PH
746 // on the Xaxis, for each calibration group
747 // Reconstruct a drift velocity
748 // A first calibration of T0 is also made using the same method
751 // Set the calibration mode
752 //const char *name = ph->GetTitle();
753 TString name = ph->GetTitle();
754 if(!SetModeCalibration(name,1)) return kFALSE;
756 //printf("Mode calibration set\n");
758 // Number of Xbins (detectors or groups of pads)
759 Int_t nbins = ph->GetNbinsX();// time
760 Int_t nybins = ph->GetNbinsY();// calibration group
761 if (!InitFit(nybins,1)) {
765 //printf("Init fit\n");
771 //printf("Init fit PH\n");
773 fStatisticMean = 0.0;
775 fNumberFitSuccess = 0;
777 // Init fCountDet and fCount
778 InitfCountDetAndfCount(1);
779 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
781 // Beginning of the loop
782 for (Int_t idect = fDect1; idect < fDect2; idect++) {
783 //printf("idect = %d\n",idect);
784 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
785 UpdatefCountDetAndfCount(idect,1);
786 ReconstructFitRowMinRowMax(idect,1);
788 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
789 projph->SetDirectory(0);
790 // Number of entries for this calibration group
791 Double_t nentries = 0;
792 for (Int_t k = 0; k < nbins; k++) {
793 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
794 nentries += ph->GetBinEntries(binnb);
799 //printf("The number of entries for the group %d is %f\n",idect,nentries);
800 // This detector has not enough statistics or was off
801 if (nentries <= fMinEntries) {
802 //printf("Not enough statistic!\n");
803 NotEnoughStatisticPH(idect,nentries);
804 if (fDebugLevel != 1) {
809 // Statistics of the histos fitted
811 fStatisticMean += nentries;
812 // Calcul of "real" coef
813 CalculVdriftCoefMean();
816 //printf("Method\n");
819 case 0: FitLagrangePoly((TH1 *) projph); break;
820 case 1: FitPente((TH1 *) projph); break;
821 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
822 default: return kFALSE;
824 // Fill the tree if end of a detector or only the pointer to the branch!!!
825 FillInfosFitPH(idect,nentries);
827 if (fDebugLevel != 1) {
832 if (fNumberFit > 0) {
833 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));
834 fStatisticMean = fStatisticMean / fNumberFit;
837 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
839 delete fDebugStreamer;
840 fDebugStreamer = 0x0;
843 //____________Functions fit Online PRF2d_______________________________________
844 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
847 // Take the 1D profiles (pad response function), projections of the 2D PRF
848 // on the Xaxis, for each calibration group
849 // Fit with a gaussian to reconstruct the sigma of the pad response function
852 // Set the calibration mode
853 //const char *name = prf->GetTitle();
854 TString name = prf->GetTitle();
855 if(!SetModeCalibration(name,2)) return kFALSE;
857 // Number of Ybins (detectors or groups of pads)
858 Int_t nybins = prf->GetNbinsY();// calibration groups
859 Int_t nbins = prf->GetNbinsX();// bins
860 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
861 if((nbg > 0) || (nbg == -1)) return kFALSE;
862 if (!InitFit(nybins,2)) {
868 fStatisticMean = 0.0;
870 fNumberFitSuccess = 0;
872 // Init fCountDet and fCount
873 InitfCountDetAndfCount(2);
874 // Beginning of the loop
875 for (Int_t idect = fDect1; idect < fDect2; idect++) {
876 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
877 UpdatefCountDetAndfCount(idect,2);
878 ReconstructFitRowMinRowMax(idect,2);
880 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
881 projprf->SetDirectory(0);
882 // Number of entries for this calibration group
883 Double_t nentries = 0;
884 for (Int_t k = 0; k < nbins; k++) {
885 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
886 nentries += prf->GetBinEntries(binnb);
888 if(nentries > 0) fNumberEnt++;
889 // This detector has not enough statistics or was off
890 if (nentries <= fMinEntries) {
891 NotEnoughStatisticPRF(idect);
892 if (fDebugLevel != 1) {
897 // Statistics of the histos fitted
899 fStatisticMean += nentries;
900 // Calcul of "real" coef
905 case 0: FitPRF((TH1 *) projprf); break;
906 case 1: RmsPRF((TH1 *) projprf); break;
907 default: return kFALSE;
909 // Fill the tree if end of a detector or only the pointer to the branch!!!
910 FillInfosFitPRF(idect);
912 if (fDebugLevel != 1) {
917 if (fNumberFit > 0) {
918 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
919 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
920 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
921 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
922 fStatisticMean = fStatisticMean / fNumberFit;
925 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
927 delete fDebugStreamer;
928 fDebugStreamer = 0x0;
931 //____________Functions fit Online PRF2d_______________________________________
932 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
935 // Take the 1D profiles (pad response function), projections of the 2D PRF
936 // on the Xaxis, for each calibration group
937 // Fit with a gaussian to reconstruct the sigma of the pad response function
940 // Set the calibration mode
941 //const char *name = prf->GetTitle();
942 TString name = prf->GetTitle();
943 if(!SetModeCalibration(name,2)) return kFALSE;
945 // Number of Ybins (detectors or groups of pads)
946 TAxis *xprf = prf->GetXaxis();
947 TAxis *yprf = prf->GetYaxis();
948 Int_t nybins = yprf->GetNbins();// calibration groups
949 Int_t nbins = xprf->GetNbins();// bins
950 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
951 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
952 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
953 if(nbg == -1) return kFALSE;
954 if(nbg > 0) fMethod = 1;
956 if (!InitFit(nybins,2)) {
962 fStatisticMean = 0.0;
964 fNumberFitSuccess = 0;
966 // Init fCountDet and fCount
967 InitfCountDetAndfCount(2);
968 // Beginning of the loop
969 for (Int_t idect = fDect1; idect < fDect2; idect++) {
970 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
971 UpdatefCountDetAndfCount(idect,2);
972 ReconstructFitRowMinRowMax(idect,2);
973 // Build the array of entries and sum
974 TArrayD arraye = TArrayD(nbins);
975 TArrayD arraym = TArrayD(nbins);
976 TArrayD arrayme = TArrayD(nbins);
977 Double_t nentries = 0;
978 //printf("nbins %d\n",nbins);
979 for (Int_t k = 0; k < nbins; k++) {
980 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
981 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
982 Double_t mean = (Double_t)prf->GetBinContent(binnb);
983 Double_t error = (Double_t)prf->GetBinError(binnb);
984 //printf("for %d we have %f\n",k,entries);
986 arraye.AddAt(entries,k);
987 arraym.AddAt(mean,k);
988 arrayme.AddAt(error,k);
990 if(nentries > 0) fNumberEnt++;
991 //printf("The number of entries for the group %d is %f\n",idect,nentries);
992 // This detector has not enough statistics or was off
993 if (nentries <= fMinEntries) {
994 NotEnoughStatisticPRF(idect);
997 // Statistics of the histos fitted
999 fStatisticMean += nentries;
1000 // Calcul of "real" coef
1001 CalculPRFCoefMean();
1005 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
1006 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
1007 default: return kFALSE;
1009 // Fill the tree if end of a detector or only the pointer to the branch!!!
1010 FillInfosFitPRF(idect);
1013 if (fNumberFit > 0) {
1014 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1015 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1016 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
1017 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1018 fStatisticMean = fStatisticMean / fNumberFit;
1021 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1023 delete fDebugStreamer;
1024 fDebugStreamer = 0x0;
1027 //____________Functions fit Online PRF2d_______________________________________
1028 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
1031 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1032 // each calibration group
1033 // Fit with a gaussian to reconstruct the sigma of the pad response function
1036 // Set the calibra mode
1037 //const char *name = calvect->GetNamePRF();
1038 TString name = calvect->GetNamePRF();
1039 if(!SetModeCalibration(name,2)) return kFALSE;
1040 //printf("test0 %s\n",name);
1042 // Number of Xbins (detectors or groups of pads)
1043 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1044 //printf("test1\n");
1047 if (!InitFitPRF()) {
1048 ///printf("test2\n");
1051 fStatisticMean = 0.0;
1053 fNumberFitSuccess = 0;
1055 // Init fCountDet and fCount
1056 InitfCountDetAndfCount(2);
1057 // Beginning of the loop
1058 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1059 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
1060 UpdatefCountDetAndfCount(idect,2);
1061 ReconstructFitRowMinRowMax(idect,2);
1063 fEntriesCurrent = 0;
1064 if(!calvect->GetPRFEntries(fCountDet)) {
1065 NotEnoughStatisticPRF(idect);
1068 TString tname("PRF");
1070 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
1071 projprf->SetDirectory(0);
1072 if(fEntriesCurrent > 0) fNumberEnt++;
1073 // This detector has not enough statistics or was off
1074 if (fEntriesCurrent <= fMinEntries) {
1075 NotEnoughStatisticPRF(idect);
1078 // Statistic of the histos fitted
1080 fStatisticMean += fEntriesCurrent;
1081 // Calcul of "real" coef
1082 CalculPRFCoefMean();
1086 case 1: FitPRF((TH1 *) projprf); break;
1087 case 2: RmsPRF((TH1 *) projprf); break;
1088 default: return kFALSE;
1090 // Fill the tree if end of a detector or only the pointer to the branch!!!
1091 FillInfosFitPRF(idect);
1094 if (fNumberFit > 0) {
1095 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));
1098 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1100 delete fDebugStreamer;
1101 fDebugStreamer = 0x0;
1104 //____________Functions fit Online PRF2d_______________________________________
1105 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
1108 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1109 // each calibration group
1110 // Fit with a gaussian to reconstruct the sigma of the pad response function
1113 // Set the calibra mode
1114 //const char *name = calvect->GetNamePRF();
1115 TString name = calvect->GetNamePRF();
1116 if(!SetModeCalibration(name,2)) return kFALSE;
1117 //printf("test0 %s\n",name);
1118 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1119 //printf("test1 %d\n",nbg);
1120 if(nbg == -1) return kFALSE;
1121 if(nbg > 0) fMethod = 1;
1123 // Number of Xbins (detectors or groups of pads)
1124 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1125 //printf("test2\n");
1128 if (!InitFitPRF()) {
1129 //printf("test3\n");
1132 fStatisticMean = 0.0;
1134 fNumberFitSuccess = 0;
1138 Double_t *arrayx = 0;
1139 Double_t *arraye = 0;
1140 Double_t *arraym = 0;
1141 Double_t *arrayme = 0;
1142 Float_t lowedge = 0.0;
1143 Float_t upedge = 0.0;
1144 // Init fCountDet and fCount
1145 InitfCountDetAndfCount(2);
1146 // Beginning of the loop
1147 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1148 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1149 UpdatefCountDetAndfCount(idect,2);
1150 ReconstructFitRowMinRowMax(idect,2);
1152 fEntriesCurrent = 0;
1153 if(!calvect->GetPRFEntries(fCountDet)) {
1154 NotEnoughStatisticPRF(idect);
1157 TString tname("PRF");
1159 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1160 nbins = projprftree->GetN();
1161 arrayx = (Double_t *)projprftree->GetX();
1162 arraye = (Double_t *)projprftree->GetEX();
1163 arraym = (Double_t *)projprftree->GetY();
1164 arrayme = (Double_t *)projprftree->GetEY();
1165 Float_t step = arrayx[1]-arrayx[0];
1166 lowedge = arrayx[0] - step/2.0;
1167 upedge = arrayx[(nbins-1)] + step/2.0;
1168 //printf("nbins est %d\n",nbins);
1169 for(Int_t k = 0; k < nbins; k++){
1170 fEntriesCurrent += (Int_t)arraye[k];
1171 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1172 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1174 if(fEntriesCurrent > 0) fNumberEnt++;
1175 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1176 // This detector has not enough statistics or was off
1177 if (fEntriesCurrent <= fMinEntries) {
1178 NotEnoughStatisticPRF(idect);
1181 // Statistic of the histos fitted
1183 fStatisticMean += fEntriesCurrent;
1184 // Calcul of "real" coef
1185 CalculPRFCoefMean();
1189 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1190 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1191 default: return kFALSE;
1193 // Fill the tree if end of a detector or only the pointer to the branch!!!
1194 FillInfosFitPRF(idect);
1197 if (fNumberFit > 0) {
1198 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));
1201 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1203 delete fDebugStreamer;
1204 fDebugStreamer = 0x0;
1207 //____________Functions fit Online CH2d________________________________________
1208 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1211 // The linear method
1214 fStatisticMean = 0.0;
1216 fNumberFitSuccess = 0;
1218 if(!InitFitLinearFitter()) return kFALSE;
1221 for(Int_t idet = 0; idet < 540; idet++){
1224 //printf("detector number %d\n",idet);
1229 fEntriesCurrent = 0;
1231 Bool_t here = calivdli->GetParam(idet,¶m);
1232 Bool_t heree = calivdli->GetError(idet,&error);
1233 //printf("here %d and heree %d\n",here, heree);
1235 fEntriesCurrent = (Int_t) error[2];
1238 //printf("Number of entries %d\n",fEntriesCurrent);
1239 // Nothing found or not enough statistic
1240 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1241 NotEnoughStatisticLinearFitter();
1248 fStatisticMean += fEntriesCurrent;
1251 if((-(param[1])) <= 0.0) {
1252 NotEnoughStatisticLinearFitter();
1256 // CalculDatabaseVdriftandTan
1257 CalculVdriftLorentzCoef();
1258 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFit detector %d, vdrift %f and %f and exB %f and %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCurrentCoef[1],fCalDetExBUsed->GetValue(idet),fCurrentCoef2[1]);
1261 fNumberFitSuccess ++;
1263 // Put the fCurrentCoef
1264 fCurrentCoef[0] = -param[1];
1265 // here the database must be the one of the reconstruction for the lorentz angle....
1266 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1267 fCurrentCoefE = error[1];
1268 fCurrentCoefE2 = error[0];
1269 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1270 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1274 FillInfosFitLinearFitter();
1279 if (fNumberFit > 0) {
1280 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));
1283 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1285 delete fDebugStreamer;
1286 fDebugStreamer = 0x0;
1290 //______________________________________________________________________________________
1291 Bool_t AliTRDCalibraFit::AnalyseExbAltFit(AliTRDCalibraExbAltFit *calivdli)
1294 // The linear method
1297 fStatisticMean = 0.0;
1299 fNumberFitSuccess = 0;
1301 if(!InitFitExbAlt()) return kFALSE;
1304 for(Int_t idet = 0; idet < 540; idet++){
1307 //printf("detector number %d\n",idet);
1312 fEntriesCurrent = 0;
1314 Bool_t here = calivdli->GetParam(idet,¶m);
1315 Bool_t heree = calivdli->GetError(idet,&error);
1316 //printf("here %d and heree %d\n",here, heree);
1318 fEntriesCurrent = (Int_t) error[2];
1321 //printf("Number of entries %d\n",fEntriesCurrent);
1322 // Nothing found or not enough statistic
1323 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1324 NotEnoughStatisticExbAlt();
1331 fStatisticMean += fEntriesCurrent;
1334 fNumberFitSuccess ++;
1336 // Put the fCurrentCoef
1337 if(TMath::Abs(param[2])>0.0001){
1338 fCurrentCoef2[0] = -param[1]/2/param[2];
1339 fCurrentCoefE2 = 0;//error[1];
1341 fCurrentCoef2[0] = 100;
1342 fCurrentCoefE2 = 0;//error[1];
1346 FillInfosFitExbAlt();
1350 if (fNumberFit > 0) {
1351 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));
1354 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1356 delete fDebugStreamer;
1357 fDebugStreamer = 0x0;
1361 //____________Functions fit Online CH2d________________________________________
1362 void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall)
1365 // The linear method
1368 // Get the mean vdrift and exb used
1369 Double_t meanvdriftused = 0.0;
1370 Double_t meanexbused = 0.0;
1371 Double_t counterdet = 0.0;
1372 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) {
1373 vdriftoverall = -100.0;
1380 TH2S *linearfitterhisto = 0x0;
1382 for(Int_t idet = 0; idet < 540; idet++){
1384 TH2S * u = calivdli->GetLinearFitterHistoForce(idet);
1385 Double_t detectorentries = u->Integral();
1386 meanvdriftused += fCalDetVdriftUsed->GetValue(idet)*detectorentries;
1387 meanexbused += fCalDetExBUsed->GetValue(idet)*detectorentries;
1388 counterdet += detectorentries;
1390 //printf("detectorentries %f\n",detectorentries);
1392 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether detector %d, vdrift %f and exB %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCalDetExBUsed->GetValue(idet));
1394 if(idet == 0) linearfitterhisto = u;
1395 else linearfitterhisto->Add(u);
1398 if(counterdet > 0.0){
1399 meanvdriftused = meanvdriftused/counterdet;
1400 meanexbused = meanexbused/counterdet;
1403 vdriftoverall = -100.0;
1409 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether MEAN vdrift %f and exB %f\n",meanvdriftused,meanexbused);
1414 TAxis *xaxis = linearfitterhisto->GetXaxis();
1415 TAxis *yaxis = linearfitterhisto->GetYaxis();
1416 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
1418 Double_t integral = linearfitterhisto->Integral();
1419 //printf("Integral is %f\n",integral);
1420 Bool_t securitybreaking = kFALSE;
1421 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
1422 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
1423 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
1424 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
1425 Double_t x = xaxis->GetBinCenter(ibinx+1);
1426 Double_t y = yaxis->GetBinCenter(ibiny+1);
1428 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
1429 if(!securitybreaking){
1430 linearfitter.AddPoint(&x,y);
1435 linearfitter.AddPoint(&x,y);
1445 //printf("AnalyseLinearFittersAllTogether::Find %d entries\n",entries);
1446 //printf("Minstats %d\n",fMinEntries);
1450 // Eval the linear fitter
1451 if(entries > fMinEntries){
1452 TVectorD par = TVectorD(2);
1454 if((linearfitter.EvalRobust(0.8)==0)) {
1455 //printf("Take the param\n");
1456 linearfitter.GetParameters(par);
1459 //printf("Finish\n");
1460 // Put the fCurrentCoef
1461 fCurrentCoef[0] = -par[1];
1462 // here the database must be the one of the reconstruction for the lorentz angle....
1463 if(fCurrentCoef[0] > 0.0) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
1464 else fCurrentCoef2[0] = 100.0;
1469 fCurrentCoef[0] = -100.0;
1470 fCurrentCoef2[0] = 100.0;
1478 fCurrentCoef[0] = -100.0;
1479 fCurrentCoef2[0] = 100.0;
1483 vdriftoverall = fCurrentCoef[0];
1484 exboverall = fCurrentCoef2[0];
1487 delete linearfitterhisto;
1488 delete fDebugStreamer;
1489 fDebugStreamer = 0x0;
1492 //____________Functions for seeing if the pad is really okey___________________
1493 //_____________________________________________________________________________
1494 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1497 // Get numberofgroupsprf
1501 const Char_t *pattern0 = "Ngp0";
1502 const Char_t *pattern1 = "Ngp1";
1503 const Char_t *pattern2 = "Ngp2";
1504 const Char_t *pattern3 = "Ngp3";
1505 const Char_t *pattern4 = "Ngp4";
1506 const Char_t *pattern5 = "Ngp5";
1507 const Char_t *pattern6 = "Ngp6";
1510 if (strstr(nametitle.Data(),pattern0)) {
1513 if (strstr(nametitle.Data(),pattern1)) {
1516 if (strstr(nametitle.Data(),pattern2)) {
1519 if (strstr(nametitle.Data(),pattern3)) {
1522 if (strstr(nametitle.Data(),pattern4)) {
1525 if (strstr(nametitle.Data(),pattern5)) {
1528 if (strstr(nametitle.Data(),pattern6)){
1535 //_____________________________________________________________________________
1536 Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1539 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1540 // corresponding to the given name
1543 if(!SetNzFromTObject(name,i)) return kFALSE;
1544 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1549 //_____________________________________________________________________________
1550 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1553 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1554 // corresponding to the given TObject
1558 const Char_t *patternrphi0 = "Nrphi0";
1559 const Char_t *patternrphi1 = "Nrphi1";
1560 const Char_t *patternrphi2 = "Nrphi2";
1561 const Char_t *patternrphi3 = "Nrphi3";
1562 const Char_t *patternrphi4 = "Nrphi4";
1563 const Char_t *patternrphi5 = "Nrphi5";
1564 const Char_t *patternrphi6 = "Nrphi6";
1567 const Char_t *patternrphi10 = "Nrphi10";
1568 const Char_t *patternrphi100 = "Nrphi100";
1569 const Char_t *patternz10 = "Nz10";
1570 const Char_t *patternz100 = "Nz100";
1573 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1574 fCalibraMode->SetAllTogether(i);
1576 if (fDebugLevel > 1) {
1577 AliInfo(Form("fNbDet %d and 100",fNbDet));
1581 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1582 fCalibraMode->SetPerSuperModule(i);
1584 if (fDebugLevel > 1) {
1585 AliInfo(Form("fNDet %d and 100",fNbDet));
1590 if (strstr(name.Data(),patternrphi0)) {
1591 fCalibraMode->SetNrphi(i ,0);
1592 if (fDebugLevel > 1) {
1593 AliInfo(Form("fNbDet %d and 0",fNbDet));
1597 if (strstr(name.Data(),patternrphi1)) {
1598 fCalibraMode->SetNrphi(i, 1);
1599 if (fDebugLevel > 1) {
1600 AliInfo(Form("fNbDet %d and 1",fNbDet));
1604 if (strstr(name.Data(),patternrphi2)) {
1605 fCalibraMode->SetNrphi(i, 2);
1606 if (fDebugLevel > 1) {
1607 AliInfo(Form("fNbDet %d and 2",fNbDet));
1611 if (strstr(name.Data(),patternrphi3)) {
1612 fCalibraMode->SetNrphi(i, 3);
1613 if (fDebugLevel > 1) {
1614 AliInfo(Form("fNbDet %d and 3",fNbDet));
1618 if (strstr(name.Data(),patternrphi4)) {
1619 fCalibraMode->SetNrphi(i, 4);
1620 if (fDebugLevel > 1) {
1621 AliInfo(Form("fNbDet %d and 4",fNbDet));
1625 if (strstr(name.Data(),patternrphi5)) {
1626 fCalibraMode->SetNrphi(i, 5);
1627 if (fDebugLevel > 1) {
1628 AliInfo(Form("fNbDet %d and 5",fNbDet));
1632 if (strstr(name.Data(),patternrphi6)) {
1633 fCalibraMode->SetNrphi(i, 6);
1634 if (fDebugLevel > 1) {
1635 AliInfo(Form("fNbDet %d and 6",fNbDet));
1640 if (fDebugLevel > 1) {
1641 AliInfo(Form("fNbDet %d and rest",fNbDet));
1643 fCalibraMode->SetNrphi(i ,0);
1647 //_____________________________________________________________________________
1648 Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1651 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1652 // corresponding to the given TObject
1656 const Char_t *patternz0 = "Nz0";
1657 const Char_t *patternz1 = "Nz1";
1658 const Char_t *patternz2 = "Nz2";
1659 const Char_t *patternz3 = "Nz3";
1660 const Char_t *patternz4 = "Nz4";
1662 const Char_t *patternrphi10 = "Nrphi10";
1663 const Char_t *patternrphi100 = "Nrphi100";
1664 const Char_t *patternz10 = "Nz10";
1665 const Char_t *patternz100 = "Nz100";
1667 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1668 fCalibraMode->SetAllTogether(i);
1670 if (fDebugLevel > 1) {
1671 AliInfo(Form("fNbDet %d and 100",fNbDet));
1675 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1676 fCalibraMode->SetPerSuperModule(i);
1678 if (fDebugLevel > 1) {
1679 AliInfo(Form("fNbDet %d and 10",fNbDet));
1683 if (strstr(name.Data(),patternz0)) {
1684 fCalibraMode->SetNz(i, 0);
1685 if (fDebugLevel > 1) {
1686 AliInfo(Form("fNbDet %d and 0",fNbDet));
1690 if (strstr(name.Data(),patternz1)) {
1691 fCalibraMode->SetNz(i ,1);
1692 if (fDebugLevel > 1) {
1693 AliInfo(Form("fNbDet %d and 1",fNbDet));
1697 if (strstr(name.Data(),patternz2)) {
1698 fCalibraMode->SetNz(i ,2);
1699 if (fDebugLevel > 1) {
1700 AliInfo(Form("fNbDet %d and 2",fNbDet));
1704 if (strstr(name.Data(),patternz3)) {
1705 fCalibraMode->SetNz(i ,3);
1706 if (fDebugLevel > 1) {
1707 AliInfo(Form("fNbDet %d and 3",fNbDet));
1711 if (strstr(name.Data(),patternz4)) {
1712 fCalibraMode->SetNz(i ,4);
1713 if (fDebugLevel > 1) {
1714 AliInfo(Form("fNbDet %d and 4",fNbDet));
1719 if (fDebugLevel > 1) {
1720 AliInfo(Form("fNbDet %d and rest",fNbDet));
1722 fCalibraMode->SetNz(i ,0);
1725 //______________________________________________________________________
1726 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1728 // Remove the results too far from the mean value and rms
1729 // type: 0 gain, 1 vdrift
1733 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1735 AliInfo("The Vector Fit is not complete!");
1738 Int_t detector = -1;
1740 Float_t value = 0.0;
1742 /////////////////////////////////
1743 // Calculate the mean values
1744 ////////////////////////////////
1746 ////////////////////////
1747 Double_t meanAll = 0.0;
1748 Double_t rmsAll = 0.0;
1753 for (Int_t k = 0; k < loop; k++) {
1754 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1755 sector = GetSector(detector);
1757 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1759 rmsAll += value*value;
1765 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1766 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1767 for (Int_t row = 0; row < rowMax; row++) {
1768 for (Int_t col = 0; col < colMax; col++) {
1769 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1771 rmsAll += value*value;
1781 meanAll = meanAll/countAll;
1782 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1784 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1785 /////////////////////////////////////////////////
1787 ////////////////////////////////////////////////
1788 Double_t defaultvalue = -1.0;
1789 if(type==1) defaultvalue = -1.5;
1790 for (Int_t k = 0; k < loop; k++) {
1791 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1792 sector = GetSector(detector);
1793 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1794 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1795 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1797 // remove the results too far away
1798 for (Int_t row = 0; row < rowMax; row++) {
1799 for (Int_t col = 0; col < colMax; col++) {
1800 value = coef[(Int_t)(col*rowMax+row)];
1801 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1802 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1808 //______________________________________________________________________
1809 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1811 // Remove the results too far from the mean and rms
1815 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1817 AliInfo("The Vector Fit is not complete!");
1820 Int_t detector = -1;
1822 Float_t value = 0.0;
1824 /////////////////////////////////
1825 // Calculate the mean values
1826 ////////////////////////////////
1828 ////////////////////////
1829 Double_t meanAll = 0.0;
1830 Double_t rmsAll = 0.0;
1835 for (Int_t k = 0; k < loop; k++) {
1836 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1837 sector = GetSector(detector);
1839 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1842 rmsAll += value*value;
1847 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1848 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1849 for (Int_t row = 0; row < rowMax; row++) {
1850 for (Int_t col = 0; col < colMax; col++) {
1851 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1853 rmsAll += value*value;
1862 meanAll = meanAll/countAll;
1863 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1865 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1866 /////////////////////////////////////////////////
1868 ////////////////////////////////////////////////
1869 for (Int_t k = 0; k < loop; k++) {
1870 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1871 sector = GetSector(detector);
1872 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1873 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1874 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1876 // remove the results too far away
1877 for (Int_t row = 0; row < rowMax; row++) {
1878 for (Int_t col = 0; col < colMax; col++) {
1879 value = coef[(Int_t)(col*rowMax+row)];
1880 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) {
1881 //printf("value outlier %f\n",value);
1882 coef[(Int_t)(col*rowMax+row)] = 100.0;
1888 //______________________________________________________________________
1889 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1891 // ofwhat is equaled to 0: mean value of all passing detectors
1892 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1895 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1897 AliInfo("The Vector Fit is not complete!");
1900 Int_t detector = -1;
1902 Float_t value = 0.0;
1904 /////////////////////////////////
1905 // Calculate the mean values
1906 ////////////////////////////////
1908 ////////////////////////
1909 Double_t meanAll = 0.0;
1910 Double_t meanSupermodule[18];
1911 Double_t meanDetector[540];
1912 Double_t rmsAll = 0.0;
1913 Double_t rmsSupermodule[18];
1914 Double_t rmsDetector[540];
1916 Int_t countSupermodule[18];
1917 Int_t countDetector[540];
1918 for(Int_t sm = 0; sm < 18; sm++){
1919 rmsSupermodule[sm] = 0.0;
1920 meanSupermodule[sm] = 0.0;
1921 countSupermodule[sm] = 0;
1923 for(Int_t det = 0; det < 540; det++){
1924 rmsDetector[det] = 0.0;
1925 meanDetector[det] = 0.0;
1926 countDetector[det] = 0;
1931 for (Int_t k = 0; k < loop; k++) {
1932 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1933 sector = GetSector(detector);
1935 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1937 rmsDetector[detector] += value*value;
1938 meanDetector[detector] += value;
1939 countDetector[detector]++;
1940 rmsSupermodule[sector] += value*value;
1941 meanSupermodule[sector] += value;
1942 countSupermodule[sector]++;
1943 rmsAll += value*value;
1949 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1950 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1951 for (Int_t row = 0; row < rowMax; row++) {
1952 for (Int_t col = 0; col < colMax; col++) {
1953 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1955 rmsDetector[detector] += value*value;
1956 meanDetector[detector] += value;
1957 countDetector[detector]++;
1958 rmsSupermodule[sector] += value*value;
1959 meanSupermodule[sector] += value;
1960 countSupermodule[sector]++;
1961 rmsAll += value*value;
1971 meanAll = meanAll/countAll;
1972 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1974 for(Int_t sm = 0; sm < 18; sm++){
1975 if(countSupermodule[sm] > 0) {
1976 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1977 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1980 for(Int_t det = 0; det < 540; det++){
1981 if(countDetector[det] > 0) {
1982 meanDetector[det] = meanDetector[det]/countDetector[det];
1983 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1986 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1987 ///////////////////////////////////////////////
1988 // Put the mean value for the no-fitted
1989 /////////////////////////////////////////////
1990 for (Int_t k = 0; k < loop; k++) {
1991 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1992 sector = GetSector(detector);
1993 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1994 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1995 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1997 for (Int_t row = 0; row < rowMax; row++) {
1998 for (Int_t col = 0; col < colMax; col++) {
1999 value = coef[(Int_t)(col*rowMax+row)];
2001 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
2003 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
2004 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
2005 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
2009 if(fDebugLevel > 1){
2011 if ( !fDebugStreamer ) {
2013 TDirectory *backup = gDirectory;
2014 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2015 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2018 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2020 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
2021 "detector="<<detector<<
2033 //______________________________________________________________________
2034 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
2036 // ofwhat is equaled to 0: mean value of all passing detectors
2037 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
2040 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
2042 AliInfo("The Vector Fit is not complete!");
2045 Int_t detector = -1;
2047 Float_t value = 0.0;
2049 /////////////////////////////////
2050 // Calculate the mean values
2051 ////////////////////////////////
2053 ////////////////////////
2054 Double_t meanAll = 0.0;
2055 Double_t rmsAll = 0.0;
2056 Double_t meanSupermodule[18];
2057 Double_t rmsSupermodule[18];
2058 Double_t meanDetector[540];
2059 Double_t rmsDetector[540];
2061 Int_t countSupermodule[18];
2062 Int_t countDetector[540];
2063 for(Int_t sm = 0; sm < 18; sm++){
2064 rmsSupermodule[sm] = 0.0;
2065 meanSupermodule[sm] = 0.0;
2066 countSupermodule[sm] = 0;
2068 for(Int_t det = 0; det < 540; det++){
2069 rmsDetector[det] = 0.0;
2070 meanDetector[det] = 0.0;
2071 countDetector[det] = 0;
2075 for (Int_t k = 0; k < loop; k++) {
2076 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2077 sector = GetSector(detector);
2079 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
2081 rmsDetector[detector] += value*value;
2082 meanDetector[detector] += value;
2083 countDetector[detector]++;
2084 rmsSupermodule[sector] += value*value;
2085 meanSupermodule[sector] += value;
2086 countSupermodule[sector]++;
2088 rmsAll += value*value;
2093 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2094 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2095 for (Int_t row = 0; row < rowMax; row++) {
2096 for (Int_t col = 0; col < colMax; col++) {
2097 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2099 rmsDetector[detector] += value*value;
2100 meanDetector[detector] += value;
2101 countDetector[detector]++;
2102 rmsSupermodule[sector] += value*value;
2103 meanSupermodule[sector] += value;
2104 countSupermodule[sector]++;
2105 rmsAll += value*value;
2115 meanAll = meanAll/countAll;
2116 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
2118 for(Int_t sm = 0; sm < 18; sm++){
2119 if(countSupermodule[sm] > 0) {
2120 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
2121 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
2124 for(Int_t det = 0; det < 540; det++){
2125 if(countDetector[det] > 0) {
2126 meanDetector[det] = meanDetector[det]/countDetector[det];
2127 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
2130 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
2131 ////////////////////////////////////////////
2132 // Put the mean value for the no-fitted
2133 /////////////////////////////////////////////
2134 for (Int_t k = 0; k < loop; k++) {
2135 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2136 sector = GetSector(detector);
2137 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2138 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2139 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
2141 for (Int_t row = 0; row < rowMax; row++) {
2142 for (Int_t col = 0; col < colMax; col++) {
2143 value = coef[(Int_t)(col*rowMax+row)];
2145 if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2147 if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
2148 else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
2149 else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2153 if(fDebugLevel > 1){
2155 if ( !fDebugStreamer ) {
2157 TDirectory *backup = gDirectory;
2158 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2159 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2162 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2164 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
2165 "detector="<<detector<<
2178 //_____________________________________________________________________________
2179 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
2182 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2183 // It takes the mean value of the coefficients per detector
2184 // This object has to be written in the database
2187 // Create the DetObject
2188 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2190 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2191 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2192 Int_t detector = -1;
2193 Float_t value = 0.0;
2196 for (Int_t k = 0; k < loop; k++) {
2197 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2200 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
2204 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2205 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2206 for (Int_t row = 0; row < rowMax; row++) {
2207 for (Int_t col = 0; col < colMax; col++) {
2208 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2209 mean += TMath::Abs(value);
2213 if(count > 0) mean = mean/count;
2215 object->SetValue(detector,mean);
2220 //_____________________________________________________________________________
2221 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
2224 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2225 // It takes the mean value of the coefficients per detector
2226 // This object has to be written in the database
2229 // Create the DetObject
2230 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2232 fScaleGain = scaleFitFactor;
2234 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2235 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2236 Int_t detector = -1;
2237 Float_t value = 0.0;
2239 for (Int_t k = 0; k < loop; k++) {
2240 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2243 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2244 if(!meanOtherBefore){
2245 if(value > 0) value = value*scaleFitFactor;
2247 else value = value*scaleFitFactor;
2248 mean = TMath::Abs(value);
2252 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2253 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2254 for (Int_t row = 0; row < rowMax; row++) {
2255 for (Int_t col = 0; col < colMax; col++) {
2256 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2257 if(!meanOtherBefore) {
2258 if(value > 0) value = value*scaleFitFactor;
2260 else value = value*scaleFitFactor;
2261 mean += TMath::Abs(value);
2265 if(count > 0) mean = mean/count;
2267 if(mean < 0.1) mean = 0.1;
2268 object->SetValue(detector,mean);
2273 //_____________________________________________________________________________
2274 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
2277 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2278 // It takes the min value of the coefficients per detector
2279 // This object has to be written in the database
2282 // Create the DetObject
2283 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2285 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2286 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2287 Int_t detector = -1;
2288 Float_t value = 0.0;
2290 for (Int_t k = 0; k < loop; k++) {
2291 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2292 Float_t min = 100.0;
2294 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2295 //printf("Create det object %f for %d\n",value,k);
2297 if(value > 70.0) value = value-100.0;
2302 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2303 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2304 for (Int_t row = 0; row < rowMax; row++) {
2305 for (Int_t col = 0; col < colMax; col++) {
2306 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2308 if(value > 70.0) value = value-100.0;
2310 if(min > value) min = value;
2314 object->SetValue(detector,min);
2320 //_____________________________________________________________________________
2321 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2324 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2325 // It takes the min value of the coefficients per detector
2326 // This object has to be written in the database
2329 // Create the DetObject
2330 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2333 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2334 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2335 Int_t detector = -1;
2336 Float_t value = 0.0;
2338 for (Int_t k = 0; k < loop; k++) {
2339 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2341 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2342 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2343 Float_t min = 100.0;
2344 for (Int_t row = 0; row < rowMax; row++) {
2345 for (Int_t col = 0; col < colMax; col++) {
2346 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2347 mean += -TMath::Abs(value);
2351 if(count > 0) mean = mean/count;
2353 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2354 if(value > 70.0) value = value-100.0;
2355 object->SetValue(detector,value);
2361 //_____________________________________________________________________________
2362 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit)
2365 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2366 // It takes the min value of the coefficients per detector
2367 // This object has to be written in the database
2370 // Create the DetObject
2371 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2374 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2375 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2376 Int_t detector = -1;
2377 Float_t value = 0.0;
2379 for (Int_t k = 0; k < loop; k++) {
2380 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2382 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2383 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2384 Float_t min = 100.0;
2385 for (Int_t row = 0; row < rowMax; row++) {
2386 for (Int_t col = 0; col < colMax; col++) {
2387 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2388 mean += -TMath::Abs(value);
2392 if(count > 0) mean = mean/count;
2394 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2395 //if(value > 70.0) value = value-100.0;
2396 object->SetValue(detector,value);
2402 //_____________________________________________________________________________
2403 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2406 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2407 // You need first to create the object for the detectors,
2408 // where the mean value is put.
2409 // This object has to be written in the database
2412 // Create the DetObject
2413 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2416 for(Int_t k = 0; k < 540; k++){
2417 AliTRDCalROC *calROC = object->GetCalROC(k);
2418 Int_t nchannels = calROC->GetNchannels();
2419 for(Int_t ch = 0; ch < nchannels; ch++){
2420 calROC->SetValue(ch,1.0);
2426 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2427 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2428 Int_t detector = -1;
2429 Float_t value = 0.0;
2431 for (Int_t k = 0; k < loop; k++) {
2432 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2433 AliTRDCalROC *calROC = object->GetCalROC(detector);
2434 Float_t mean = detobject->GetValue(detector);
2435 if(TMath::Abs(mean) <= 0.0000000001) continue;
2436 Int_t rowMax = calROC->GetNrows();
2437 Int_t colMax = calROC->GetNcols();
2438 for (Int_t row = 0; row < rowMax; row++) {
2439 for (Int_t col = 0; col < colMax; col++) {
2440 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2441 if(value > 0) value = value*scaleFitFactor;
2442 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2450 //_____________________________________________________________________________
2451 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2454 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2455 // You need first to create the object for the detectors,
2456 // where the mean value is put.
2457 // This object has to be written in the database
2460 // Create the DetObject
2461 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2464 for(Int_t k = 0; k < 540; k++){
2465 AliTRDCalROC *calROC = object->GetCalROC(k);
2466 Int_t nchannels = calROC->GetNchannels();
2467 for(Int_t ch = 0; ch < nchannels; ch++){
2468 calROC->SetValue(ch,1.0);
2474 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2475 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2476 Int_t detector = -1;
2477 Float_t value = 0.0;
2479 for (Int_t k = 0; k < loop; k++) {
2480 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2481 AliTRDCalROC *calROC = object->GetCalROC(detector);
2482 Float_t mean = detobject->GetValue(detector);
2483 if(mean == 0) continue;
2484 Int_t rowMax = calROC->GetNrows();
2485 Int_t colMax = calROC->GetNcols();
2486 for (Int_t row = 0; row < rowMax; row++) {
2487 for (Int_t col = 0; col < colMax; col++) {
2488 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2489 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2497 //_____________________________________________________________________________
2498 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2501 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2502 // You need first to create the object for the detectors,
2503 // where the mean value is put.
2504 // This object has to be written in the database
2507 // Create the DetObject
2508 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2511 for(Int_t k = 0; k < 540; k++){
2512 AliTRDCalROC *calROC = object->GetCalROC(k);
2513 Int_t nchannels = calROC->GetNchannels();
2514 for(Int_t ch = 0; ch < nchannels; ch++){
2515 calROC->SetValue(ch,0.0);
2521 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2522 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2523 Int_t detector = -1;
2524 Float_t value = 0.0;
2526 for (Int_t k = 0; k < loop; k++) {
2527 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2528 AliTRDCalROC *calROC = object->GetCalROC(detector);
2529 Float_t min = detobject->GetValue(detector);
2530 Int_t rowMax = calROC->GetNrows();
2531 Int_t colMax = calROC->GetNcols();
2532 for (Int_t row = 0; row < rowMax; row++) {
2533 for (Int_t col = 0; col < colMax; col++) {
2534 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2536 if(value > 70.0) value = value - 100.0;
2538 calROC->SetValue(col,row,value-min);
2546 //_____________________________________________________________________________
2547 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2550 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2551 // This object has to be written in the database
2554 // Create the DetObject
2555 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2557 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2558 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2559 Int_t detector = -1;
2560 Float_t value = 0.0;
2562 for (Int_t k = 0; k < loop; k++) {
2563 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2564 AliTRDCalROC *calROC = object->GetCalROC(detector);
2565 Int_t rowMax = calROC->GetNrows();
2566 Int_t colMax = calROC->GetNcols();
2567 for (Int_t row = 0; row < rowMax; row++) {
2568 for (Int_t col = 0; col < colMax; col++) {
2569 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2570 calROC->SetValue(col,row,TMath::Abs(value));
2578 //_____________________________________________________________________________
2579 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2582 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2583 // 0 successful fit 1 not successful fit
2584 // mean is the mean value over the successful fit
2585 // do not use it for t0: no meaning
2588 // Create the CalObject
2589 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2593 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2595 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2596 for(Int_t k = 0; k < 540; k++){
2597 object->SetValue(k,1.0);
2600 Int_t detector = -1;
2601 Float_t value = 0.0;
2603 for (Int_t k = 0; k < loop; k++) {
2604 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2605 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2606 if(value <= 0) object->SetValue(detector,1.0);
2608 object->SetValue(detector,0.0);
2613 if(count > 0) mean /= count;
2616 //_____________________________________________________________________________
2617 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2620 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2621 // 0 not successful fit 1 successful fit
2622 // mean mean value over the successful fit
2625 // Create the CalObject
2626 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2630 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2632 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2633 for(Int_t k = 0; k < 540; k++){
2634 AliTRDCalROC *calROC = object->GetCalROC(k);
2635 Int_t nchannels = calROC->GetNchannels();
2636 for(Int_t ch = 0; ch < nchannels; ch++){
2637 calROC->SetValue(ch,1.0);
2641 Int_t detector = -1;
2642 Float_t value = 0.0;
2644 for (Int_t k = 0; k < loop; k++) {
2645 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2646 AliTRDCalROC *calROC = object->GetCalROC(detector);
2647 Int_t nchannels = calROC->GetNchannels();
2648 for (Int_t ch = 0; ch < nchannels; ch++) {
2649 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2650 if(value <= 0) calROC->SetValue(ch,1.0);
2652 calROC->SetValue(ch,0.0);
2658 if(count > 0) mean /= count;
2661 //_____________________________________________________________________________
2662 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2665 // Set FitPH if 1 then each detector will be fitted
2668 if (periodeFitPH > 0) {
2669 fFitPHPeriode = periodeFitPH;
2672 AliInfo("periodeFitPH must be higher than 0!");
2676 //_____________________________________________________________________________
2677 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2680 // The fit of the deposited charge distribution begins at
2681 // histo->Mean()/beginFitCharge
2682 // You can here set beginFitCharge
2685 if (beginFitCharge > 0) {
2686 fBeginFitCharge = beginFitCharge;
2689 AliInfo("beginFitCharge must be strict positif!");
2694 //_____________________________________________________________________________
2695 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2698 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2699 // You can here set t0Shift0
2703 fT0Shift0 = t0Shift;
2706 AliInfo("t0Shift0 must be strict positif!");
2711 //_____________________________________________________________________________
2712 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2715 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2716 // You can here set t0Shift1
2720 fT0Shift1 = t0Shift;
2723 AliInfo("t0Shift must be strict positif!");
2728 //_____________________________________________________________________________
2729 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2732 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2733 // You can here set rangeFitPRF
2736 if ((rangeFitPRF > 0) &&
2737 (rangeFitPRF <= 1.5)) {
2738 fRangeFitPRF = rangeFitPRF;
2741 AliInfo("rangeFitPRF must be between 0 and 1.0");
2746 //_____________________________________________________________________________
2747 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2750 // Minimum entries for fitting
2753 if (minEntries > 0) {
2754 fMinEntries = minEntries;
2757 AliInfo("fMinEntries must be >= 0.");
2762 //_____________________________________________________________________________
2763 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2766 // Rebin with rebin time less bins the Ch histo
2767 // You can set here rebin that should divide the number of bins of CH histo
2772 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2775 AliInfo("You have to choose a positiv value!");
2779 //_____________________________________________________________________________
2780 Bool_t AliTRDCalibraFit::FillVectorFit()
2783 // For the Fit functions fill the vector Fit
2786 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2789 if (GetStack(fCountDet) == 2) {
2796 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2797 Float_t *coef = new Float_t[ntotal];
2798 for (Int_t i = 0; i < ntotal; i++) {
2799 coef[i] = fCurrentCoefDetector[i];
2802 Int_t detector = fCountDet;
2804 fitInfo->SetCoef(coef);
2805 fitInfo->SetDetector(detector);
2806 fVectorFit.Add((TObject *) fitInfo);
2811 //_____________________________________________________________________________
2812 Bool_t AliTRDCalibraFit::FillVectorFit2()
2815 // For the Fit functions fill the vector Fit
2818 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2821 if (GetStack(fCountDet) == 2) {
2828 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2829 Float_t *coef = new Float_t[ntotal];
2830 for (Int_t i = 0; i < ntotal; i++) {
2831 coef[i] = fCurrentCoefDetector2[i];
2834 Int_t detector = fCountDet;
2836 fitInfo->SetCoef(coef);
2837 fitInfo->SetDetector(detector);
2838 fVectorFit2.Add((TObject *) fitInfo);
2843 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2844 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2847 // Init the number of expected bins and fDect1[i] fDect2[i]
2850 gStyle->SetPalette(1);
2851 gStyle->SetOptStat(1111);
2852 gStyle->SetPadBorderMode(0);
2853 gStyle->SetCanvasColor(10);
2854 gStyle->SetPadLeftMargin(0.13);
2855 gStyle->SetPadRightMargin(0.01);
2857 // Mode groups of pads: the total number of bins!
2858 CalculNumberOfBinsExpected(i);
2860 // Quick verification that we have the good pad calibration mode!
2861 if (fNumberOfBinsExpected != nbins) {
2862 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2866 // Security for fDebug 3 and 4
2867 if ((fDebugLevel >= 3) &&
2871 AliInfo("This detector doesn't exit!");
2875 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2876 CalculDect1Dect2(i);
2881 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2882 Bool_t AliTRDCalibraFit::InitFitCH()
2885 // Init the fVectorFitCH for normalisation
2886 // Init the histo for debugging
2891 fScaleFitFactor = 0.0;
2892 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2893 fCurrentCoefDetector = new Float_t[2304];
2894 for (Int_t k = 0; k < 2304; k++) {
2895 fCurrentCoefDetector[k] = 0.0;
2897 fVectorFit.SetName("gainfactorscoefficients");
2899 // fDebug == 0 nothing
2900 // fDebug == 1 and fFitVoir no histo
2901 if (fDebugLevel == 1) {
2902 if(!CheckFitVoir()) return kFALSE;
2904 //Get the CalDet object
2906 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2908 AliInfo("Could not get calibDB");
2911 if(fCalDet) delete fCalDet;
2912 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2915 Float_t devalue = 1.0;
2916 if(fCalDet) delete fCalDet;
2917 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2918 for(Int_t k = 0; k < 540; k++){
2919 fCalDet->SetValue(k,devalue);
2925 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2926 Bool_t AliTRDCalibraFit::InitFitPH()
2929 // Init the arrays of results
2930 // Init the histos for debugging
2934 fVectorFit.SetName("driftvelocitycoefficients");
2935 fVectorFit2.SetName("t0coefficients");
2937 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2938 fCurrentCoefDetector = new Float_t[2304];
2939 for (Int_t k = 0; k < 2304; k++) {
2940 fCurrentCoefDetector[k] = 0.0;
2942 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
2943 fCurrentCoefDetector2 = new Float_t[2304];
2944 for (Int_t k = 0; k < 2304; k++) {
2945 fCurrentCoefDetector2[k] = 0.0;
2948 //fDebug == 0 nothing
2949 // fDebug == 1 and fFitVoir no histo
2950 if (fDebugLevel == 1) {
2951 if(!CheckFitVoir()) return kFALSE;
2953 //Get the CalDet object
2955 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2957 AliInfo("Could not get calibDB");
2960 if(fCalDet) delete fCalDet;
2961 if(fCalDet2) delete fCalDet2;
2962 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2963 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2966 Float_t devalue = 1.5;
2967 Float_t devalue2 = 0.0;
2968 if(fCalDet) delete fCalDet;
2969 if(fCalDet2) delete fCalDet2;
2970 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2971 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2972 for(Int_t k = 0; k < 540; k++){
2973 fCalDet->SetValue(k,devalue);
2974 fCalDet2->SetValue(k,devalue2);
2979 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2980 Bool_t AliTRDCalibraFit::InitFitPRF()
2983 // Init the calibration mode (Nz, Nrphi), the histograms for
2984 // debugging the fit methods if fDebug > 0,
2988 fVectorFit.SetName("prfwidthcoefficients");
2990 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2991 fCurrentCoefDetector = new Float_t[2304];
2992 for (Int_t k = 0; k < 2304; k++) {
2993 fCurrentCoefDetector[k] = 0.0;
2996 // fDebug == 0 nothing
2997 // fDebug == 1 and fFitVoir no histo
2998 if (fDebugLevel == 1) {
2999 if(!CheckFitVoir()) return kFALSE;
3003 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3004 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
3007 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3012 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
3013 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3014 fCurrentCoefDetector = new Float_t[2304];
3015 fCurrentCoefDetector2 = new Float_t[2304];
3016 for (Int_t k = 0; k < 2304; k++) {
3017 fCurrentCoefDetector[k] = 0.0;
3018 fCurrentCoefDetector2[k] = 0.0;
3021 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) return kFALSE;
3025 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3026 Bool_t AliTRDCalibraFit::InitFitExbAlt()
3029 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3034 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3035 fCurrentCoefDetector2 = new Float_t[2304];
3036 for (Int_t k = 0; k < 2304; k++) {
3037 fCurrentCoefDetector2[k] = 0.0;
3042 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3043 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
3046 // Init the current detector where we are fCountDet and the
3047 // next fCount for the functions Fit...
3050 // Loop on the Xbins of ch!!
3051 fCountDet = -1; // Current detector
3052 fCount = 0; // To find the next detector
3055 if (fDebugLevel >= 3) {
3056 // Set countdet to the detector
3057 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3058 // Set counter to write at the end of the detector
3060 // Get the right calib objects
3063 if(fDebugLevel == 1) {
3065 fCalibraMode->CalculXBins(fCountDet,i);
3066 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
3067 while(fCalibraMode->GetXbins(i) <=fFitVoir){
3069 fCalibraMode->CalculXBins(fCountDet,i);
3070 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
3076 fCount = fCalibraMode->GetXbins(i);
3078 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3079 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3080 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3081 ,(Int_t) GetStack(fCountDet)
3082 ,(Int_t) GetSector(fCountDet),i);
3085 //_______________________________________________________________________________
3086 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
3089 // Calculate the number of bins expected (calibration groups)
3092 fNumberOfBinsExpected = 0;
3094 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
3095 fNumberOfBinsExpected = 1;
3099 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
3100 fNumberOfBinsExpected = 18;
3104 fCalibraMode->ModePadCalibration(2,i);
3105 fCalibraMode->ModePadFragmentation(0,2,0,i);
3106 fCalibraMode->SetDetChamb2(i);
3107 if (fDebugLevel > 1) {
3108 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
3110 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
3111 fCalibraMode->ModePadCalibration(0,i);
3112 fCalibraMode->ModePadFragmentation(0,0,0,i);
3113 fCalibraMode->SetDetChamb0(i);
3114 if (fDebugLevel > 1) {
3115 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
3117 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
3120 //_______________________________________________________________________________
3121 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
3124 // Calculate the range of fits
3129 if (fDebugLevel == 1) {
3133 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
3135 fDect2 = fNumberOfBinsExpected;
3137 if (fDebugLevel >= 3) {
3138 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3139 fCalibraMode->CalculXBins(fCountDet,i);
3140 fDect1 = fCalibraMode->GetXbins(i);
3141 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3142 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3143 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3144 ,(Int_t) GetStack(fCountDet)
3145 ,(Int_t) GetSector(fCountDet),i);
3146 // Set for the next detector
3147 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3150 //_______________________________________________________________________________
3151 Bool_t AliTRDCalibraFit::CheckFitVoir()
3154 // Check if fFitVoir is in the range
3157 if (fFitVoir < fNumberOfBinsExpected) {
3158 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
3161 AliInfo("fFitVoir is out of range of the histo!");
3166 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3167 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
3170 // See if we are in a new detector and update the
3171 // variables fNfragZ and fNfragRphi if yes
3172 // Will never happen for only one detector (3 and 4)
3173 // Doesn't matter for 2
3175 if (fCount == idect) {
3176 // On en est au detector (or first detector in the group)
3178 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
3179 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3180 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3181 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3182 ,(Int_t) GetStack(fCountDet)
3183 ,(Int_t) GetSector(fCountDet),i);
3184 // Set for the next detector
3185 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3190 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3191 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
3194 // Reconstruct the min pad row, max pad row, min pad col and
3195 // max pad col of the calibration group for the Fit functions
3196 // idect is the calibration group inside the detector
3198 if (fDebugLevel != 1) {
3199 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
3201 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
3202 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
3204 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3205 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
3208 // For the case where there are not enough entries in the histograms
3209 // of the calibration group, the value present in the choosen database
3210 // will be put. A negativ sign enables to know that a fit was not possible.
3213 if (fDebugLevel == 1) {
3214 AliInfo("The element has not enough statistic to be fitted");
3216 else if (fNbDet > 0){
3217 Int_t firstdetector = fCountDet;
3218 Int_t lastdetector = fCountDet+fNbDet;
3219 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3220 // loop over detectors
3221 for(Int_t det = firstdetector; det < lastdetector; det++){
3223 //Set the calibration object again
3227 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3229 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3230 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3231 ,(Int_t) GetStack(fCountDet)
3232 ,(Int_t) GetSector(fCountDet),0);
3233 // Reconstruct row min row max
3234 ReconstructFitRowMinRowMax(idect,0);
3236 // Calcul the coef from the database choosen for the detector
3237 CalculChargeCoefMean(kFALSE);
3239 //stack 2, not stack 2
3241 if(GetStack(fCountDet) == 2) factor = 12;
3244 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3245 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3246 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3247 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3251 //Put default value negative
3252 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3253 fCurrentCoefE = 0.0;
3258 if(fDebugLevel > 1){
3260 if ( !fDebugStreamer ) {
3262 TDirectory *backup = gDirectory;
3263 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3264 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3267 Int_t detector = fCountDet;
3268 Int_t caligroup = idect;
3269 Short_t rowmin = fCalibraMode->GetRowMin(0);
3270 Short_t rowmax = fCalibraMode->GetRowMax(0);
3271 Short_t colmin = fCalibraMode->GetColMin(0);
3272 Short_t colmax = fCalibraMode->GetColMax(0);
3273 Float_t gf = fCurrentCoef[0];
3274 Float_t gfs = fCurrentCoef[1];
3275 Float_t gfE = fCurrentCoefE;
3277 (*fDebugStreamer) << "FillFillCH" <<
3278 "detector=" << detector <<
3279 "caligroup=" << caligroup <<
3280 "rowmin=" << rowmin <<
3281 "rowmax=" << rowmax <<
3282 "colmin=" << colmin <<
3283 "colmax=" << colmax <<
3291 for (Int_t k = 0; k < 2304; k++) {
3292 fCurrentCoefDetector[k] = 0.0;
3296 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3300 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
3302 // Calcul the coef from the database choosen
3303 CalculChargeCoefMean(kFALSE);
3305 //stack 2, not stack 2
3307 if(GetStack(fCountDet) == 2) factor = 12;
3310 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3311 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3312 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3313 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3317 //Put default value negative
3318 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3319 fCurrentCoefE = 0.0;
3328 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3329 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
3332 // For the case where there are not enough entries in the histograms
3333 // of the calibration group, the value present in the choosen database
3334 // will be put. A negativ sign enables to know that a fit was not possible.
3336 if (fDebugLevel == 1) {
3337 AliInfo("The element has not enough statistic to be fitted");
3339 else if (fNbDet > 0) {
3341 Int_t firstdetector = fCountDet;
3342 Int_t lastdetector = fCountDet+fNbDet;
3343 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3344 // loop over detectors
3345 for(Int_t det = firstdetector; det < lastdetector; det++){
3347 //Set the calibration object again
3351 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3353 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3354 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3355 ,(Int_t) GetStack(fCountDet)
3356 ,(Int_t) GetSector(fCountDet),1);
3357 // Reconstruct row min row max
3358 ReconstructFitRowMinRowMax(idect,1);
3360 // Calcul the coef from the database choosen for the detector
3361 CalculVdriftCoefMean();
3364 //stack 2, not stack 2
3366 if(GetStack(fCountDet) == 2) factor = 12;
3369 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3370 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3371 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3372 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3373 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3377 //Put default value negative
3378 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3379 fCurrentCoefE = 0.0;
3380 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3381 fCurrentCoefE2 = 0.0;
3387 if(fDebugLevel > 1){
3389 if ( !fDebugStreamer ) {
3391 TDirectory *backup = gDirectory;
3392 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3393 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3397 Int_t detector = fCountDet;
3398 Int_t caligroup = idect;
3399 Short_t rowmin = fCalibraMode->GetRowMin(1);
3400 Short_t rowmax = fCalibraMode->GetRowMax(1);
3401 Short_t colmin = fCalibraMode->GetColMin(1);
3402 Short_t colmax = fCalibraMode->GetColMax(1);
3403 Float_t vf = fCurrentCoef[0];
3404 Float_t vs = fCurrentCoef[1];
3405 Float_t vfE = fCurrentCoefE;
3406 Float_t t0f = fCurrentCoef2[0];
3407 Float_t t0s = fCurrentCoef2[1];
3408 Float_t t0E = fCurrentCoefE2;
3412 (* fDebugStreamer) << "FillFillPH"<<
3413 "detector="<<detector<<
3414 "nentries="<<nentries<<
3415 "caligroup="<<caligroup<<
3429 for (Int_t k = 0; k < 2304; k++) {
3430 fCurrentCoefDetector[k] = 0.0;
3431 fCurrentCoefDetector2[k] = 0.0;
3435 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3439 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3441 CalculVdriftCoefMean();
3444 //stack 2 and not stack 2
3446 if(GetStack(fCountDet) == 2) factor = 12;
3450 // Fill the fCurrentCoefDetector 2
3451 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3452 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3453 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3454 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3458 // Put the default value
3459 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3460 fCurrentCoefE = 0.0;
3461 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3462 fCurrentCoefE2 = 0.0;
3464 FillFillPH(idect,nentries);
3473 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3474 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3477 // For the case where there are not enough entries in the histograms
3478 // of the calibration group, the value present in the choosen database
3479 // will be put. A negativ sign enables to know that a fit was not possible.
3482 if (fDebugLevel == 1) {
3483 AliInfo("The element has not enough statistic to be fitted");
3485 else if (fNbDet > 0){
3487 Int_t firstdetector = fCountDet;
3488 Int_t lastdetector = fCountDet+fNbDet;
3489 // AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3491 // loop over detectors
3492 for(Int_t det = firstdetector; det < lastdetector; det++){
3494 //Set the calibration object again
3498 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3500 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3501 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3502 ,(Int_t) GetStack(fCountDet)
3503 ,(Int_t) GetSector(fCountDet),2);
3504 // Reconstruct row min row max
3505 ReconstructFitRowMinRowMax(idect,2);
3507 // Calcul the coef from the database choosen for the detector
3508 CalculPRFCoefMean();
3510 //stack 2, not stack 2
3512 if(GetStack(fCountDet) == 2) factor = 12;
3515 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3516 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3517 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3518 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3522 //Put default value negative
3523 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3524 fCurrentCoefE = 0.0;
3529 if(fDebugLevel > 1){
3531 if ( !fDebugStreamer ) {
3533 TDirectory *backup = gDirectory;
3534 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3535 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3538 Int_t detector = fCountDet;
3539 Int_t layer = GetLayer(fCountDet);
3540 Int_t caligroup = idect;
3541 Short_t rowmin = fCalibraMode->GetRowMin(2);
3542 Short_t rowmax = fCalibraMode->GetRowMax(2);
3543 Short_t colmin = fCalibraMode->GetColMin(2);
3544 Short_t colmax = fCalibraMode->GetColMax(2);
3545 Float_t widf = fCurrentCoef[0];
3546 Float_t wids = fCurrentCoef[1];
3547 Float_t widfE = fCurrentCoefE;
3549 (* fDebugStreamer) << "FillFillPRF"<<
3550 "detector="<<detector<<
3552 "caligroup="<<caligroup<<
3563 for (Int_t k = 0; k < 2304; k++) {
3564 fCurrentCoefDetector[k] = 0.0;
3568 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3572 // AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3574 CalculPRFCoefMean();
3576 // stack 2 and not stack 2
3578 if(GetStack(fCountDet) == 2) factor = 12;
3582 // Fill the fCurrentCoefDetector
3583 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3584 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3585 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3589 // Put the default value
3590 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3591 fCurrentCoefE = 0.0;
3599 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3600 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3603 // For the case where there are not enough entries in the histograms
3604 // of the calibration group, the value present in the choosen database
3605 // will be put. A negativ sign enables to know that a fit was not possible.
3608 // Calcul the coef from the database choosen
3609 CalculVdriftLorentzCoef();
3612 if(GetStack(fCountDet) == 2) factor = 1728;
3616 // Fill the fCurrentCoefDetector
3617 for (Int_t k = 0; k < factor; k++) {
3618 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3619 // should be negative
3620 fCurrentCoefDetector2[k] = fCurrentCoef2[1]+100.0;
3624 //Put default opposite sign only for vdrift
3625 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3626 fCurrentCoefE = 0.0;
3627 fCurrentCoef2[0] = fCurrentCoef2[1]+100.0;
3628 fCurrentCoefE2 = 0.0;
3630 FillFillLinearFitter();
3635 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3636 Bool_t AliTRDCalibraFit::NotEnoughStatisticExbAlt()
3639 // For the case where there are not enough entries in the histograms
3640 // of the calibration group, the value present in the choosen database
3641 // will be put. A negativ sign enables to know that a fit was not possible.
3645 if(GetStack(fCountDet) == 2) factor = 1728;
3649 // Fill the fCurrentCoefDetector
3650 for (Int_t k = 0; k < factor; k++) {
3651 fCurrentCoefDetector2[k] = 100.0;
3654 fCurrentCoef2[0] = 100.0;
3655 fCurrentCoefE2 = 0.0;
3662 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3663 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3666 // Fill the coefficients found with the fits or other
3667 // methods from the Fit functions
3670 if (fDebugLevel != 1) {
3672 Int_t firstdetector = fCountDet;
3673 Int_t lastdetector = fCountDet+fNbDet;
3674 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3675 // loop over detectors
3676 for(Int_t det = firstdetector; det < lastdetector; det++){
3678 //Set the calibration object again
3682 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3684 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3685 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3686 ,(Int_t) GetStack(fCountDet)
3687 ,(Int_t) GetSector(fCountDet),0);
3688 // Reconstruct row min row max
3689 ReconstructFitRowMinRowMax(idect,0);
3691 // Calcul the coef from the database choosen for the detector
3692 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3693 else CalculChargeCoefMean(kTRUE);
3695 //stack 2, not stack 2
3697 if(GetStack(fCountDet) == 2) factor = 12;
3700 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3701 Double_t coeftoput = 1.0;
3702 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3703 else coeftoput = fCurrentCoef[0];
3704 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3705 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3706 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3713 if(fDebugLevel > 1){
3715 if ( !fDebugStreamer ) {
3717 TDirectory *backup = gDirectory;
3718 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3719 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3722 Int_t detector = fCountDet;
3723 Int_t caligroup = idect;
3724 Short_t rowmin = fCalibraMode->GetRowMin(0);
3725 Short_t rowmax = fCalibraMode->GetRowMax(0);
3726 Short_t colmin = fCalibraMode->GetColMin(0);
3727 Short_t colmax = fCalibraMode->GetColMax(0);
3728 Float_t gf = fCurrentCoef[0];
3729 Float_t gfs = fCurrentCoef[1];
3730 Float_t gfE = fCurrentCoefE;
3732 (*fDebugStreamer) << "FillFillCH" <<
3733 "detector=" << detector <<
3734 "caligroup=" << caligroup <<
3735 "rowmin=" << rowmin <<
3736 "rowmax=" << rowmax <<
3737 "colmin=" << colmin <<
3738 "colmax=" << colmax <<
3746 for (Int_t k = 0; k < 2304; k++) {
3747 fCurrentCoefDetector[k] = 0.0;
3751 //printf("Check the count now: fCountDet %d\n",fCountDet);
3756 if(GetStack(fCountDet) == 2) factor = 12;
3759 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3760 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3761 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3772 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3773 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3776 // Fill the coefficients found with the fits or other
3777 // methods from the Fit functions
3780 if (fDebugLevel != 1) {
3783 Int_t firstdetector = fCountDet;
3784 Int_t lastdetector = fCountDet+fNbDet;
3785 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3787 // loop over detectors
3788 for(Int_t det = firstdetector; det < lastdetector; det++){
3790 //Set the calibration object again
3794 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3796 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3797 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3798 ,(Int_t) GetStack(fCountDet)
3799 ,(Int_t) GetSector(fCountDet),1);
3800 // Reconstruct row min row max
3801 ReconstructFitRowMinRowMax(idect,1);
3803 // Calcul the coef from the database choosen for the detector
3804 CalculVdriftCoefMean();
3807 //stack 2, not stack 2
3809 if(GetStack(fCountDet) == 2) factor = 12;
3812 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3813 Double_t coeftoput = 1.5;
3814 Double_t coeftoput2 = 0.0;
3816 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3817 else coeftoput = fCurrentCoef[0];
3819 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3820 else coeftoput2 = fCurrentCoef2[0];
3822 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3823 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3824 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3825 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3833 if(fDebugLevel > 1){
3835 if ( !fDebugStreamer ) {
3837 TDirectory *backup = gDirectory;
3838 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3839 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3843 Int_t detector = fCountDet;
3844 Int_t caligroup = idect;
3845 Short_t rowmin = fCalibraMode->GetRowMin(1);
3846 Short_t rowmax = fCalibraMode->GetRowMax(1);
3847 Short_t colmin = fCalibraMode->GetColMin(1);
3848 Short_t colmax = fCalibraMode->GetColMax(1);
3849 Float_t vf = fCurrentCoef[0];
3850 Float_t vs = fCurrentCoef[1];
3851 Float_t vfE = fCurrentCoefE;
3852 Float_t t0f = fCurrentCoef2[0];
3853 Float_t t0s = fCurrentCoef2[1];
3854 Float_t t0E = fCurrentCoefE2;
3858 (* fDebugStreamer) << "FillFillPH"<<
3859 "detector="<<detector<<
3860 "nentries="<<nentries<<
3861 "caligroup="<<caligroup<<
3875 for (Int_t k = 0; k < 2304; k++) {
3876 fCurrentCoefDetector[k] = 0.0;
3877 fCurrentCoefDetector2[k] = 0.0;
3881 //printf("Check the count now: fCountDet %d\n",fCountDet);
3886 if(GetStack(fCountDet) == 2) factor = 12;
3889 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3890 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3891 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3892 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3896 FillFillPH(idect,nentries);
3901 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3902 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3905 // Fill the coefficients found with the fits or other
3906 // methods from the Fit functions
3909 if (fDebugLevel != 1) {
3912 Int_t firstdetector = fCountDet;
3913 Int_t lastdetector = fCountDet+fNbDet;
3914 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3916 // loop over detectors
3917 for(Int_t det = firstdetector; det < lastdetector; det++){
3919 //Set the calibration object again
3923 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3925 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3926 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3927 ,(Int_t) GetStack(fCountDet)
3928 ,(Int_t) GetSector(fCountDet),2);
3929 // Reconstruct row min row max
3930 ReconstructFitRowMinRowMax(idect,2);
3932 // Calcul the coef from the database choosen for the detector
3933 CalculPRFCoefMean();
3935 //stack 2, not stack 2
3937 if(GetStack(fCountDet) == 2) factor = 12;
3940 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3941 Double_t coeftoput = 1.0;
3942 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3943 else coeftoput = fCurrentCoef[0];
3944 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3945 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3946 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3953 if(fDebugLevel > 1){
3955 if ( !fDebugStreamer ) {
3957 TDirectory *backup = gDirectory;
3958 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3959 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3962 Int_t detector = fCountDet;
3963 Int_t layer = GetLayer(fCountDet);
3964 Int_t caligroup = idect;
3965 Short_t rowmin = fCalibraMode->GetRowMin(2);
3966 Short_t rowmax = fCalibraMode->GetRowMax(2);
3967 Short_t colmin = fCalibraMode->GetColMin(2);
3968 Short_t colmax = fCalibraMode->GetColMax(2);
3969 Float_t widf = fCurrentCoef[0];
3970 Float_t wids = fCurrentCoef[1];
3971 Float_t widfE = fCurrentCoefE;
3973 (* fDebugStreamer) << "FillFillPRF"<<
3974 "detector="<<detector<<
3976 "caligroup="<<caligroup<<
3987 for (Int_t k = 0; k < 2304; k++) {
3988 fCurrentCoefDetector[k] = 0.0;
3992 //printf("Check the count now: fCountDet %d\n",fCountDet);
3997 if(GetStack(fCountDet) == 2) factor = 12;
4000 // Pointer to the branch
4001 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
4002 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
4003 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
4013 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4014 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
4017 // Fill the coefficients found with the fits or other
4018 // methods from the Fit functions
4022 if(GetStack(fCountDet) == 2) factor = 1728;
4025 // Pointer to the branch
4026 for (Int_t k = 0; k < factor; k++) {
4027 fCurrentCoefDetector[k] = fCurrentCoef[0];
4028 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4031 FillFillLinearFitter();
4036 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4037 Bool_t AliTRDCalibraFit::FillInfosFitExbAlt()
4040 // Fill the coefficients found with the fits or other
4041 // methods from the Fit functions
4045 if(GetStack(fCountDet) == 2) factor = 1728;
4048 // Pointer to the branch
4049 for (Int_t k = 0; k < factor; k++) {
4050 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4058 //________________________________________________________________________________
4059 void AliTRDCalibraFit::FillFillCH(Int_t idect)
4062 // DebugStream and fVectorFit
4065 // End of one detector
4066 if (idect == (fCount-1)) {
4069 for (Int_t k = 0; k < 2304; k++) {
4070 fCurrentCoefDetector[k] = 0.0;
4074 if(fDebugLevel > 1){
4076 if ( !fDebugStreamer ) {
4078 TDirectory *backup = gDirectory;
4079 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
4080 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4083 Int_t detector = fCountDet;
4084 Int_t caligroup = idect;
4085 Short_t rowmin = fCalibraMode->GetRowMin(0);
4086 Short_t rowmax = fCalibraMode->GetRowMax(0);
4087 Short_t colmin = fCalibraMode->GetColMin(0);
4088 Short_t colmax = fCalibraMode->GetColMax(0);
4089 Float_t gf = fCurrentCoef[0];
4090 Float_t gfs = fCurrentCoef[1];
4091 Float_t gfE = fCurrentCoefE;
4093 (*fDebugStreamer) << "FillFillCH" <<
4094 "detector=" << detector <<
4095 "caligroup=" << caligroup <<
4096 "rowmin=" << rowmin <<
4097 "rowmax=" << rowmax <<
4098 "colmin=" << colmin <<
4099 "colmax=" << colmax <<
4107 //________________________________________________________________________________
4108 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
4111 // DebugStream and fVectorFit and fVectorFit2
4114 // End of one detector
4115 if (idect == (fCount-1)) {
4119 for (Int_t k = 0; k < 2304; k++) {
4120 fCurrentCoefDetector[k] = 0.0;
4121 fCurrentCoefDetector2[k] = 0.0;
4125 if(fDebugLevel > 1){
4127 if ( !fDebugStreamer ) {
4129 TDirectory *backup = gDirectory;
4130 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
4131 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4135 Int_t detector = fCountDet;
4136 Int_t caligroup = idect;
4137 Short_t rowmin = fCalibraMode->GetRowMin(1);
4138 Short_t rowmax = fCalibraMode->GetRowMax(1);
4139 Short_t colmin = fCalibraMode->GetColMin(1);
4140 Short_t colmax = fCalibraMode->GetColMax(1);
4141 Float_t vf = fCurrentCoef[0];
4142 Float_t vs = fCurrentCoef[1];
4143 Float_t vfE = fCurrentCoefE;
4144 Float_t t0f = fCurrentCoef2[0];
4145 Float_t t0s = fCurrentCoef2[1];
4146 Float_t t0E = fCurrentCoefE2;
4150 (* fDebugStreamer) << "FillFillPH"<<
4151 "detector="<<detector<<
4152 "nentries="<<nentries<<
4153 "caligroup="<<caligroup<<
4168 //________________________________________________________________________________
4169 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
4172 // DebugStream and fVectorFit
4175 // End of one detector
4176 if (idect == (fCount-1)) {
4179 for (Int_t k = 0; k < 2304; k++) {
4180 fCurrentCoefDetector[k] = 0.0;
4185 if(fDebugLevel > 1){
4187 if ( !fDebugStreamer ) {
4189 TDirectory *backup = gDirectory;
4190 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4191 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4194 Int_t detector = fCountDet;
4195 Int_t layer = GetLayer(fCountDet);
4196 Int_t caligroup = idect;
4197 Short_t rowmin = fCalibraMode->GetRowMin(2);
4198 Short_t rowmax = fCalibraMode->GetRowMax(2);
4199 Short_t colmin = fCalibraMode->GetColMin(2);
4200 Short_t colmax = fCalibraMode->GetColMax(2);
4201 Float_t widf = fCurrentCoef[0];
4202 Float_t wids = fCurrentCoef[1];
4203 Float_t widfE = fCurrentCoefE;
4205 (* fDebugStreamer) << "FillFillPRF"<<
4206 "detector="<<detector<<
4208 "caligroup="<<caligroup<<
4220 //________________________________________________________________________________
4221 void AliTRDCalibraFit::FillFillLinearFitter()
4224 // DebugStream and fVectorFit
4227 // End of one detector
4233 for (Int_t k = 0; k < 2304; k++) {
4234 fCurrentCoefDetector[k] = 0.0;
4235 fCurrentCoefDetector2[k] = 0.0;
4239 if(fDebugLevel > 1){
4241 if ( !fDebugStreamer ) {
4243 TDirectory *backup = gDirectory;
4244 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
4245 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4248 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4249 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4250 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4251 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4252 Float_t tiltangle = padplane->GetTiltingAngle();
4253 Int_t detector = fCountDet;
4254 Int_t stack = GetStack(fCountDet);
4255 Int_t layer = GetLayer(fCountDet);
4256 Float_t vf = fCurrentCoef[0];
4257 Float_t vs = fCurrentCoef[1];
4258 Float_t vfE = fCurrentCoefE;
4259 Float_t lorentzangler = fCurrentCoef2[0];
4260 Float_t elorentzangler = fCurrentCoefE2;
4261 Float_t lorentzangles = fCurrentCoef2[1];
4263 (* fDebugStreamer) << "FillFillLinearFitter"<<
4264 "detector="<<detector<<
4269 "tiltangle="<<tiltangle<<
4273 "lorentzangler="<<lorentzangler<<
4274 "Elorentzangler="<<elorentzangler<<
4275 "lorentzangles="<<lorentzangles<<
4280 //________________________________________________________________________________
4281 void AliTRDCalibraFit::FillFillExbAlt()
4284 // DebugStream and fVectorFit
4287 // End of one detector
4292 for (Int_t k = 0; k < 2304; k++) {
4293 fCurrentCoefDetector2[k] = 0.0;
4297 if(fDebugLevel > 1){
4299 if ( !fDebugStreamer ) {
4301 TDirectory *backup = gDirectory;
4302 fDebugStreamer = new TTreeSRedirector("TRDDebugFitExbAlt.root");
4303 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4306 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4307 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4308 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4309 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4310 Float_t tiltangle = padplane->GetTiltingAngle();
4311 Int_t detector = fCountDet;
4312 Int_t stack = GetStack(fCountDet);
4313 Int_t layer = GetLayer(fCountDet);
4314 Float_t vf = fCurrentCoef2[0];
4315 Float_t vfE = fCurrentCoefE2;
4317 (* fDebugStreamer) << "FillFillLinearFitter"<<
4318 "detector="<<detector<<
4323 "tiltangle="<<tiltangle<<
4331 //____________Calcul Coef Mean_________________________________________________
4333 //_____________________________________________________________________________
4334 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
4337 // For the detector Dect calcul the mean time 0
4338 // for the calibration group idect from the choosen database
4341 fCurrentCoef2[1] = 0.0;
4342 if(fDebugLevel != 1){
4343 if(((fCalibraMode->GetNz(1) > 0) ||
4344 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4346 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4347 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4348 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4352 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4358 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4362 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
4363 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
4364 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4367 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
4375 //_____________________________________________________________________________
4376 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
4379 // For the detector Dect calcul the mean gain factor
4380 // for the calibration group idect from the choosen database
4383 fCurrentCoef[1] = 0.0;
4384 if(fDebugLevel != 1){
4385 if (((fCalibraMode->GetNz(0) > 0) ||
4386 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
4387 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4388 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4389 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4390 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4393 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4397 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4398 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
4403 //_____________________________________________________________________________
4404 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
4407 // For the detector Dect calcul the mean sigma of pad response
4408 // function for the calibration group idect from the choosen database
4411 fCurrentCoef[1] = 0.0;
4412 if(fDebugLevel != 1){
4413 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4414 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4415 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
4418 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
4422 //_____________________________________________________________________________
4423 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
4426 // For the detector dect calcul the mean drift velocity for the
4427 // calibration group idect from the choosen database
4430 fCurrentCoef[1] = 0.0;
4431 if(fDebugLevel != 1){
4432 if (((fCalibraMode->GetNz(1) > 0) ||
4433 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4435 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4436 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4437 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4441 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4446 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4451 //_____________________________________________________________________________
4452 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4455 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4458 fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet);
4459 fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet);
4463 //_____________________________________________________________________________
4464 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4467 // Default width of the PRF if there is no database as reference
4472 //case 0: return 0.515;
4473 //case 1: return 0.502;
4474 //case 2: return 0.491;
4475 //case 3: return 0.481;
4476 //case 4: return 0.471;
4477 //case 5: return 0.463;
4478 //default: return 0.0;
4481 case 0: return 0.538429;
4482 case 1: return 0.524302;
4483 case 2: return 0.511591;
4484 case 3: return 0.500140;
4485 case 4: return 0.489821;
4486 case 5: return 0.480524;
4487 default: return 0.0;
4490 //________________________________________________________________________________
4491 void AliTRDCalibraFit::SetCalROC(Int_t i)
4494 // Set the calib object for fCountDet
4497 Float_t value = 0.0;
4499 //Get the CalDet object
4501 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4503 AliInfo("Could not get calibDB");
4510 fCalROC->~AliTRDCalROC();
4511 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4512 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4516 fCalROC->~AliTRDCalROC();
4517 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4518 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4520 fCalROC2->~AliTRDCalROC();
4521 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4522 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4526 fCalROC->~AliTRDCalROC();
4527 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4528 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4537 if(fCalROC) delete fCalROC;
4538 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4539 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4540 fCalROC->SetValue(k,1.0);
4544 if(fCalROC) delete fCalROC;
4545 if(fCalROC2) delete fCalROC2;
4546 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4547 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4548 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4549 fCalROC->SetValue(k,1.0);
4550 fCalROC2->SetValue(k,0.0);
4554 if(fCalROC) delete fCalROC;
4555 value = GetPRFDefault(GetLayer(fCountDet));
4556 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4557 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4558 fCalROC->SetValue(k,value);
4566 //____________Fit Methods______________________________________________________
4568 //_____________________________________________________________________________
4569 void AliTRDCalibraFit::FitPente(TH1* projPH)
4572 // Slope methode for the drift velocity
4576 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4583 fCurrentCoefE = 0.0;
4584 fCurrentCoefE2 = 0.0;
4585 fCurrentCoef[0] = 0.0;
4586 fCurrentCoef2[0] = 0.0;
4587 TLine *line = new TLine();
4590 TAxis *xpph = projPH->GetXaxis();
4591 Int_t nbins = xpph->GetNbins();
4592 Double_t lowedge = xpph->GetBinLowEdge(1);
4593 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4594 Double_t widbins = (upedge - lowedge) / nbins;
4595 Double_t limit = upedge + 0.5 * widbins;
4598 // Beginning of the signal
4599 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4600 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4601 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4603 binmax = (Int_t) pentea->GetMaximumBin();
4606 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4608 if (binmax >= nbins) {
4611 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4613 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4614 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4615 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4616 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4617 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4618 if (TMath::Abs(l3P2am) > 0.00000001) {
4619 fPhd[0] = -(l3P1am / (2 * l3P2am));
4622 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4623 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4626 // Amplification region
4629 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4630 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4637 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4639 if (binmax >= nbins) {
4642 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4644 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4645 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4646 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4647 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4648 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4649 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4650 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4652 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4653 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4656 fCurrentCoefE2 = fCurrentCoefE;
4659 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4660 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4661 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4664 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4667 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4669 if (binmin >= nbins) {
4672 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4677 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4678 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4679 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4680 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4681 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4682 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4683 if (TMath::Abs(l3P2dr) > 0.00000001) {
4684 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4686 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4687 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4689 Float_t fPhdt0 = 0.0;
4690 Float_t t0Shift = 0.0;
4693 t0Shift = fT0Shift1;
4697 t0Shift = fT0Shift0;
4700 if ((fPhd[2] > fPhd[0]) &&
4701 (fPhd[2] > fPhd[1]) &&
4702 (fPhd[1] > fPhd[0]) &&
4704 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4705 fNumberFitSuccess++;
4707 if (fPhdt0 >= 0.0) {
4708 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4709 if (fCurrentCoef2[0] < -3.0) {
4710 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4714 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4719 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4720 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4723 if (fDebugLevel == 1) {
4724 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4727 line->SetLineColor(2);
4728 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4729 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4730 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4731 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4732 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4733 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4734 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4735 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4738 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4747 //_____________________________________________________________________________
4748 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4751 // Slope methode but with polynomes de Lagrange
4755 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4758 //Double_t *x = new Double_t[5];
4759 //Double_t *y = new Double_t[5];
4776 fCurrentCoefE = 0.0;
4777 fCurrentCoefE2 = 1.0;
4778 fCurrentCoef[0] = 0.0;
4779 fCurrentCoef2[0] = 0.0;
4780 TLine *line = new TLine();
4781 TF1 * polynome = 0x0;
4782 TF1 * polynomea = 0x0;
4783 TF1 * polynomeb = 0x0;
4791 TAxis *xpph = projPH->GetXaxis();
4792 Int_t nbins = xpph->GetNbins();
4793 Double_t lowedge = xpph->GetBinLowEdge(1);
4794 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4795 Double_t widbins = (upedge - lowedge) / nbins;
4796 Double_t limit = upedge + 0.5 * widbins;
4801 // Beginning of the signal
4802 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4803 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4804 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4807 binmax = (Int_t) pentea->GetMaximumBin();
4809 Double_t minnn = 0.0;
4810 Double_t maxxx = 0.0;
4812 Int_t kase = nbins-binmax;
4820 minnn = pentea->GetBinCenter(binmax-2);
4821 maxxx = pentea->GetBinCenter(binmax);
4822 x[0] = pentea->GetBinCenter(binmax-2);
4823 x[1] = pentea->GetBinCenter(binmax-1);
4824 x[2] = pentea->GetBinCenter(binmax);
4825 y[0] = pentea->GetBinContent(binmax-2);
4826 y[1] = pentea->GetBinContent(binmax-1);
4827 y[2] = pentea->GetBinContent(binmax);
4828 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4829 //AliInfo("At the limit for beginning!");
4832 minnn = pentea->GetBinCenter(binmax-2);
4833 maxxx = pentea->GetBinCenter(binmax+1);
4834 x[0] = pentea->GetBinCenter(binmax-2);
4835 x[1] = pentea->GetBinCenter(binmax-1);
4836 x[2] = pentea->GetBinCenter(binmax);
4837 x[3] = pentea->GetBinCenter(binmax+1);
4838 y[0] = pentea->GetBinContent(binmax-2);
4839 y[1] = pentea->GetBinContent(binmax-1);
4840 y[2] = pentea->GetBinContent(binmax);
4841 y[3] = pentea->GetBinContent(binmax+1);
4842 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4850 minnn = pentea->GetBinCenter(binmax);
4851 maxxx = pentea->GetBinCenter(binmax+2);
4852 x[0] = pentea->GetBinCenter(binmax);
4853 x[1] = pentea->GetBinCenter(binmax+1);
4854 x[2] = pentea->GetBinCenter(binmax+2);
4855 y[0] = pentea->GetBinContent(binmax);
4856 y[1] = pentea->GetBinContent(binmax+1);
4857 y[2] = pentea->GetBinContent(binmax+2);
4858 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4861 minnn = pentea->GetBinCenter(binmax-1);
4862 maxxx = pentea->GetBinCenter(binmax+2);
4863 x[0] = pentea->GetBinCenter(binmax-1);
4864 x[1] = pentea->GetBinCenter(binmax);
4865 x[2] = pentea->GetBinCenter(binmax+1);
4866 x[3] = pentea->GetBinCenter(binmax+2);
4867 y[0] = pentea->GetBinContent(binmax-1);
4868 y[1] = pentea->GetBinContent(binmax);
4869 y[2] = pentea->GetBinContent(binmax+1);
4870 y[3] = pentea->GetBinContent(binmax+2);
4871 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4874 minnn = pentea->GetBinCenter(binmax-2);
4875 maxxx = pentea->GetBinCenter(binmax+2);
4876 x[0] = pentea->GetBinCenter(binmax-2);
4877 x[1] = pentea->GetBinCenter(binmax-1);
4878 x[2] = pentea->GetBinCenter(binmax);
4879 x[3] = pentea->GetBinCenter(binmax+1);
4880 x[4] = pentea->GetBinCenter(binmax+2);
4881 y[0] = pentea->GetBinContent(binmax-2);
4882 y[1] = pentea->GetBinContent(binmax-1);
4883 y[2] = pentea->GetBinContent(binmax);
4884 y[3] = pentea->GetBinContent(binmax+1);
4885 y[4] = pentea->GetBinContent(binmax+2);
4886 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4894 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4895 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4897 Double_t step = (maxxx-minnn)/10000;
4899 Double_t maxvalue = 0.0;
4900 Double_t placemaximum = minnn;
4901 for(Int_t o = 0; o < 10000; o++){
4902 if(o == 0) maxvalue = polynomeb->Eval(l);
4903 if(maxvalue < (polynomeb->Eval(l))){
4904 maxvalue = polynomeb->Eval(l);
4909 fPhd[0] = placemaximum;
4912 // Amplification region
4915 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4916 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4922 Double_t minn = 0.0;
4923 Double_t maxx = 0.0;
4935 Int_t kase1 = nbins - binmax;
4937 //Determination of minn and maxx
4938 //case binmax = nbins
4943 minn = projPH->GetBinCenter(binmax-2);
4944 maxx = projPH->GetBinCenter(binmax);
4945 x[0] = projPH->GetBinCenter(binmax-2);
4946 x[1] = projPH->GetBinCenter(binmax-1);
4947 x[2] = projPH->GetBinCenter(binmax);
4948 y[0] = projPH->GetBinContent(binmax-2);
4949 y[1] = projPH->GetBinContent(binmax-1);
4950 y[2] = projPH->GetBinContent(binmax);
4951 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4952 //AliInfo("At the limit for the drift!");
4955 minn = projPH->GetBinCenter(binmax-2);
4956 maxx = projPH->GetBinCenter(binmax+1);
4957 x[0] = projPH->GetBinCenter(binmax-2);
4958 x[1] = projPH->GetBinCenter(binmax-1);
4959 x[2] = projPH->GetBinCenter(binmax);
4960 x[3] = projPH->GetBinCenter(binmax+1);
4961 y[0] = projPH->GetBinContent(binmax-2);
4962 y[1] = projPH->GetBinContent(binmax-1);
4963 y[2] = projPH->GetBinContent(binmax);
4964 y[3] = projPH->GetBinContent(binmax+1);
4965 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4974 minn = projPH->GetBinCenter(binmax);
4975 maxx = projPH->GetBinCenter(binmax+2);
4976 x[0] = projPH->GetBinCenter(binmax);
4977 x[1] = projPH->GetBinCenter(binmax+1);
4978 x[2] = projPH->GetBinCenter(binmax+2);
4979 y[0] = projPH->GetBinContent(binmax);
4980 y[1] = projPH->GetBinContent(binmax+1);
4981 y[2] = projPH->GetBinContent(binmax+2);
4982 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4985 minn = projPH->GetBinCenter(binmax-1);
4986 maxx = projPH->GetBinCenter(binmax+2);
4987 x[0] = projPH->GetBinCenter(binmax-1);
4988 x[1] = projPH->GetBinCenter(binmax);
4989 x[2] = projPH->GetBinCenter(binmax+1);
4990 x[3] = projPH->GetBinCenter(binmax+2);
4991 y[0] = projPH->GetBinContent(binmax-1);
4992 y[1] = projPH->GetBinContent(binmax);
4993 y[2] = projPH->GetBinContent(binmax+1);
4994 y[3] = projPH->GetBinContent(binmax+2);
4995 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4998 minn = projPH->GetBinCenter(binmax-2);
4999 maxx = projPH->GetBinCenter(binmax+2);
5000 x[0] = projPH->GetBinCenter(binmax-2);
5001 x[1] = projPH->GetBinCenter(binmax-1);
5002 x[2] = projPH->GetBinCenter(binmax);
5003 x[3] = projPH->GetBinCenter(binmax+1);
5004 x[4] = projPH->GetBinCenter(binmax+2);
5005 y[0] = projPH->GetBinContent(binmax-2);
5006 y[1] = projPH->GetBinContent(binmax-1);
5007 y[2] = projPH->GetBinContent(binmax);
5008 y[3] = projPH->GetBinContent(binmax+1);
5009 y[4] = projPH->GetBinContent(binmax+2);
5010 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5017 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
5018 polynomea->SetParameters(c0,c1,c2,c3,c4);
5020 Double_t step = (maxx-minn)/1000;
5022 Double_t maxvalue = 0.0;
5023 Double_t placemaximum = minn;
5024 for(Int_t o = 0; o < 1000; o++){
5025 if(o == 0) maxvalue = polynomea->Eval(l);
5026 if(maxvalue < (polynomea->Eval(l))){
5027 maxvalue = polynomea->Eval(l);
5032 fPhd[1] = placemaximum;
5036 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
5037 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5038 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5041 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5047 //AliInfo("Put the binmax from 1 to 2 to enable the fit");
5051 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
5052 //AliInfo("Too many fluctuations at the end!");
5055 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
5056 //AliInfo("Too many fluctuations at the end!");
5059 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
5060 //AliInfo("No entries for the next bin!");
5061 pente->SetBinContent(binmin,0);
5062 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5078 Bool_t case1 = kFALSE;
5079 Bool_t case2 = kFALSE;
5080 Bool_t case4 = kFALSE;
5082 //Determination of min and max
5083 //case binmin <= nbins-3
5085 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5086 min = pente->GetBinCenter(binmin-2);
5087 max = pente->GetBinCenter(binmin+2);
5088 x[0] = pente->GetBinCenter(binmin-2);
5089 x[1] = pente->GetBinCenter(binmin-1);
5090 x[2] = pente->GetBinCenter(binmin);
5091 x[3] = pente->GetBinCenter(binmin+1);
5092 x[4] = pente->GetBinCenter(binmin+2);
5093 y[0] = pente->GetBinContent(binmin-2);
5094 y[1] = pente->GetBinContent(binmin-1);
5095 y[2] = pente->GetBinContent(binmin);
5096 y[3] = pente->GetBinContent(binmin+1);
5097 y[4] = pente->GetBinContent(binmin+2);
5098 //Calcul the polynome de Lagrange
5099 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5101 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5102 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5103 //AliInfo("polynome 4 false 1");
5106 if(((binmin+3) <= (nbins-1)) &&
5107 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
5108 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
5109 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
5110 //AliInfo("polynome 4 false 2");
5114 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5115 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
5116 //AliInfo("polynome 4 case 1");
5119 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
5120 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5121 //AliInfo("polynome 4 case 4");
5126 //case binmin = nbins-2
5128 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5130 min = pente->GetBinCenter(binmin-2);
5131 max = pente->GetBinCenter(binmin+1);
5132 x[0] = pente->GetBinCenter(binmin-2);
5133 x[1] = pente->GetBinCenter(binmin-1);
5134 x[2] = pente->GetBinCenter(binmin);
5135 x[3] = pente->GetBinCenter(binmin+1);
5136 y[0] = pente->GetBinContent(binmin-2);
5137 y[1] = pente->GetBinContent(binmin-1);
5138 y[2] = pente->GetBinContent(binmin);
5139 y[3] = pente->GetBinContent(binmin+1);
5140 //Calcul the polynome de Lagrange
5141 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5142 //richtung +: nothing
5144 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5145 //AliInfo("polynome 3- case 2");
5150 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5152 min = pente->GetBinCenter(binmin-1);
5153 max = pente->GetBinCenter(binmin+2);
5154 x[0] = pente->GetBinCenter(binmin-1);
5155 x[1] = pente->GetBinCenter(binmin);
5156 x[2] = pente->GetBinCenter(binmin+1);
5157 x[3] = pente->GetBinCenter(binmin+2);
5158 y[0] = pente->GetBinContent(binmin-1);
5159 y[1] = pente->GetBinContent(binmin);
5160 y[2] = pente->GetBinContent(binmin+1);
5161 y[3] = pente->GetBinContent(binmin+2);
5162 //Calcul the polynome de Lagrange
5163 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5165 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5166 //AliInfo("polynome 3+ case 2");
5171 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
5172 min = pente->GetBinCenter(binmin);
5173 max = pente->GetBinCenter(binmin+2);
5174 x[0] = pente->GetBinCenter(binmin);
5175 x[1] = pente->GetBinCenter(binmin+1);
5176 x[2] = pente->GetBinCenter(binmin+2);
5177 y[0] = pente->GetBinContent(binmin);
5178 y[1] = pente->GetBinContent(binmin+1);
5179 y[2] = pente->GetBinContent(binmin+2);
5180 //Calcul the polynome de Lagrange
5181 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5183 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5184 //AliInfo("polynome 2+ false");
5189 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5191 min = pente->GetBinCenter(binmin-1);
5192 max = pente->GetBinCenter(binmin+1);
5193 x[0] = pente->GetBinCenter(binmin-1);
5194 x[1] = pente->GetBinCenter(binmin);
5195 x[2] = pente->GetBinCenter(binmin+1);
5196 y[0] = pente->GetBinContent(binmin-1);
5197 y[1] = pente->GetBinContent(binmin);
5198 y[2] = pente->GetBinContent(binmin+1);
5199 //Calcul the polynome de Lagrange
5200 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5201 //richtung +: nothing
5202 //richtung -: nothing
5204 //case binmin = nbins-1
5206 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5207 min = pente->GetBinCenter(binmin-2);
5208 max = pente->GetBinCenter(binmin);
5209 x[0] = pente->GetBinCenter(binmin-2);
5210 x[1] = pente->GetBinCenter(binmin-1);
5211 x[2] = pente->GetBinCenter(binmin);
5212 y[0] = pente->GetBinContent(binmin-2);
5213 y[1] = pente->GetBinContent(binmin-1);
5214 y[2] = pente->GetBinContent(binmin);
5215 //Calcul the polynome de Lagrange
5216 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5217 //AliInfo("At the limit for the drift!");
5218 //fluctuation too big!
5219 //richtung +: nothing
5221 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5222 //AliInfo("polynome 2- false ");
5226 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
5228 //AliInfo("At the limit for the drift and not usable!");
5232 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5234 //AliInfo("For the drift...problem!");
5236 //pass but should not happen
5237 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
5239 //AliInfo("For the drift...problem!");
5243 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5244 polynome->SetParameters(c0,c1,c2,c3,c4);
5245 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5246 Double_t step = (max-min)/1000;
5248 Double_t minvalue = 0.0;
5249 Double_t placeminimum = min;
5250 for(Int_t o = 0; o < 1000; o++){
5251 if(o == 0) minvalue = polynome->Eval(l);
5252 if(minvalue > (polynome->Eval(l))){
5253 minvalue = polynome->Eval(l);
5258 fPhd[2] = placeminimum;
5260 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
5261 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
5262 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
5264 Float_t fPhdt0 = 0.0;
5265 Float_t t0Shift = 0.0;
5268 t0Shift = fT0Shift1;
5272 t0Shift = fT0Shift0;
5275 if ((fPhd[2] > fPhd[0]) &&
5276 (fPhd[2] > fPhd[1]) &&
5277 (fPhd[1] > fPhd[0]) &&
5279 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5280 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5281 else fNumberFitSuccess++;
5282 if (fPhdt0 >= 0.0) {
5283 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5284 //printf("Value of timeoffset %f\n",fCurrentCoef2[0]);
5285 if (fCurrentCoef2[0] < -3.0) {
5286 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5290 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5294 ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
5295 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5297 if((fPhd[1] > fPhd[0]) &&
5299 if (fPhdt0 >= 0.0) {
5300 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5301 if (fCurrentCoef2[0] < -3.0) {
5302 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5306 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5310 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5311 //printf("Fit failed!\n");
5315 if (fDebugLevel == 1) {
5316 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5319 line->SetLineColor(2);
5320 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5321 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5322 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5323 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5324 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5325 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5326 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
5327 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5330 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5337 if(polynome) delete polynome;
5338 if(polynomea) delete polynomea;
5339 if(polynomeb) delete polynomeb;
5340 //if(x) delete [] x;
5341 //if(y) delete [] y;
5342 if(line) delete line;
5347 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5348 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5349 //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]);
5350 projPH->SetDirectory(0);
5354 //_____________________________________________________________________________
5355 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5358 // Fit methode for the drift velocity
5362 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5365 TAxis *xpph = projPH->GetXaxis();
5366 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5368 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5369 fPH->SetParameter(0,0.469); // Scaling
5370 fPH->SetParameter(1,0.18); // Start
5371 fPH->SetParameter(2,0.0857325); // AR
5372 fPH->SetParameter(3,1.89); // DR
5373 fPH->SetParameter(4,0.08); // QA/QD
5374 fPH->SetParameter(5,0.0); // Baseline
5376 TLine *line = new TLine();
5378 fCurrentCoef[0] = 0.0;
5379 fCurrentCoef2[0] = 0.0;
5380 fCurrentCoefE = 0.0;
5381 fCurrentCoefE2 = 0.0;
5383 if (idect%fFitPHPeriode == 0) {
5385 AliInfo(Form("The detector %d will be fitted",idect));
5386 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5387 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5388 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5389 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5390 fPH->SetParameter(4,0.225); // QA/QD
5391 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5393 if (fDebugLevel != 1) {
5394 projPH->Fit(fPH,"0M","",0.0,upedge);
5397 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5399 projPH->Fit(fPH,"M+","",0.0,upedge);
5401 line->SetLineColor(4);
5402 line->DrawLine(fPH->GetParameter(1)
5404 ,fPH->GetParameter(1)
5405 ,projPH->GetMaximum());
5406 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5408 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5409 ,projPH->GetMaximum());
5410 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5412 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5413 ,projPH->GetMaximum());
5416 if (fPH->GetParameter(3) != 0) {
5417 fNumberFitSuccess++;
5418 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
5419 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
5420 fCurrentCoef2[0] = fPH->GetParameter(1);
5421 fCurrentCoefE2 = fPH->GetParError(1);
5424 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5425 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5431 // Put the default value
5432 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5433 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5436 if (fDebugLevel != 1) {
5441 //_____________________________________________________________________________
5442 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5445 // Fit methode for the sigma of the pad response function
5450 fCurrentCoef[0] = 0.0;
5451 fCurrentCoefE = 0.0;
5453 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5455 if(TMath::Abs(ret+4) <= 0.000000001){
5456 fCurrentCoef[0] = -fCurrentCoef[1];
5460 fNumberFitSuccess++;
5461 fCurrentCoef[0] = param[2];
5462 fCurrentCoefE = ret;
5466 //_____________________________________________________________________________
5467 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)
5470 // Fit methode for the sigma of the pad response function
5473 //We should have at least 3 points
5474 if(nBins <=3) return -4.0;
5476 TLinearFitter fitter(3,"pol2");
5477 fitter.StoreData(kFALSE);
5478 fitter.ClearPoints();
5480 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5481 Float_t entries = 0;
5482 Int_t nbbinwithentries = 0;
5483 for (Int_t i=0; i<nBins; i++){
5485 if(arraye[i] > 15) nbbinwithentries++;
5486 //printf("entries for i %d: %f\n",i,arraye[i]);
5488 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5489 //printf("entries %f\n",entries);
5490 //printf("nbbinwithentries %d\n",nbbinwithentries);
5493 Float_t errorm = 0.0;
5494 Float_t errorn = 0.0;
5495 Float_t error = 0.0;
5498 for (Int_t ibin=0;ibin<nBins; ibin++){
5499 Float_t entriesI = arraye[ibin];
5500 Float_t valueI = arraym[ibin];
5501 Double_t xcenter = 0.0;
5503 if ((entriesI>15) && (valueI>0.0)){
5504 xcenter = xMin+(ibin+0.5)*binWidth;
5509 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5510 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5511 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5514 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5515 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5516 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5518 if(TMath::Abs(error) < 0.000000001) continue;
5519 val = TMath::Log(Float_t(valueI));
5520 fitter.AddPoint(&xcenter,val,error);
5524 if(fDebugLevel > 1){
5526 if ( !fDebugStreamer ) {
5528 TDirectory *backup = gDirectory;
5529 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5530 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5533 Int_t detector = fCountDet;
5534 Int_t layer = GetLayer(fCountDet);
5537 (* fDebugStreamer) << "FitGausMIFill"<<
5538 "detector="<<detector<<
5542 "entriesI="<<entriesI<<
5545 "xcenter="<<xcenter<<
5555 if(npoints <=3) return -4.0;
5560 fitter.GetParameters(par);
5561 chi2 = fitter.GetChisquare()/Float_t(npoints);
5564 if (!param) param = new TVectorD(3);
5565 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5566 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5567 Double_t deltax = (fitter.GetParError(2))/x;
5568 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5571 (*param)[1] = par[1]/(-2.*par[2]);
5572 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5573 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5574 if ( lnparam0>307 ) return -4;
5575 (*param)[0] = TMath::Exp(lnparam0);
5577 if(fDebugLevel > 1){
5579 if ( !fDebugStreamer ) {
5581 TDirectory *backup = gDirectory;
5582 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5583 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5586 Int_t detector = fCountDet;
5587 Int_t layer = GetLayer(fCountDet);
5590 (* fDebugStreamer) << "FitGausMIFit"<<
5591 "detector="<<detector<<
5594 "errorsigma="<<chi2<<
5595 "mean="<<(*param)[1]<<
5596 "sigma="<<(*param)[2]<<
5597 "constant="<<(*param)[0]<<
5602 if((chi2/(*param)[2]) > 0.1){
5604 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5609 if(fDebugLevel == 1){
5610 TString name("PRF");
5611 name += (Int_t)xMin;
5612 name += (Int_t)xMax;
5613 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5616 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5617 for(Int_t k = 0; k < nBins; k++){
5618 histo->SetBinContent(k+1,arraym[k]);
5619 histo->SetBinError(k+1,arrayme[k]);
5622 name += "functionf";
5623 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5624 f1->SetParameter(0, (*param)[0]);
5625 f1->SetParameter(1, (*param)[1]);
5626 f1->SetParameter(2, (*param)[2]);
5634 //_____________________________________________________________________________
5635 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5638 // Fit methode for the sigma of the pad response function
5641 fCurrentCoef[0] = 0.0;
5642 fCurrentCoefE = 0.0;
5644 if (fDebugLevel != 1) {
5645 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5648 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5650 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5653 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5654 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5655 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5657 fNumberFitSuccess++;
5660 //_____________________________________________________________________________
5661 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5664 // Fit methode for the sigma of the pad response function
5666 fCurrentCoef[0] = 0.0;
5667 fCurrentCoefE = 0.0;
5668 if (fDebugLevel == 1) {
5669 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5673 fCurrentCoef[0] = projPRF->GetRMS();
5674 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5676 fNumberFitSuccess++;
5679 //_____________________________________________________________________________
5680 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5683 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5686 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5689 Int_t nbins = (Int_t)(nybins/(2*nbg));
5690 Float_t lowedge = -3.0*nbg;
5691 Float_t upedge = lowedge + 3.0;
5694 Double_t xvalues = -0.2*nbg+0.1;
5696 Int_t total = 2*nbg;
5699 for(Int_t k = 0; k < total; k++){
5700 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5702 y = fCurrentCoef[0]*fCurrentCoef[0];
5703 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5706 if(fDebugLevel > 1){
5708 if ( !fDebugStreamer ) {
5710 TDirectory *backup = gDirectory;
5711 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5712 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5715 Int_t detector = fCountDet;
5716 Int_t layer = GetLayer(fCountDet);
5717 Int_t nbtotal = total;
5719 Float_t low = lowedge;
5720 Float_t up = upedge;
5721 Float_t tnp = xvalues;
5722 Float_t wid = fCurrentCoef[0];
5723 Float_t widfE = fCurrentCoefE;
5725 (* fDebugStreamer) << "FitTnpRange0"<<
5726 "detector="<<detector<<
5728 "nbtotal="<<nbtotal<<
5746 fCurrentCoefE = 0.0;
5747 fCurrentCoef[0] = 0.0;
5749 //printf("npoints\n",npoints);
5752 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5757 linearfitter.Eval();
5758 linearfitter.GetParameters(pars0);
5759 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5760 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5761 Double_t min0 = 0.0;
5762 Double_t ermin0 = 0.0;
5763 //Double_t prfe0 = 0.0;
5764 Double_t prf0 = 0.0;
5765 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5766 min0 = -pars0[1]/(2*pars0[2]);
5767 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5768 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5771 prfe0 = linearfitter->GetParError(0)*pointError0
5772 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5773 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5774 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5775 fCurrentCoefE = (Float_t) prfe0;
5777 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5780 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5784 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5787 if(fDebugLevel > 1){
5789 if ( !fDebugStreamer ) {
5791 TDirectory *backup = gDirectory;
5792 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5793 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5796 Int_t detector = fCountDet;
5797 Int_t layer = GetLayer(fCountDet);
5798 Int_t nbtotal = total;
5799 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5800 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5802 (* fDebugStreamer) << "FitTnpRange1"<<
5803 "detector="<<detector<<
5805 "nbtotal="<<nbtotal<<
5809 "npoints="<<npoints<<
5812 "sigmaprf="<<fCurrentCoef[0]<<
5813 "sigprf="<<fCurrentCoef[1]<<
5820 //_____________________________________________________________________________
5821 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5824 // Only mean methode for the gain factor
5827 fCurrentCoef[0] = mean;
5828 fCurrentCoefE = 0.0;
5829 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5830 if (fDebugLevel == 1) {
5831 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5835 CalculChargeCoefMean(kTRUE);
5836 fNumberFitSuccess++;
5838 //_____________________________________________________________________________
5839 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5842 // mean w methode for the gain factor
5846 Int_t nybins = projch->GetNbinsX();
5848 //The weight function
5849 Double_t a = 0.00228515;
5850 Double_t b = -0.00231487;
5851 Double_t c = 0.00044298;
5852 Double_t d = -0.00379239;
5853 Double_t e = 0.00338349;
5863 //A arbitrary error for the moment
5864 fCurrentCoefE = 0.0;
5865 fCurrentCoef[0] = 0.0;
5868 Double_t sumw = 0.0;
5870 Float_t sumAll = (Float_t) nentries;
5871 Int_t sumCurrent = 0;
5872 for(Int_t k = 0; k <nybins; k++){
5873 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5874 if (fraction>0.95) break;
5875 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5876 e*fraction*fraction*fraction*fraction;
5877 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5878 sum += weight*projch->GetBinContent(k+1);
5879 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5880 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5882 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5884 if (fDebugLevel == 1) {
5885 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5889 fNumberFitSuccess++;
5890 CalculChargeCoefMean(kTRUE);
5892 //_____________________________________________________________________________
5893 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5896 // mean w methode for the gain factor
5900 Int_t nybins = projch->GetNbinsX();
5902 //The weight function
5903 Double_t a = 0.00228515;
5904 Double_t b = -0.00231487;
5905 Double_t c = 0.00044298;
5906 Double_t d = -0.00379239;
5907 Double_t e = 0.00338349;
5917 //A arbitrary error for the moment
5918 fCurrentCoefE = 0.0;
5919 fCurrentCoef[0] = 0.0;
5922 Double_t sumw = 0.0;
5924 Int_t sumCurrent = 0;
5925 for(Int_t k = 0; k <nybins; k++){
5926 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5927 if (fraction>0.95) break;
5928 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5929 e*fraction*fraction*fraction*fraction;
5930 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5931 sum += weight*projch->GetBinContent(k+1);
5932 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5933 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5935 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5937 if (fDebugLevel == 1) {
5938 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5942 fNumberFitSuccess++;
5944 //_____________________________________________________________________________
5945 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5948 // Fit methode for the gain factor
5951 fCurrentCoef[0] = 0.0;
5952 fCurrentCoefE = 0.0;
5953 Double_t chisqrl = 0.0;
5954 Double_t chisqrg = 0.0;
5955 Double_t chisqr = 0.0;
5956 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5958 projch->Fit("landau","0",""
5959 ,(Double_t) mean/fBeginFitCharge
5960 ,projch->GetBinCenter(projch->GetNbinsX()));
5961 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5962 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5963 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5964 chisqrl = projch->GetFunction("landau")->GetChisquare();
5966 projch->Fit("gaus","0",""
5967 ,(Double_t) mean/fBeginFitCharge
5968 ,projch->GetBinCenter(projch->GetNbinsX()));
5969 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5970 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5971 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5973 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5974 if (fDebugLevel != 1) {
5975 projch->Fit("fLandauGaus","0",""
5976 ,(Double_t) mean/fBeginFitCharge
5977 ,projch->GetBinCenter(projch->GetNbinsX()));
5978 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5981 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5983 projch->Fit("fLandauGaus","+",""
5984 ,(Double_t) mean/fBeginFitCharge
5985 ,projch->GetBinCenter(projch->GetNbinsX()));
5986 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5988 fLandauGaus->Draw("same");
5991 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5992 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5993 fNumberFitSuccess++;
5994 CalculChargeCoefMean(kTRUE);
5995 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5996 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5999 CalculChargeCoefMean(kFALSE);
6000 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6003 if (fDebugLevel != 1) {
6008 //_____________________________________________________________________________
6009 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
6012 // Fit methode for the gain factor more time consuming
6016 //Some parameters to initialise
6017 Double_t widthLandau, widthGaus, mPV, integral;
6018 Double_t chisquarel = 0.0;
6019 Double_t chisquareg = 0.0;
6020 projch->Fit("landau","0M+",""
6022 ,projch->GetBinCenter(projch->GetNbinsX()));
6023 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6024 chisquarel = projch->GetFunction("landau")->GetChisquare();
6025 projch->Fit("gaus","0M+",""
6027 ,projch->GetBinCenter(projch->GetNbinsX()));
6028 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6029 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6031 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6032 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
6034 // Setting fit range and start values
6036 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
6037 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
6038 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
6039 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
6040 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
6041 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
6042 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
6045 fCurrentCoef[0] = 0.0;
6046 fCurrentCoefE = 0.0;
6050 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
6055 Double_t projchPeak;
6056 Double_t projchFWHM;
6057 LanGauPro(fp,projchPeak,projchFWHM);
6059 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6060 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6061 fNumberFitSuccess++;
6062 CalculChargeCoefMean(kTRUE);
6063 fCurrentCoef[0] = fp[1];
6064 fCurrentCoefE = fpe[1];
6065 //chargeCoefE2 = chisqr;
6068 CalculChargeCoefMean(kFALSE);
6069 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6071 if (fDebugLevel == 1) {
6072 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
6073 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6076 fitsnr->Draw("same");
6082 //_____________________________________________________________________________
6083 void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6086 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
6088 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
6089 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
6090 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
6095 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
6096 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
6100 //_____________________________________________________________________________
6101 void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6104 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
6106 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
6107 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
6108 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
6109 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
6113 c2 = -(x0*(x[1]+x[2]+x[3])
6114 +x1*(x[0]+x[2]+x[3])
6115 +x2*(x[0]+x[1]+x[3])
6116 +x3*(x[0]+x[1]+x[2]));
6117 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
6118 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
6119 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
6120 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
6122 c0 = -(x0*x[1]*x[2]*x[3]
6125 +x3*x[0]*x[1]*x[2]);
6130 //_____________________________________________________________________________
6131 void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6134 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
6136 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
6137 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
6138 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
6139 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
6140 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
6143 c4 = x0+x1+x2+x3+x4;
6144 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
6145 +x1*(x[0]+x[2]+x[3]+x[4])
6146 +x2*(x[0]+x[1]+x[3]+x[4])
6147 +x3*(x[0]+x[1]+x[2]+x[4])
6148 +x4*(x[0]+x[1]+x[2]+x[3]));
6149 c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
6150 +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])
6151 +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])
6152 +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])
6153 +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]));
6155 c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
6156 +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])
6157 +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])
6158 +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])
6159 +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]));
6161 c0 = (x0*x[1]*x[2]*x[3]*x[4]
6162 +x1*x[0]*x[2]*x[3]*x[4]
6163 +x2*x[0]*x[1]*x[3]*x[4]
6164 +x3*x[0]*x[1]*x[2]*x[4]
6165 +x4*x[0]*x[1]*x[2]*x[3]);
6168 //_____________________________________________________________________________
6169 void AliTRDCalibraFit::NormierungCharge()
6172 // Normalisation of the gain factor resulting for the fits
6175 // Calcul of the mean of choosen method by fFitChargeNDB
6177 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
6178 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
6180 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
6181 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
6182 //printf("detector %d coef[0] %f\n",detector,coef[0]);
6183 if (GetStack(detector) == 2) {
6186 if (GetStack(detector) != 2) {
6189 for (Int_t j = 0; j < total; j++) {
6197 fScaleFitFactor = fScaleFitFactor / sum;
6200 fScaleFitFactor = 1.0;
6203 //methode de boeuf mais bon...
6204 Double_t scalefactor = fScaleFitFactor;
6206 if(fDebugLevel > 1){
6208 if ( !fDebugStreamer ) {
6210 TDirectory *backup = gDirectory;
6211 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
6212 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
6214 (* fDebugStreamer) << "NormierungCharge"<<
6215 "scalefactor="<<scalefactor<<
6219 //_____________________________________________________________________________
6220 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
6223 // Rebin of the 1D histo for the gain calibration if needed.
6224 // you have to choose fRebin, divider of fNumberBinCharge
6227 TAxis *xhist = hist->GetXaxis();
6228 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6229 ,xhist->GetBinLowEdge(1)
6230 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6232 AliInfo(Form("fRebin: %d",fRebin));
6234 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6236 for (Int_t ji = i; ji < i+fRebin; ji++) {
6237 sum += hist->GetBinContent(ji);
6240 rehist->SetBinContent(k,sum);
6248 //_____________________________________________________________________________
6249 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
6252 // Rebin of the 1D histo for the gain calibration if needed
6253 // you have to choose fRebin divider of fNumberBinCharge
6256 TAxis *xhist = hist->GetXaxis();
6257 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6258 ,xhist->GetBinLowEdge(1)
6259 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6261 AliInfo(Form("fRebin: %d",fRebin));
6263 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6265 for (Int_t ji = i; ji < i+fRebin; ji++) {
6266 sum += hist->GetBinContent(ji);
6269 rehist->SetBinContent(k,sum);
6277 //____________Some basic geometry function_____________________________________
6280 //_____________________________________________________________________________
6281 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
6284 // Reconstruct the plane number from the detector number
6287 return ((Int_t) (d % 6));
6291 //_____________________________________________________________________________
6292 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
6295 // Reconstruct the stack number from the detector number
6297 const Int_t kNlayer = 6;
6299 return ((Int_t) (d % 30) / kNlayer);
6303 //_____________________________________________________________________________
6304 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6307 // Reconstruct the sector number from the detector number
6311 return ((Int_t) (d / fg));
6316 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6318 //_______________________________________________________________________________
6319 void AliTRDCalibraFit::ResetVectorFit()
6322 // Reset the VectorFits
6325 fVectorFit.SetOwner();
6327 fVectorFit2.SetOwner();
6328 fVectorFit2.Clear();
6332 //____________Private Functions________________________________________________
6335 //_____________________________________________________________________________
6336 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
6339 // Function for the fit
6342 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6344 //PARAMETERS FOR FIT PH
6346 //fAsymmGauss->SetParameter(0,0.113755);
6347 //fAsymmGauss->SetParameter(1,0.350706);
6348 //fAsymmGauss->SetParameter(2,0.0604244);
6349 //fAsymmGauss->SetParameter(3,7.65596);
6350 //fAsymmGauss->SetParameter(4,1.00124);
6351 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6359 Double_t dx = 0.005;
6360 Double_t xs = par[1];
6362 Double_t paras[2] = { 0.0, 0.0 };
6365 if ((xs >= par[1]) &&
6366 (xs < (par[1]+par[2]))) {
6367 //fAsymmGauss->SetParameter(0,par[0]);
6368 //fAsymmGauss->SetParameter(1,xs);
6369 //ss += fAsymmGauss->Eval(xx);
6372 ss += AsymmGauss(&xx,paras);
6374 if ((xs >= (par[1]+par[2])) &&
6375 (xs < (par[1]+par[2]+par[3]))) {
6376 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6377 //fAsymmGauss->SetParameter(1,xs);
6378 //ss += fAsymmGauss->Eval(xx);
6379 paras[0] = par[0]*par[4];
6381 ss += AsymmGauss(&xx,paras);
6390 //_____________________________________________________________________________
6391 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
6394 // Function for the fit
6397 //par[0] = normalization
6405 Double_t par1save = par[1];
6406 //Double_t par2save = par[2];
6407 Double_t par2save = 0.0604244;
6408 //Double_t par3save = par[3];
6409 Double_t par3save = 7.65596;
6410 //Double_t par5save = par[5];
6411 Double_t par5save = 0.870597;
6412 Double_t dx = x[0] - par1save;
6414 Double_t sigma2 = par2save*par2save;
6415 Double_t sqrt2 = TMath::Sqrt(2.0);
6416 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6417 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6418 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6419 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6421 //return par[0]*(exp1+par[4]*exp2);
6422 return par[0] * (exp1 + 1.00124 * exp2);
6426 //_____________________________________________________________________________
6427 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
6430 // Sum Landau + Gaus with identical mean
6433 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6434 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6435 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6436 Double_t val = valLandau + valGaus;
6442 //_____________________________________________________________________________
6443 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6446 // Function for the fit
6449 // par[0]=Width (scale) parameter of Landau density
6450 // par[1]=Most Probable (MP, location) parameter of Landau density
6451 // par[2]=Total area (integral -inf to inf, normalization constant)
6452 // par[3]=Width (sigma) of convoluted Gaussian function
6454 // In the Landau distribution (represented by the CERNLIB approximation),
6455 // the maximum is located at x=-0.22278298 with the location parameter=0.
6456 // This shift is corrected within this function, so that the actual
6457 // maximum is identical to the MP parameter.
6460 // Numeric constants
6461 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6462 Double_t mpshift = -0.22278298; // Landau maximum location
6464 // Control constants
6465 Double_t np = 100.0; // Number of convolution steps
6466 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6478 // MP shift correction
6479 mpc = par[1] - mpshift * par[0];
6481 // Range of convolution integral
6482 xlow = x[0] - sc * par[3];
6483 xupp = x[0] + sc * par[3];
6485 step = (xupp - xlow) / np;
6487 // Convolution integral of Landau and Gaussian by sum
6488 for (i = 1.0; i <= np/2; i++) {
6490 xx = xlow + (i-.5) * step;
6491 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6492 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6494 xx = xupp - (i-.5) * step;
6495 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6496 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6500 return (par[2] * step * sum * invsq2pi / par[3]);
6503 //_____________________________________________________________________________
6504 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6505 , const Double_t *parlimitslo, const Double_t *parlimitshi
6506 , Double_t *fitparams, Double_t *fiterrors
6507 , Double_t *chiSqr, Int_t *ndf) const
6510 // Function for the fit
6514 Char_t funname[100];
6516 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6521 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6522 ffit->SetParameters(startvalues);
6523 ffit->SetParNames("Width","MP","Area","GSigma");
6525 for (i = 0; i < 4; i++) {
6526 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6529 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6531 ffit->GetParameters(fitparams); // Obtain fit parameters
6532 for (i = 0; i < 4; i++) {
6533 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6535 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6536 ndf[0] = ffit->GetNDF(); // Obtain ndf
6538 return (ffit); // Return fit function
6542 //_____________________________________________________________________________
6543 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6546 // Function for the fit
6559 Int_t maxcalls = 10000;
6561 // Search for maximum
6562 p = params[1] - 0.1 * params[0];
6563 step = 0.05 * params[0];
6567 while ((l != lold) && (i < maxcalls)) {
6571 l = LanGauFun(&x,params);
6573 step = -step / 10.0;
6578 if (i == maxcalls) {
6584 // Search for right x location of fy
6585 p = maxx + params[0];
6591 while ( (l != lold) && (i < maxcalls) ) {
6596 l = TMath::Abs(LanGauFun(&x,params) - fy);
6610 // Search for left x location of fy
6612 p = maxx - 0.5 * params[0];
6618 while ((l != lold) && (i < maxcalls)) {
6622 l = TMath::Abs(LanGauFun(&x,params) - fy);
6624 step = -step / 10.0;
6629 if (i == maxcalls) {
6638 //_____________________________________________________________________________
6639 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6642 // Gaus with identical mean
6645 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];