1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
63 #include <TTreeStream.h>
64 #include <TLinearFitter.h>
69 #include "AliMathBase.h"
71 #include "AliTRDCalibraFit.h"
72 #include "AliTRDCalibraMode.h"
73 #include "AliTRDCalibraVector.h"
74 #include "AliTRDCalibraVdriftLinearFit.h"
75 #include "AliTRDcalibDB.h"
76 #include "AliTRDgeometry.h"
77 #include "AliTRDpadPlane.h"
78 #include "AliTRDgeometry.h"
79 #include "AliTRDCommonParam.h"
80 #include "./Cal/AliTRDCalROC.h"
81 #include "./Cal/AliTRDCalPad.h"
82 #include "./Cal/AliTRDCalDet.h"
85 ClassImp(AliTRDCalibraFit)
87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
90 //_____________singleton implementation_________________________________________________
91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
94 // Singleton implementation
97 if (fgTerminated != kFALSE) {
101 if (fgInstance == 0) {
102 fgInstance = new AliTRDCalibraFit();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFit::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
124 //______________________________________________________________________________________
125 AliTRDCalibraFit::AliTRDCalibraFit()
128 ,fNumberOfBinsExpected(0)
130 ,fBeginFitCharge(3.5)
132 ,fTakeTheMaxPH(kTRUE)
140 ,fNumberFitSuccess(0)
147 ,fCalibraMode(new AliTRDCalibraMode())
152 ,fScaleFitFactor(0.0)
161 ,fCurrentCoefDetector(0x0)
162 ,fCurrentCoefDetector2(0x0)
167 // Default constructor
170 fGeo = new AliTRDgeometry();
172 // Current variables initialised
173 for (Int_t k = 0; k < 2; k++) {
174 fCurrentCoef[k] = 0.0;
175 fCurrentCoef2[k] = 0.0;
177 for (Int_t i = 0; i < 3; i++) {
183 //______________________________________________________________________________________
184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
189 ,fBeginFitCharge(c.fBeginFitCharge)
190 ,fFitPHPeriode(c.fFitPHPeriode)
191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
192 ,fT0Shift0(c.fT0Shift0)
193 ,fT0Shift1(c.fT0Shift1)
194 ,fRangeFitPRF(c.fRangeFitPRF)
196 ,fMinEntries(c.fMinEntries)
198 ,fNumberFit(c.fNumberFit)
199 ,fNumberFitSuccess(c.fNumberFitSuccess)
200 ,fNumberEnt(c.fNumberEnt)
201 ,fStatisticMean(c.fStatisticMean)
203 ,fDebugLevel(c.fDebugLevel)
204 ,fFitVoir(c.fFitVoir)
205 ,fMagneticField(c.fMagneticField)
207 ,fCurrentCoefE(c.fCurrentCoefE)
208 ,fCurrentCoefE2(c.fCurrentCoefE2)
211 ,fScaleFitFactor(c.fScaleFitFactor)
212 ,fEntriesCurrent(c.fEntriesCurrent)
213 ,fCountDet(c.fCountDet)
220 ,fCurrentCoefDetector(0x0)
221 ,fCurrentCoefDetector2(0x0)
229 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
231 //Current variables initialised
232 for (Int_t k = 0; k < 2; k++) {
233 fCurrentCoef[k] = 0.0;
234 fCurrentCoef2[k] = 0.0;
236 for (Int_t i = 0; i < 3; i++) {
240 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
241 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
243 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
244 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
246 fVectorFit.SetName(c.fVectorFit.GetName());
247 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
248 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
249 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
251 if (GetStack(detector) == 2) {
257 Float_t *coef = new Float_t[ntotal];
258 for (Int_t i = 0; i < ntotal; i++) {
259 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
261 fitInfo->SetCoef(coef);
262 fitInfo->SetDetector(detector);
263 fVectorFit.Add((TObject *) fitInfo);
265 fVectorFit.SetName(c.fVectorFit.GetName());
266 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
267 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
268 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
270 if (GetStack(detector) == 2) {
276 Float_t *coef = new Float_t[ntotal];
277 for (Int_t i = 0; i < ntotal; i++) {
278 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
280 fitInfo->SetCoef(coef);
281 fitInfo->SetDetector(detector);
282 fVectorFit2.Add((TObject *) fitInfo);
287 fGeo = new AliTRDgeometry();
290 //____________________________________________________________________________________
291 AliTRDCalibraFit::~AliTRDCalibraFit()
294 // AliTRDCalibraFit destructor
296 if ( fDebugStreamer ) delete fDebugStreamer;
297 if ( fCalDet ) delete fCalDet;
298 if ( fCalDet2 ) delete fCalDet2;
299 if ( fCalROC ) delete fCalROC;
300 if ( fCalROC2 ) delete fCalROC2;
301 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
302 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
304 fVectorFit2.Delete();
310 //_____________________________________________________________________________
311 void AliTRDCalibraFit::Destroy()
323 //_____________________________________________________________________________
324 void AliTRDCalibraFit::DestroyDebugStreamer()
327 // Delete DebugStreamer
330 if ( fDebugStreamer ) delete fDebugStreamer;
331 fDebugStreamer = 0x0;
334 //__________________________________________________________________________________
335 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
338 // From the drift velocity and t0
339 // return the position of the peak and maximum negative slope
342 const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
343 Double_t widbins = 0.1; // 0.1 mus
345 //peak and maxnegslope in mus
346 Double_t begind = t0*widbins + fT0Shift0;
347 Double_t peakd = t0*widbins + fT0Shift1;
348 Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
350 // peak and maxnegslope in timebin
351 begin = TMath::Nint(begind*widbins);
352 peak = TMath::Nint(peakd*widbins);
353 end = TMath::Nint(maxnegslope*widbins);
356 //____________Functions fit Online CH2d________________________________________
357 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
360 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
361 // calibration group normalized the resulted coefficients (to 1 normally)
364 // Set the calibration mode
365 const char *name = ch->GetTitle();
366 if(!SetModeCalibration(name,0)) return kFALSE;
368 // Number of Ybins (detectors or groups of pads)
369 Int_t nbins = ch->GetNbinsX();// charge
370 Int_t nybins = ch->GetNbinsY();// groups number
371 if (!InitFit(nybins,0)) {
377 fStatisticMean = 0.0;
379 fNumberFitSuccess = 0;
381 // Init fCountDet and fCount
382 InitfCountDetAndfCount(0);
383 // Beginning of the loop betwwen dect1 and dect2
384 for (Int_t idect = fDect1; idect < fDect2; idect++) {
385 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
386 UpdatefCountDetAndfCount(idect,0);
387 ReconstructFitRowMinRowMax(idect, 0);
389 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
390 projch->SetDirectory(0);
391 // Number of entries for this calibration group
392 Double_t nentries = 0.0;
394 for (Int_t k = 0; k < nbins; k++) {
395 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
396 nentries += ch->GetBinContent(binnb);
397 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
398 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
400 projch->SetEntries(nentries);
401 //printf("The number of entries for the group %d is %f\n",idect,nentries);
406 // Rebin and statistic stuff
408 projch = ReBin((TH1I *) projch);
410 // This detector has not enough statistics or was off
411 if (nentries <= fMinEntries) {
412 NotEnoughStatisticCH(idect);
413 if (fDebugLevel != 1) {
418 // Statistics of the group fitted
419 fStatisticMean += nentries;
424 case 0: FitMeanW((TH1 *) projch, nentries); break;
425 case 1: FitMean((TH1 *) projch, nentries, mean); break;
426 case 2: FitCH((TH1 *) projch, mean); break;
427 case 3: FitBisCH((TH1 *) projch, mean); break;
428 default: return kFALSE;
431 FillInfosFitCH(idect);
433 if (fDebugLevel != 1) {
438 if (fDebugLevel != 1) {
442 if (fNumberFit > 0) {
443 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
444 fStatisticMean = fStatisticMean / fNumberFit;
447 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
449 delete fDebugStreamer;
450 fDebugStreamer = 0x0;
454 //____________Functions fit Online CH2d________________________________________
455 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
458 // Reconstruct a 1D histo from the vectorCH for each calibration group,
459 // fit the histo, normalized the resulted coefficients (to 1 normally)
462 // Set the calibraMode
463 const char *name = calvect->GetNameCH();
464 if(!SetModeCalibration(name,0)) return kFALSE;
466 // Number of Xbins (detectors or groups of pads)
467 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
473 fStatisticMean = 0.0;
475 fNumberFitSuccess = 0;
477 // Init fCountDet and fCount
478 InitfCountDetAndfCount(0);
479 // Beginning of the loop between dect1 and dect2
480 for (Int_t idect = fDect1; idect < fDect2; idect++) {
481 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
482 UpdatefCountDetAndfCount(idect,0);
483 ReconstructFitRowMinRowMax(idect,0);
485 Double_t nentries = 0.0;
487 if(!calvect->GetCHEntries(fCountDet)) {
488 NotEnoughStatisticCH(idect);
494 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
495 projch->SetDirectory(0);
496 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
497 nentries += projch->GetBinContent(k+1);
498 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
504 //printf("The number of entries for the group %d is %f\n",idect,nentries);
507 projch = ReBin((TH1F *) projch);
509 // This detector has not enough statistics or was not found in VectorCH
510 if (nentries <= fMinEntries) {
511 NotEnoughStatisticCH(idect);
514 // Statistic of the histos fitted
515 fStatisticMean += nentries;
520 case 0: FitMeanW((TH1 *) projch, nentries); break;
521 case 1: FitMean((TH1 *) projch, nentries, mean); break;
522 case 2: FitCH((TH1 *) projch, mean); break;
523 case 3: FitBisCH((TH1 *) projch, mean); break;
524 default: return kFALSE;
527 FillInfosFitCH(idect);
530 if (fDebugLevel != 1) {
534 if (fNumberFit > 0) {
535 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
536 fStatisticMean = fStatisticMean / fNumberFit;
539 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
541 delete fDebugStreamer;
542 fDebugStreamer = 0x0;
545 //________________functions fit Online PH2d____________________________________
546 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
549 // Take the 1D profiles (average pulse height), projections of the 2D PH
550 // on the Xaxis, for each calibration group
551 // Reconstruct a drift velocity
552 // A first calibration of T0 is also made using the same method
555 // Set the calibration mode
556 const char *name = ph->GetTitle();
557 if(!SetModeCalibration(name,1)) return kFALSE;
559 //printf("Mode calibration set\n");
561 // Number of Xbins (detectors or groups of pads)
562 Int_t nbins = ph->GetNbinsX();// time
563 Int_t nybins = ph->GetNbinsY();// calibration group
564 if (!InitFit(nybins,1)) {
568 //printf("Init fit\n");
574 //printf("Init fit PH\n");
576 fStatisticMean = 0.0;
578 fNumberFitSuccess = 0;
580 // Init fCountDet and fCount
581 InitfCountDetAndfCount(1);
582 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
584 // Beginning of the loop
585 for (Int_t idect = fDect1; idect < fDect2; idect++) {
586 //printf("idect = %d\n",idect);
587 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
588 UpdatefCountDetAndfCount(idect,1);
589 ReconstructFitRowMinRowMax(idect,1);
591 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
592 projph->SetDirectory(0);
593 // Number of entries for this calibration group
594 Double_t nentries = 0;
595 for (Int_t k = 0; k < nbins; k++) {
596 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
597 nentries += ph->GetBinEntries(binnb);
602 //printf("The number of entries for the group %d is %f\n",idect,nentries);
603 // This detector has not enough statistics or was off
604 if (nentries <= fMinEntries) {
605 //printf("Not enough statistic!\n");
606 NotEnoughStatisticPH(idect,nentries);
607 if (fDebugLevel != 1) {
612 // Statistics of the histos fitted
614 fStatisticMean += nentries;
615 // Calcul of "real" coef
616 CalculVdriftCoefMean();
619 //printf("Method\n");
622 case 0: FitLagrangePoly((TH1 *) projph); break;
623 case 1: FitPente((TH1 *) projph); break;
624 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
625 default: return kFALSE;
627 // Fill the tree if end of a detector or only the pointer to the branch!!!
628 FillInfosFitPH(idect,nentries);
630 if (fDebugLevel != 1) {
635 if (fNumberFit > 0) {
636 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));
637 fStatisticMean = fStatisticMean / fNumberFit;
640 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
642 delete fDebugStreamer;
643 fDebugStreamer = 0x0;
646 //____________Functions fit Online PH2d________________________________________
647 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
650 // Reconstruct the average pulse height from the vectorPH for each
652 // Reconstruct a drift velocity
653 // A first calibration of T0 is also made using the same method (slope method)
656 // Set the calibration mode
657 const char *name = calvect->GetNamePH();
658 if(!SetModeCalibration(name,1)) return kFALSE;
660 // Number of Xbins (detectors or groups of pads)
661 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
667 fStatisticMean = 0.0;
669 fNumberFitSuccess = 0;
671 // Init fCountDet and fCount
672 InitfCountDetAndfCount(1);
673 // Beginning of the loop
674 for (Int_t idect = fDect1; idect < fDect2; idect++) {
675 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
676 UpdatefCountDetAndfCount(idect,1);
677 ReconstructFitRowMinRowMax(idect,1);
680 if(!calvect->GetPHEntries(fCountDet)) {
681 NotEnoughStatisticPH(idect,fEntriesCurrent);
686 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
687 projph->SetDirectory(0);
688 if(fEntriesCurrent > 0) fNumberEnt++;
689 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
690 // This detector has not enough statistics or was off
691 if (fEntriesCurrent <= fMinEntries) {
692 //printf("Not enough stat!\n");
693 NotEnoughStatisticPH(idect,fEntriesCurrent);
696 // Statistic of the histos fitted
698 fStatisticMean += fEntriesCurrent;
699 // Calcul of "real" coef
700 CalculVdriftCoefMean();
705 case 0: FitLagrangePoly((TH1 *) projph); break;
706 case 1: FitPente((TH1 *) projph); break;
707 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
708 default: return kFALSE;
710 // Fill the tree if end of a detector or only the pointer to the branch!!!
711 FillInfosFitPH(idect,fEntriesCurrent);
715 if (fNumberFit > 0) {
716 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));
717 fStatisticMean = fStatisticMean / fNumberFit;
720 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
722 delete fDebugStreamer;
723 fDebugStreamer = 0x0;
726 //____________Functions fit Online PRF2d_______________________________________
727 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
730 // Take the 1D profiles (pad response function), projections of the 2D PRF
731 // on the Xaxis, for each calibration group
732 // Fit with a gaussian to reconstruct the sigma of the pad response function
735 // Set the calibration mode
736 const char *name = prf->GetTitle();
737 if(!SetModeCalibration(name,2)) return kFALSE;
739 // Number of Ybins (detectors or groups of pads)
740 Int_t nybins = prf->GetNbinsY();// calibration groups
741 Int_t nbins = prf->GetNbinsX();// bins
742 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
743 if((nbg > 0) || (nbg == -1)) return kFALSE;
744 if (!InitFit(nybins,2)) {
750 fStatisticMean = 0.0;
752 fNumberFitSuccess = 0;
754 // Init fCountDet and fCount
755 InitfCountDetAndfCount(2);
756 // Beginning of the loop
757 for (Int_t idect = fDect1; idect < fDect2; idect++) {
758 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
759 UpdatefCountDetAndfCount(idect,2);
760 ReconstructFitRowMinRowMax(idect,2);
762 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
763 projprf->SetDirectory(0);
764 // Number of entries for this calibration group
765 Double_t nentries = 0;
766 for (Int_t k = 0; k < nbins; k++) {
767 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
768 nentries += prf->GetBinEntries(binnb);
770 if(nentries > 0) fNumberEnt++;
771 // This detector has not enough statistics or was off
772 if (nentries <= fMinEntries) {
773 NotEnoughStatisticPRF(idect);
774 if (fDebugLevel != 1) {
779 // Statistics of the histos fitted
781 fStatisticMean += nentries;
782 // Calcul of "real" coef
787 case 0: FitPRF((TH1 *) projprf); break;
788 case 1: RmsPRF((TH1 *) projprf); break;
789 default: return kFALSE;
791 // Fill the tree if end of a detector or only the pointer to the branch!!!
792 FillInfosFitPRF(idect);
794 if (fDebugLevel != 1) {
799 if (fNumberFit > 0) {
800 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
801 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
802 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
803 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
804 fStatisticMean = fStatisticMean / fNumberFit;
807 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
809 delete fDebugStreamer;
810 fDebugStreamer = 0x0;
813 //____________Functions fit Online PRF2d_______________________________________
814 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
817 // Take the 1D profiles (pad response function), projections of the 2D PRF
818 // on the Xaxis, for each calibration group
819 // Fit with a gaussian to reconstruct the sigma of the pad response function
822 // Set the calibration mode
823 const char *name = prf->GetTitle();
824 if(!SetModeCalibration(name,2)) return kFALSE;
826 // Number of Ybins (detectors or groups of pads)
827 TAxis *xprf = prf->GetXaxis();
828 TAxis *yprf = prf->GetYaxis();
829 Int_t nybins = yprf->GetNbins();// calibration groups
830 Int_t nbins = xprf->GetNbins();// bins
831 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
832 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
833 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
834 if(nbg == -1) return kFALSE;
835 if(nbg > 0) fMethod = 1;
837 if (!InitFit(nybins,2)) {
843 fStatisticMean = 0.0;
845 fNumberFitSuccess = 0;
847 // Init fCountDet and fCount
848 InitfCountDetAndfCount(2);
849 // Beginning of the loop
850 for (Int_t idect = fDect1; idect < fDect2; idect++) {
851 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
852 UpdatefCountDetAndfCount(idect,2);
853 ReconstructFitRowMinRowMax(idect,2);
854 // Build the array of entries and sum
855 TArrayD arraye = TArrayD(nbins);
856 TArrayD arraym = TArrayD(nbins);
857 TArrayD arrayme = TArrayD(nbins);
858 Double_t nentries = 0;
859 //printf("nbins %d\n",nbins);
860 for (Int_t k = 0; k < nbins; k++) {
861 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
862 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
863 Double_t mean = (Double_t)prf->GetBinContent(binnb);
864 Double_t error = (Double_t)prf->GetBinError(binnb);
865 //printf("for %d we have %f\n",k,entries);
867 arraye.AddAt(entries,k);
868 arraym.AddAt(mean,k);
869 arrayme.AddAt(error,k);
871 if(nentries > 0) fNumberEnt++;
872 //printf("The number of entries for the group %d is %f\n",idect,nentries);
873 // This detector has not enough statistics or was off
874 if (nentries <= fMinEntries) {
875 NotEnoughStatisticPRF(idect);
878 // Statistics of the histos fitted
880 fStatisticMean += nentries;
881 // Calcul of "real" coef
886 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
887 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
888 default: return kFALSE;
890 // Fill the tree if end of a detector or only the pointer to the branch!!!
891 FillInfosFitPRF(idect);
894 if (fNumberFit > 0) {
895 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
896 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
897 AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
898 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
899 fStatisticMean = fStatisticMean / fNumberFit;
902 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
904 delete fDebugStreamer;
905 fDebugStreamer = 0x0;
908 //____________Functions fit Online PRF2d_______________________________________
909 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
912 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
913 // each calibration group
914 // Fit with a gaussian to reconstruct the sigma of the pad response function
917 // Set the calibra mode
918 const char *name = calvect->GetNamePRF();
919 if(!SetModeCalibration(name,2)) return kFALSE;
920 //printf("test0 %s\n",name);
922 // Number of Xbins (detectors or groups of pads)
923 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
928 ///printf("test2\n");
931 fStatisticMean = 0.0;
933 fNumberFitSuccess = 0;
935 // Init fCountDet and fCount
936 InitfCountDetAndfCount(2);
937 // Beginning of the loop
938 for (Int_t idect = fDect1; idect < fDect2; idect++) {
939 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
940 UpdatefCountDetAndfCount(idect,2);
941 ReconstructFitRowMinRowMax(idect,2);
944 if(!calvect->GetPRFEntries(fCountDet)) {
945 NotEnoughStatisticPRF(idect);
948 TString tname("PRF");
950 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
951 projprf->SetDirectory(0);
952 if(fEntriesCurrent > 0) fNumberEnt++;
953 // This detector has not enough statistics or was off
954 if (fEntriesCurrent <= fMinEntries) {
955 NotEnoughStatisticPRF(idect);
958 // Statistic of the histos fitted
960 fStatisticMean += fEntriesCurrent;
961 // Calcul of "real" coef
966 case 1: FitPRF((TH1 *) projprf); break;
967 case 2: RmsPRF((TH1 *) projprf); break;
968 default: return kFALSE;
970 // Fill the tree if end of a detector or only the pointer to the branch!!!
971 FillInfosFitPRF(idect);
974 if (fNumberFit > 0) {
975 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));
978 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
980 delete fDebugStreamer;
981 fDebugStreamer = 0x0;
984 //____________Functions fit Online PRF2d_______________________________________
985 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
988 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
989 // each calibration group
990 // Fit with a gaussian to reconstruct the sigma of the pad response function
993 // Set the calibra mode
994 const char *name = calvect->GetNamePRF();
995 if(!SetModeCalibration(name,2)) return kFALSE;
996 //printf("test0 %s\n",name);
997 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
998 //printf("test1 %d\n",nbg);
999 if(nbg == -1) return kFALSE;
1000 if(nbg > 0) fMethod = 1;
1002 // Number of Xbins (detectors or groups of pads)
1003 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1004 //printf("test2\n");
1007 if (!InitFitPRF()) {
1008 //printf("test3\n");
1011 fStatisticMean = 0.0;
1013 fNumberFitSuccess = 0;
1017 Double_t *arrayx = 0;
1018 Double_t *arraye = 0;
1019 Double_t *arraym = 0;
1020 Double_t *arrayme = 0;
1021 Float_t lowedge = 0.0;
1022 Float_t upedge = 0.0;
1023 // Init fCountDet and fCount
1024 InitfCountDetAndfCount(2);
1025 // Beginning of the loop
1026 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1027 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1028 UpdatefCountDetAndfCount(idect,2);
1029 ReconstructFitRowMinRowMax(idect,2);
1031 fEntriesCurrent = 0;
1032 if(!calvect->GetPRFEntries(fCountDet)) {
1033 NotEnoughStatisticPRF(idect);
1036 TString tname("PRF");
1038 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1039 nbins = projprftree->GetN();
1040 arrayx = (Double_t *)projprftree->GetX();
1041 arraye = (Double_t *)projprftree->GetEX();
1042 arraym = (Double_t *)projprftree->GetY();
1043 arrayme = (Double_t *)projprftree->GetEY();
1044 Float_t step = arrayx[1]-arrayx[0];
1045 lowedge = arrayx[0] - step/2.0;
1046 upedge = arrayx[(nbins-1)] + step/2.0;
1047 //printf("nbins est %d\n",nbins);
1048 for(Int_t k = 0; k < nbins; k++){
1049 fEntriesCurrent += (Int_t)arraye[k];
1050 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1051 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1053 if(fEntriesCurrent > 0) fNumberEnt++;
1054 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1055 // This detector has not enough statistics or was off
1056 if (fEntriesCurrent <= fMinEntries) {
1057 NotEnoughStatisticPRF(idect);
1060 // Statistic of the histos fitted
1062 fStatisticMean += fEntriesCurrent;
1063 // Calcul of "real" coef
1064 CalculPRFCoefMean();
1068 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1069 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1070 default: return kFALSE;
1072 // Fill the tree if end of a detector or only the pointer to the branch!!!
1073 FillInfosFitPRF(idect);
1076 if (fNumberFit > 0) {
1077 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));
1080 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1082 delete fDebugStreamer;
1083 fDebugStreamer = 0x0;
1086 //____________Functions fit Online CH2d________________________________________
1087 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1090 // The linear method
1093 fStatisticMean = 0.0;
1095 fNumberFitSuccess = 0;
1097 if(!InitFitLinearFitter()) return kFALSE;
1100 for(Int_t idet = 0; idet < 540; idet++){
1103 //printf("detector number %d\n",idet);
1108 fEntriesCurrent = 0;
1110 Bool_t here = calivdli->GetParam(idet,¶m);
1111 Bool_t heree = calivdli->GetError(idet,&error);
1112 //printf("here %d and heree %d\n",here, heree);
1114 fEntriesCurrent = (Int_t) error[2];
1117 //printf("Number of entries %d\n",fEntriesCurrent);
1118 // Nothing found or not enough statistic
1119 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1120 NotEnoughStatisticLinearFitter();
1127 fStatisticMean += fEntriesCurrent;
1130 if((-(param[1])) <= 0.0) {
1131 NotEnoughStatisticLinearFitter();
1135 // CalculDatabaseVdriftandTan
1136 CalculVdriftLorentzCoef();
1139 fNumberFitSuccess ++;
1141 // Put the fCurrentCoef
1142 fCurrentCoef[0] = -param[1];
1143 // here the database must be the one of the reconstruction for the lorentz angle....
1144 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1145 fCurrentCoefE = error[1];
1146 fCurrentCoefE2 = error[0];
1147 if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1148 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1152 FillInfosFitLinearFitter();
1157 if (fNumberFit > 0) {
1158 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));
1161 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1163 delete fDebugStreamer;
1164 fDebugStreamer = 0x0;
1168 //____________Functions for seeing if the pad is really okey___________________
1169 //_____________________________________________________________________________
1170 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1173 // Get numberofgroupsprf
1177 const Char_t *pattern0 = "Ngp0";
1178 const Char_t *pattern1 = "Ngp1";
1179 const Char_t *pattern2 = "Ngp2";
1180 const Char_t *pattern3 = "Ngp3";
1181 const Char_t *pattern4 = "Ngp4";
1182 const Char_t *pattern5 = "Ngp5";
1183 const Char_t *pattern6 = "Ngp6";
1186 if (strstr(nametitle,pattern0)) {
1189 if (strstr(nametitle,pattern1)) {
1192 if (strstr(nametitle,pattern2)) {
1195 if (strstr(nametitle,pattern3)) {
1198 if (strstr(nametitle,pattern4)) {
1201 if (strstr(nametitle,pattern5)) {
1204 if (strstr(nametitle,pattern6)){
1211 //_____________________________________________________________________________
1212 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1215 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1216 // corresponding to the given name
1219 if(!SetNzFromTObject(name,i)) return kFALSE;
1220 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1225 //_____________________________________________________________________________
1226 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1229 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1230 // corresponding to the given TObject
1234 const Char_t *patternrphi0 = "Nrphi0";
1235 const Char_t *patternrphi1 = "Nrphi1";
1236 const Char_t *patternrphi2 = "Nrphi2";
1237 const Char_t *patternrphi3 = "Nrphi3";
1238 const Char_t *patternrphi4 = "Nrphi4";
1239 const Char_t *patternrphi5 = "Nrphi5";
1240 const Char_t *patternrphi6 = "Nrphi6";
1243 const Char_t *patternrphi10 = "Nrphi10";
1244 const Char_t *patternrphi100 = "Nrphi100";
1245 const Char_t *patternz10 = "Nz10";
1246 const Char_t *patternz100 = "Nz100";
1249 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1250 fCalibraMode->SetAllTogether(i);
1252 if (fDebugLevel > 1) {
1253 AliInfo(Form("fNbDet %d and 100",fNbDet));
1257 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1258 fCalibraMode->SetPerSuperModule(i);
1260 if (fDebugLevel > 1) {
1261 AliInfo(Form("fNDet %d and 100",fNbDet));
1266 if (strstr(name,patternrphi0)) {
1267 fCalibraMode->SetNrphi(i ,0);
1268 if (fDebugLevel > 1) {
1269 AliInfo(Form("fNbDet %d and 0",fNbDet));
1273 if (strstr(name,patternrphi1)) {
1274 fCalibraMode->SetNrphi(i, 1);
1275 if (fDebugLevel > 1) {
1276 AliInfo(Form("fNbDet %d and 1",fNbDet));
1280 if (strstr(name,patternrphi2)) {
1281 fCalibraMode->SetNrphi(i, 2);
1282 if (fDebugLevel > 1) {
1283 AliInfo(Form("fNbDet %d and 2",fNbDet));
1287 if (strstr(name,patternrphi3)) {
1288 fCalibraMode->SetNrphi(i, 3);
1289 if (fDebugLevel > 1) {
1290 AliInfo(Form("fNbDet %d and 3",fNbDet));
1294 if (strstr(name,patternrphi4)) {
1295 fCalibraMode->SetNrphi(i, 4);
1296 if (fDebugLevel > 1) {
1297 AliInfo(Form("fNbDet %d and 4",fNbDet));
1301 if (strstr(name,patternrphi5)) {
1302 fCalibraMode->SetNrphi(i, 5);
1303 if (fDebugLevel > 1) {
1304 AliInfo(Form("fNbDet %d and 5",fNbDet));
1308 if (strstr(name,patternrphi6)) {
1309 fCalibraMode->SetNrphi(i, 6);
1310 if (fDebugLevel > 1) {
1311 AliInfo(Form("fNbDet %d and 6",fNbDet));
1316 if (fDebugLevel > 1) {
1317 AliInfo(Form("fNbDet %d and rest",fNbDet));
1319 fCalibraMode->SetNrphi(i ,0);
1323 //_____________________________________________________________________________
1324 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1327 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1328 // corresponding to the given TObject
1332 const Char_t *patternz0 = "Nz0";
1333 const Char_t *patternz1 = "Nz1";
1334 const Char_t *patternz2 = "Nz2";
1335 const Char_t *patternz3 = "Nz3";
1336 const Char_t *patternz4 = "Nz4";
1338 const Char_t *patternrphi10 = "Nrphi10";
1339 const Char_t *patternrphi100 = "Nrphi100";
1340 const Char_t *patternz10 = "Nz10";
1341 const Char_t *patternz100 = "Nz100";
1343 if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1344 fCalibraMode->SetAllTogether(i);
1346 if (fDebugLevel > 1) {
1347 AliInfo(Form("fNbDet %d and 100",fNbDet));
1351 if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1352 fCalibraMode->SetPerSuperModule(i);
1354 if (fDebugLevel > 1) {
1355 AliInfo(Form("fNbDet %d and 10",fNbDet));
1359 if (strstr(name,patternz0)) {
1360 fCalibraMode->SetNz(i, 0);
1361 if (fDebugLevel > 1) {
1362 AliInfo(Form("fNbDet %d and 0",fNbDet));
1366 if (strstr(name,patternz1)) {
1367 fCalibraMode->SetNz(i ,1);
1368 if (fDebugLevel > 1) {
1369 AliInfo(Form("fNbDet %d and 1",fNbDet));
1373 if (strstr(name,patternz2)) {
1374 fCalibraMode->SetNz(i ,2);
1375 if (fDebugLevel > 1) {
1376 AliInfo(Form("fNbDet %d and 2",fNbDet));
1380 if (strstr(name,patternz3)) {
1381 fCalibraMode->SetNz(i ,3);
1382 if (fDebugLevel > 1) {
1383 AliInfo(Form("fNbDet %d and 3",fNbDet));
1387 if (strstr(name,patternz4)) {
1388 fCalibraMode->SetNz(i ,4);
1389 if (fDebugLevel > 1) {
1390 AliInfo(Form("fNbDet %d and 4",fNbDet));
1395 if (fDebugLevel > 1) {
1396 AliInfo(Form("fNbDet %d and rest",fNbDet));
1398 fCalibraMode->SetNz(i ,0);
1401 //______________________________________________________________________
1402 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1404 // Remove the results too far from the mean value and rms
1405 // type: 0 gain, 1 vdrift
1409 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1411 AliInfo("The Vector Fit is not complete!");
1414 Int_t detector = -1;
1416 Float_t value = 0.0;
1418 /////////////////////////////////
1419 // Calculate the mean values
1420 ////////////////////////////////
1422 ////////////////////////
1423 Double_t meanAll = 0.0;
1424 Double_t rmsAll = 0.0;
1429 for (Int_t k = 0; k < loop; k++) {
1430 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1431 sector = GetSector(detector);
1433 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1435 rmsAll += value*value;
1441 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1442 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1443 for (Int_t row = 0; row < rowMax; row++) {
1444 for (Int_t col = 0; col < colMax; col++) {
1445 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1447 rmsAll += value*value;
1457 meanAll = meanAll/countAll;
1458 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1460 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1461 /////////////////////////////////////////////////
1463 ////////////////////////////////////////////////
1464 Double_t defaultvalue = -1.0;
1465 if(type==1) defaultvalue = -1.5;
1466 for (Int_t k = 0; k < loop; k++) {
1467 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1468 sector = GetSector(detector);
1469 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1470 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1471 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1473 // remove the results too far away
1474 for (Int_t row = 0; row < rowMax; row++) {
1475 for (Int_t col = 0; col < colMax; col++) {
1476 value = coef[(Int_t)(col*rowMax+row)];
1477 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1478 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1484 //______________________________________________________________________
1485 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1487 // Remove the results too far from the mean and rms
1491 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1493 AliInfo("The Vector Fit is not complete!");
1496 Int_t detector = -1;
1498 Float_t value = 0.0;
1500 /////////////////////////////////
1501 // Calculate the mean values
1502 ////////////////////////////////
1504 ////////////////////////
1505 Double_t meanAll = 0.0;
1506 Double_t rmsAll = 0.0;
1511 for (Int_t k = 0; k < loop; k++) {
1512 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1513 sector = GetSector(detector);
1515 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1518 rmsAll += value*value;
1523 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1524 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1525 for (Int_t row = 0; row < rowMax; row++) {
1526 for (Int_t col = 0; col < colMax; col++) {
1527 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1529 rmsAll += value*value;
1538 meanAll = meanAll/countAll;
1539 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1541 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1542 /////////////////////////////////////////////////
1544 ////////////////////////////////////////////////
1545 for (Int_t k = 0; k < loop; k++) {
1546 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1547 sector = GetSector(detector);
1548 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1549 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1550 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1552 // remove the results too far away
1553 for (Int_t row = 0; row < rowMax; row++) {
1554 for (Int_t col = 0; col < colMax; col++) {
1555 value = coef[(Int_t)(col*rowMax+row)];
1556 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) coef[(Int_t)(col*rowMax+row)] = 100.0;
1561 //______________________________________________________________________
1562 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1564 // ofwhat is equaled to 0: mean value of all passing detectors
1565 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1568 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1570 AliInfo("The Vector Fit is not complete!");
1573 Int_t detector = -1;
1575 Float_t value = 0.0;
1577 /////////////////////////////////
1578 // Calculate the mean values
1579 ////////////////////////////////
1581 ////////////////////////
1582 Double_t meanAll = 0.0;
1583 Double_t meanSupermodule[18];
1584 Double_t meanDetector[540];
1585 Double_t rmsAll = 0.0;
1586 Double_t rmsSupermodule[18];
1587 Double_t rmsDetector[540];
1589 Int_t countSupermodule[18];
1590 Int_t countDetector[540];
1591 for(Int_t sm = 0; sm < 18; sm++){
1592 rmsSupermodule[sm] = 0.0;
1593 meanSupermodule[sm] = 0.0;
1594 countSupermodule[sm] = 0;
1596 for(Int_t det = 0; det < 540; det++){
1597 rmsDetector[det] = 0.0;
1598 meanDetector[det] = 0.0;
1599 countDetector[det] = 0;
1604 for (Int_t k = 0; k < loop; k++) {
1605 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1606 sector = GetSector(detector);
1608 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1610 rmsDetector[detector] += value*value;
1611 meanDetector[detector] += value;
1612 countDetector[detector]++;
1613 rmsSupermodule[sector] += value*value;
1614 meanSupermodule[sector] += value;
1615 countSupermodule[sector]++;
1616 rmsAll += value*value;
1622 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1623 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1624 for (Int_t row = 0; row < rowMax; row++) {
1625 for (Int_t col = 0; col < colMax; col++) {
1626 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1628 rmsDetector[detector] += value*value;
1629 meanDetector[detector] += value;
1630 countDetector[detector]++;
1631 rmsSupermodule[sector] += value*value;
1632 meanSupermodule[sector] += value;
1633 countSupermodule[sector]++;
1634 rmsAll += value*value;
1644 meanAll = meanAll/countAll;
1645 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1647 for(Int_t sm = 0; sm < 18; sm++){
1648 if(countSupermodule[sm] > 0) {
1649 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1650 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1653 for(Int_t det = 0; det < 540; det++){
1654 if(countDetector[det] > 0) {
1655 meanDetector[det] = meanDetector[det]/countDetector[det];
1656 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1659 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1660 ///////////////////////////////////////////////
1661 // Put the mean value for the no-fitted
1662 /////////////////////////////////////////////
1663 for (Int_t k = 0; k < loop; k++) {
1664 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1665 sector = GetSector(detector);
1666 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1667 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1668 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1670 for (Int_t row = 0; row < rowMax; row++) {
1671 for (Int_t col = 0; col < colMax; col++) {
1672 value = coef[(Int_t)(col*rowMax+row)];
1674 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1676 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1677 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1678 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1682 if(fDebugLevel > 1){
1684 if ( !fDebugStreamer ) {
1686 TDirectory *backup = gDirectory;
1687 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1688 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1691 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1693 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1694 "detector="<<detector<<
1706 //______________________________________________________________________
1707 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1709 // ofwhat is equaled to 0: mean value of all passing detectors
1710 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1713 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1715 AliInfo("The Vector Fit is not complete!");
1718 Int_t detector = -1;
1720 Float_t value = 0.0;
1722 /////////////////////////////////
1723 // Calculate the mean values
1724 ////////////////////////////////
1726 ////////////////////////
1727 Double_t meanAll = 0.0;
1728 Double_t rmsAll = 0.0;
1729 Double_t meanSupermodule[18];
1730 Double_t rmsSupermodule[18];
1731 Double_t meanDetector[540];
1732 Double_t rmsDetector[540];
1734 Int_t countSupermodule[18];
1735 Int_t countDetector[540];
1736 for(Int_t sm = 0; sm < 18; sm++){
1737 rmsSupermodule[sm] = 0.0;
1738 meanSupermodule[sm] = 0.0;
1739 countSupermodule[sm] = 0;
1741 for(Int_t det = 0; det < 540; det++){
1742 rmsDetector[det] = 0.0;
1743 meanDetector[det] = 0.0;
1744 countDetector[det] = 0;
1748 for (Int_t k = 0; k < loop; k++) {
1749 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1750 sector = GetSector(detector);
1752 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1754 rmsDetector[detector] += value*value;
1755 meanDetector[detector] += value;
1756 countDetector[detector]++;
1757 rmsSupermodule[sector] += value*value;
1758 meanSupermodule[sector] += value;
1759 countSupermodule[sector]++;
1761 rmsAll += value*value;
1766 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1767 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1768 for (Int_t row = 0; row < rowMax; row++) {
1769 for (Int_t col = 0; col < colMax; col++) {
1770 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1772 rmsDetector[detector] += value*value;
1773 meanDetector[detector] += value;
1774 countDetector[detector]++;
1775 rmsSupermodule[sector] += value*value;
1776 meanSupermodule[sector] += value;
1777 countSupermodule[sector]++;
1778 rmsAll += value*value;
1788 meanAll = meanAll/countAll;
1789 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1791 for(Int_t sm = 0; sm < 18; sm++){
1792 if(countSupermodule[sm] > 0) {
1793 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1794 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1797 for(Int_t det = 0; det < 540; det++){
1798 if(countDetector[det] > 0) {
1799 meanDetector[det] = meanDetector[det]/countDetector[det];
1800 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1803 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1804 ////////////////////////////////////////////
1805 // Put the mean value for the no-fitted
1806 /////////////////////////////////////////////
1807 for (Int_t k = 0; k < loop; k++) {
1808 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1809 sector = GetSector(detector);
1810 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1811 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1812 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1814 for (Int_t row = 0; row < rowMax; row++) {
1815 for (Int_t col = 0; col < colMax; col++) {
1816 value = coef[(Int_t)(col*rowMax+row)];
1818 if((ofwhat == 0) && (meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1820 if((meanDetector[detector] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1821 else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1822 else if((meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1826 if(fDebugLevel > 1){
1828 if ( !fDebugStreamer ) {
1830 TDirectory *backup = gDirectory;
1831 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1832 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1835 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
1837 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1838 "detector="<<detector<<
1851 //_____________________________________________________________________________
1852 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1855 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1856 // It takes the mean value of the coefficients per detector
1857 // This object has to be written in the database
1860 // Create the DetObject
1861 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1863 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1864 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1865 Int_t detector = -1;
1866 Float_t value = 0.0;
1869 for (Int_t k = 0; k < loop; k++) {
1870 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1873 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1877 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1878 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1879 for (Int_t row = 0; row < rowMax; row++) {
1880 for (Int_t col = 0; col < colMax; col++) {
1881 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1882 mean += TMath::Abs(value);
1886 if(count > 0) mean = mean/count;
1888 object->SetValue(detector,mean);
1893 //_____________________________________________________________________________
1894 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1897 // It creates the AliTRDCalDet object from the AliTRDFitInfo
1898 // It takes the mean value of the coefficients per detector
1899 // This object has to be written in the database
1902 // Create the DetObject
1903 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1906 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1907 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1908 Int_t detector = -1;
1909 Float_t value = 0.0;
1911 for (Int_t k = 0; k < loop; k++) {
1912 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1915 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1916 if(!meanOtherBefore){
1917 if(value > 0) value = value*scaleFitFactor;
1919 else value = value*scaleFitFactor;
1920 mean = TMath::Abs(value);
1924 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1925 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1926 for (Int_t row = 0; row < rowMax; row++) {
1927 for (Int_t col = 0; col < colMax; col++) {
1928 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1929 if(!meanOtherBefore) {
1930 if(value > 0) value = value*scaleFitFactor;
1932 else value = value*scaleFitFactor;
1933 mean += TMath::Abs(value);
1937 if(count > 0) mean = mean/count;
1939 object->SetValue(detector,mean);
1944 //_____________________________________________________________________________
1945 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1948 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1949 // It takes the min value of the coefficients per detector
1950 // This object has to be written in the database
1953 // Create the DetObject
1954 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1956 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1957 if(loop != 540) AliInfo("The Vector Fit is not complete!");
1958 Int_t detector = -1;
1959 Float_t value = 0.0;
1961 for (Int_t k = 0; k < loop; k++) {
1962 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1963 Float_t min = 100.0;
1965 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1967 if(value > 70.0) value = value-100.0;
1972 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1973 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1974 for (Int_t row = 0; row < rowMax; row++) {
1975 for (Int_t col = 0; col < colMax; col++) {
1976 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1978 if(value > 70.0) value = value-100.0;
1980 if(min > value) min = value;
1984 object->SetValue(detector,min);
1990 //_____________________________________________________________________________
1991 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
1994 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1995 // It takes the min value of the coefficients per detector
1996 // This object has to be written in the database
1999 // Create the DetObject
2000 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2003 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2004 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2005 Int_t detector = -1;
2006 Float_t value = 0.0;
2008 for (Int_t k = 0; k < loop; k++) {
2009 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2011 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2012 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2013 Float_t min = 100.0;
2014 for (Int_t row = 0; row < rowMax; row++) {
2015 for (Int_t col = 0; col < colMax; col++) {
2016 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2017 mean += -TMath::Abs(value);
2021 if(count > 0) mean = mean/count;
2023 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2024 object->SetValue(detector,-TMath::Abs(value));
2030 //_____________________________________________________________________________
2031 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2034 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2035 // You need first to create the object for the detectors,
2036 // where the mean value is put.
2037 // This object has to be written in the database
2040 // Create the DetObject
2041 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2044 for(Int_t k = 0; k < 540; k++){
2045 AliTRDCalROC *calROC = object->GetCalROC(k);
2046 Int_t nchannels = calROC->GetNchannels();
2047 for(Int_t ch = 0; ch < nchannels; ch++){
2048 calROC->SetValue(ch,1.0);
2054 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2055 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2056 Int_t detector = -1;
2057 Float_t value = 0.0;
2059 for (Int_t k = 0; k < loop; k++) {
2060 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2061 AliTRDCalROC *calROC = object->GetCalROC(detector);
2062 Float_t mean = detobject->GetValue(detector);
2063 if(mean == 0) continue;
2064 Int_t rowMax = calROC->GetNrows();
2065 Int_t colMax = calROC->GetNcols();
2066 for (Int_t row = 0; row < rowMax; row++) {
2067 for (Int_t col = 0; col < colMax; col++) {
2068 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2069 if(value > 0) value = value*scaleFitFactor;
2070 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2078 //_____________________________________________________________________________
2079 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2082 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2083 // You need first to create the object for the detectors,
2084 // where the mean value is put.
2085 // This object has to be written in the database
2088 // Create the DetObject
2089 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2092 for(Int_t k = 0; k < 540; k++){
2093 AliTRDCalROC *calROC = object->GetCalROC(k);
2094 Int_t nchannels = calROC->GetNchannels();
2095 for(Int_t ch = 0; ch < nchannels; ch++){
2096 calROC->SetValue(ch,1.0);
2102 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2103 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2104 Int_t detector = -1;
2105 Float_t value = 0.0;
2107 for (Int_t k = 0; k < loop; k++) {
2108 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2109 AliTRDCalROC *calROC = object->GetCalROC(detector);
2110 Float_t mean = detobject->GetValue(detector);
2111 if(mean == 0) continue;
2112 Int_t rowMax = calROC->GetNrows();
2113 Int_t colMax = calROC->GetNcols();
2114 for (Int_t row = 0; row < rowMax; row++) {
2115 for (Int_t col = 0; col < colMax; col++) {
2116 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2117 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2125 //_____________________________________________________________________________
2126 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2129 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2130 // You need first to create the object for the detectors,
2131 // where the mean value is put.
2132 // This object has to be written in the database
2135 // Create the DetObject
2136 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2139 for(Int_t k = 0; k < 540; k++){
2140 AliTRDCalROC *calROC = object->GetCalROC(k);
2141 Int_t nchannels = calROC->GetNchannels();
2142 for(Int_t ch = 0; ch < nchannels; ch++){
2143 calROC->SetValue(ch,0.0);
2149 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2150 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2151 Int_t detector = -1;
2152 Float_t value = 0.0;
2154 for (Int_t k = 0; k < loop; k++) {
2155 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2156 AliTRDCalROC *calROC = object->GetCalROC(detector);
2157 Float_t min = detobject->GetValue(detector);
2158 Int_t rowMax = calROC->GetNrows();
2159 Int_t colMax = calROC->GetNcols();
2160 for (Int_t row = 0; row < rowMax; row++) {
2161 for (Int_t col = 0; col < colMax; col++) {
2162 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2164 if(value > 70.0) value = value - 100.0;
2166 calROC->SetValue(col,row,value-min);
2174 //_____________________________________________________________________________
2175 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2178 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2179 // This object has to be written in the database
2182 // Create the DetObject
2183 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2185 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2186 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2187 Int_t detector = -1;
2188 Float_t value = 0.0;
2190 for (Int_t k = 0; k < loop; k++) {
2191 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2192 AliTRDCalROC *calROC = object->GetCalROC(detector);
2193 Int_t rowMax = calROC->GetNrows();
2194 Int_t colMax = calROC->GetNcols();
2195 for (Int_t row = 0; row < rowMax; row++) {
2196 for (Int_t col = 0; col < colMax; col++) {
2197 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2198 calROC->SetValue(col,row,TMath::Abs(value));
2206 //_____________________________________________________________________________
2207 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2210 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2211 // 0 successful fit 1 not successful fit
2212 // mean is the mean value over the successful fit
2213 // do not use it for t0: no meaning
2216 // Create the CalObject
2217 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2221 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2223 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2224 for(Int_t k = 0; k < 540; k++){
2225 object->SetValue(k,1.0);
2228 Int_t detector = -1;
2229 Float_t value = 0.0;
2231 for (Int_t k = 0; k < loop; k++) {
2232 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2233 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2234 if(value <= 0) object->SetValue(detector,1.0);
2236 object->SetValue(detector,0.0);
2241 if(count > 0) mean /= count;
2244 //_____________________________________________________________________________
2245 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2248 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2249 // 0 not successful fit 1 successful fit
2250 // mean mean value over the successful fit
2253 // Create the CalObject
2254 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2258 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2260 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2261 for(Int_t k = 0; k < 540; k++){
2262 AliTRDCalROC *calROC = object->GetCalROC(k);
2263 Int_t nchannels = calROC->GetNchannels();
2264 for(Int_t ch = 0; ch < nchannels; ch++){
2265 calROC->SetValue(ch,1.0);
2269 Int_t detector = -1;
2270 Float_t value = 0.0;
2272 for (Int_t k = 0; k < loop; k++) {
2273 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2274 AliTRDCalROC *calROC = object->GetCalROC(detector);
2275 Int_t nchannels = calROC->GetNchannels();
2276 for (Int_t ch = 0; ch < nchannels; ch++) {
2277 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2278 if(value <= 0) calROC->SetValue(ch,1.0);
2280 calROC->SetValue(ch,0.0);
2286 if(count > 0) mean /= count;
2289 //_____________________________________________________________________________
2290 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2293 // Set FitPH if 1 then each detector will be fitted
2296 if (periodeFitPH > 0) {
2297 fFitPHPeriode = periodeFitPH;
2300 AliInfo("periodeFitPH must be higher than 0!");
2304 //_____________________________________________________________________________
2305 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2308 // The fit of the deposited charge distribution begins at
2309 // histo->Mean()/beginFitCharge
2310 // You can here set beginFitCharge
2313 if (beginFitCharge > 0) {
2314 fBeginFitCharge = beginFitCharge;
2317 AliInfo("beginFitCharge must be strict positif!");
2322 //_____________________________________________________________________________
2323 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2326 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2327 // You can here set t0Shift0
2331 fT0Shift0 = t0Shift;
2334 AliInfo("t0Shift0 must be strict positif!");
2339 //_____________________________________________________________________________
2340 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2343 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2344 // You can here set t0Shift1
2348 fT0Shift1 = t0Shift;
2351 AliInfo("t0Shift must be strict positif!");
2356 //_____________________________________________________________________________
2357 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2360 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2361 // You can here set rangeFitPRF
2364 if ((rangeFitPRF > 0) &&
2365 (rangeFitPRF <= 1.5)) {
2366 fRangeFitPRF = rangeFitPRF;
2369 AliInfo("rangeFitPRF must be between 0 and 1.0");
2374 //_____________________________________________________________________________
2375 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2378 // Minimum entries for fitting
2381 if (minEntries > 0) {
2382 fMinEntries = minEntries;
2385 AliInfo("fMinEntries must be >= 0.");
2390 //_____________________________________________________________________________
2391 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2394 // Rebin with rebin time less bins the Ch histo
2395 // You can set here rebin that should divide the number of bins of CH histo
2400 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2403 AliInfo("You have to choose a positiv value!");
2407 //_____________________________________________________________________________
2408 Bool_t AliTRDCalibraFit::FillVectorFit()
2411 // For the Fit functions fill the vector Fit
2414 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2417 if (GetStack(fCountDet) == 2) {
2424 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2425 Float_t *coef = new Float_t[ntotal];
2426 for (Int_t i = 0; i < ntotal; i++) {
2427 coef[i] = fCurrentCoefDetector[i];
2430 Int_t detector = fCountDet;
2432 fitInfo->SetCoef(coef);
2433 fitInfo->SetDetector(detector);
2434 fVectorFit.Add((TObject *) fitInfo);
2439 //_____________________________________________________________________________
2440 Bool_t AliTRDCalibraFit::FillVectorFit2()
2443 // For the Fit functions fill the vector Fit
2446 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2449 if (GetStack(fCountDet) == 2) {
2456 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2457 Float_t *coef = new Float_t[ntotal];
2458 for (Int_t i = 0; i < ntotal; i++) {
2459 coef[i] = fCurrentCoefDetector2[i];
2462 Int_t detector = fCountDet;
2464 fitInfo->SetCoef(coef);
2465 fitInfo->SetDetector(detector);
2466 fVectorFit2.Add((TObject *) fitInfo);
2471 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2472 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2475 // Init the number of expected bins and fDect1[i] fDect2[i]
2478 gStyle->SetPalette(1);
2479 gStyle->SetOptStat(1111);
2480 gStyle->SetPadBorderMode(0);
2481 gStyle->SetCanvasColor(10);
2482 gStyle->SetPadLeftMargin(0.13);
2483 gStyle->SetPadRightMargin(0.01);
2485 // Mode groups of pads: the total number of bins!
2486 CalculNumberOfBinsExpected(i);
2488 // Quick verification that we have the good pad calibration mode!
2489 if (fNumberOfBinsExpected != nbins) {
2490 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2494 // Security for fDebug 3 and 4
2495 if ((fDebugLevel >= 3) &&
2499 AliInfo("This detector doesn't exit!");
2503 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2504 CalculDect1Dect2(i);
2509 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2510 Bool_t AliTRDCalibraFit::InitFitCH()
2513 // Init the fVectorFitCH for normalisation
2514 // Init the histo for debugging
2519 fScaleFitFactor = 0.0;
2520 fCurrentCoefDetector = new Float_t[2304];
2521 for (Int_t k = 0; k < 2304; k++) {
2522 fCurrentCoefDetector[k] = 0.0;
2524 fVectorFit.SetName("gainfactorscoefficients");
2526 // fDebug == 0 nothing
2527 // fDebug == 1 and fFitVoir no histo
2528 if (fDebugLevel == 1) {
2529 if(!CheckFitVoir()) return kFALSE;
2531 //Get the CalDet object
2533 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2535 AliInfo("Could not get calibDB");
2538 if(fCalDet) delete fCalDet;
2539 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2542 Float_t devalue = 1.0;
2543 if(fCalDet) delete fCalDet;
2544 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2545 for(Int_t k = 0; k < 540; k++){
2546 fCalDet->SetValue(k,devalue);
2552 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2553 Bool_t AliTRDCalibraFit::InitFitPH()
2556 // Init the arrays of results
2557 // Init the histos for debugging
2561 fVectorFit.SetName("driftvelocitycoefficients");
2562 fVectorFit2.SetName("t0coefficients");
2564 fCurrentCoefDetector = new Float_t[2304];
2565 for (Int_t k = 0; k < 2304; k++) {
2566 fCurrentCoefDetector[k] = 0.0;
2569 fCurrentCoefDetector2 = new Float_t[2304];
2570 for (Int_t k = 0; k < 2304; k++) {
2571 fCurrentCoefDetector2[k] = 0.0;
2574 //fDebug == 0 nothing
2575 // fDebug == 1 and fFitVoir no histo
2576 if (fDebugLevel == 1) {
2577 if(!CheckFitVoir()) return kFALSE;
2579 //Get the CalDet object
2581 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2583 AliInfo("Could not get calibDB");
2586 if(fCalDet) delete fCalDet;
2587 if(fCalDet2) delete fCalDet2;
2588 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2589 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2592 Float_t devalue = 1.5;
2593 Float_t devalue2 = 0.0;
2594 if(fCalDet) delete fCalDet;
2595 if(fCalDet2) delete fCalDet2;
2596 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2597 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2598 for(Int_t k = 0; k < 540; k++){
2599 fCalDet->SetValue(k,devalue);
2600 fCalDet2->SetValue(k,devalue2);
2605 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2606 Bool_t AliTRDCalibraFit::InitFitPRF()
2609 // Init the calibration mode (Nz, Nrphi), the histograms for
2610 // debugging the fit methods if fDebug > 0,
2614 fVectorFit.SetName("prfwidthcoefficients");
2616 fCurrentCoefDetector = new Float_t[2304];
2617 for (Int_t k = 0; k < 2304; k++) {
2618 fCurrentCoefDetector[k] = 0.0;
2621 // fDebug == 0 nothing
2622 // fDebug == 1 and fFitVoir no histo
2623 if (fDebugLevel == 1) {
2624 if(!CheckFitVoir()) return kFALSE;
2628 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2629 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2632 // Init the fCalDet, fVectorFit fCurrentCoefDetector
2637 fCurrentCoefDetector = new Float_t[2304];
2638 fCurrentCoefDetector2 = new Float_t[2304];
2639 for (Int_t k = 0; k < 2304; k++) {
2640 fCurrentCoefDetector[k] = 0.0;
2641 fCurrentCoefDetector2[k] = 0.0;
2644 //printf("test0\n");
2646 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2648 AliInfo("Could not get calibDB");
2652 //Get the CalDet object
2654 if(fCalDet) delete fCalDet;
2655 if(fCalDet2) delete fCalDet2;
2656 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2657 //printf("test1\n");
2658 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2659 //printf("test2\n");
2660 for(Int_t k = 0; k < 540; k++){
2661 fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2663 //printf("test3\n");
2666 Float_t devalue = 1.5;
2667 Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2668 if(fCalDet) delete fCalDet;
2669 if(fCalDet2) delete fCalDet2;
2670 //printf("test1\n");
2671 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2672 fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2673 //printf("test2\n");
2674 for(Int_t k = 0; k < 540; k++){
2675 fCalDet->SetValue(k,devalue);
2676 fCalDet2->SetValue(k,devalue2);
2678 //printf("test3\n");
2683 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2684 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2687 // Init the current detector where we are fCountDet and the
2688 // next fCount for the functions Fit...
2691 // Loop on the Xbins of ch!!
2692 fCountDet = -1; // Current detector
2693 fCount = 0; // To find the next detector
2696 if (fDebugLevel >= 3) {
2697 // Set countdet to the detector
2698 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2699 // Set counter to write at the end of the detector
2701 // Get the right calib objects
2704 if(fDebugLevel == 1) {
2706 fCalibraMode->CalculXBins(fCountDet,i);
2707 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2708 while(fCalibraMode->GetXbins(i) <=fFitVoir){
2710 fCalibraMode->CalculXBins(fCountDet,i);
2711 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2717 fCount = fCalibraMode->GetXbins(i);
2719 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2720 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2721 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2722 ,(Int_t) GetStack(fCountDet)
2723 ,(Int_t) GetSector(fCountDet),i);
2726 //_______________________________________________________________________________
2727 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2730 // Calculate the number of bins expected (calibration groups)
2733 fNumberOfBinsExpected = 0;
2735 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2736 fNumberOfBinsExpected = 1;
2740 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2741 fNumberOfBinsExpected = 18;
2745 fCalibraMode->ModePadCalibration(2,i);
2746 fCalibraMode->ModePadFragmentation(0,2,0,i);
2747 fCalibraMode->SetDetChamb2(i);
2748 if (fDebugLevel > 1) {
2749 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2751 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2752 fCalibraMode->ModePadCalibration(0,i);
2753 fCalibraMode->ModePadFragmentation(0,0,0,i);
2754 fCalibraMode->SetDetChamb0(i);
2755 if (fDebugLevel > 1) {
2756 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2758 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2761 //_______________________________________________________________________________
2762 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2765 // Calculate the range of fits
2770 if (fDebugLevel == 1) {
2774 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2776 fDect2 = fNumberOfBinsExpected;
2778 if (fDebugLevel >= 3) {
2779 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2780 fCalibraMode->CalculXBins(fCountDet,i);
2781 fDect1 = fCalibraMode->GetXbins(i);
2782 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2783 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2784 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2785 ,(Int_t) GetStack(fCountDet)
2786 ,(Int_t) GetSector(fCountDet),i);
2787 // Set for the next detector
2788 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2791 //_______________________________________________________________________________
2792 Bool_t AliTRDCalibraFit::CheckFitVoir()
2795 // Check if fFitVoir is in the range
2798 if (fFitVoir < fNumberOfBinsExpected) {
2799 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2802 AliInfo("fFitVoir is out of range of the histo!");
2807 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2808 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2811 // See if we are in a new detector and update the
2812 // variables fNfragZ and fNfragRphi if yes
2813 // Will never happen for only one detector (3 and 4)
2814 // Doesn't matter for 2
2816 if (fCount == idect) {
2817 // On en est au detector (or first detector in the group)
2819 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2820 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2821 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2822 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2823 ,(Int_t) GetStack(fCountDet)
2824 ,(Int_t) GetSector(fCountDet),i);
2825 // Set for the next detector
2826 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2831 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2832 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2835 // Reconstruct the min pad row, max pad row, min pad col and
2836 // max pad col of the calibration group for the Fit functions
2837 // idect is the calibration group inside the detector
2839 if (fDebugLevel != 1) {
2840 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2842 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2843 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2845 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2846 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2849 // For the case where there are not enough entries in the histograms
2850 // of the calibration group, the value present in the choosen database
2851 // will be put. A negativ sign enables to know that a fit was not possible.
2854 if (fDebugLevel == 1) {
2855 AliInfo("The element has not enough statistic to be fitted");
2857 else if (fNbDet > 0){
2858 Int_t firstdetector = fCountDet;
2859 Int_t lastdetector = fCountDet+fNbDet;
2860 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2861 ,idect,firstdetector,lastdetector));
2862 // loop over detectors
2863 for(Int_t det = firstdetector; det < lastdetector; det++){
2865 //Set the calibration object again
2869 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2871 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2872 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2873 ,(Int_t) GetStack(fCountDet)
2874 ,(Int_t) GetSector(fCountDet),0);
2875 // Reconstruct row min row max
2876 ReconstructFitRowMinRowMax(idect,0);
2878 // Calcul the coef from the database choosen for the detector
2879 CalculChargeCoefMean(kFALSE);
2881 //stack 2, not stack 2
2883 if(GetStack(fCountDet) == 2) factor = 12;
2886 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2887 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2888 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2889 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2893 //Put default value negative
2894 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2895 fCurrentCoefE = 0.0;
2900 if(fDebugLevel > 1){
2902 if ( !fDebugStreamer ) {
2904 TDirectory *backup = gDirectory;
2905 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2906 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2909 Int_t detector = fCountDet;
2910 Int_t caligroup = idect;
2911 Short_t rowmin = fCalibraMode->GetRowMin(0);
2912 Short_t rowmax = fCalibraMode->GetRowMax(0);
2913 Short_t colmin = fCalibraMode->GetColMin(0);
2914 Short_t colmax = fCalibraMode->GetColMax(0);
2915 Float_t gf = fCurrentCoef[0];
2916 Float_t gfs = fCurrentCoef[1];
2917 Float_t gfE = fCurrentCoefE;
2919 (*fDebugStreamer) << "FillFillCH" <<
2920 "detector=" << detector <<
2921 "caligroup=" << caligroup <<
2922 "rowmin=" << rowmin <<
2923 "rowmax=" << rowmax <<
2924 "colmin=" << colmin <<
2925 "colmax=" << colmax <<
2933 for (Int_t k = 0; k < 2304; k++) {
2934 fCurrentCoefDetector[k] = 0.0;
2938 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2942 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2943 ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2945 // Calcul the coef from the database choosen
2946 CalculChargeCoefMean(kFALSE);
2948 //stack 2, not stack 2
2950 if(GetStack(fCountDet) == 2) factor = 12;
2953 // Fill the fCurrentCoefDetector with negative value to say: not fitted
2954 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2955 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2956 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2960 //Put default value negative
2961 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2962 fCurrentCoefE = 0.0;
2971 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2972 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2975 // For the case where there are not enough entries in the histograms
2976 // of the calibration group, the value present in the choosen database
2977 // will be put. A negativ sign enables to know that a fit was not possible.
2979 if (fDebugLevel == 1) {
2980 AliInfo("The element has not enough statistic to be fitted");
2982 else if (fNbDet > 0) {
2984 Int_t firstdetector = fCountDet;
2985 Int_t lastdetector = fCountDet+fNbDet;
2986 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2987 ,idect,firstdetector,lastdetector));
2988 // loop over detectors
2989 for(Int_t det = firstdetector; det < lastdetector; det++){
2991 //Set the calibration object again
2995 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2997 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
2998 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2999 ,(Int_t) GetStack(fCountDet)
3000 ,(Int_t) GetSector(fCountDet),1);
3001 // Reconstruct row min row max
3002 ReconstructFitRowMinRowMax(idect,1);
3004 // Calcul the coef from the database choosen for the detector
3005 CalculVdriftCoefMean();
3008 //stack 2, not stack 2
3010 if(GetStack(fCountDet) == 2) factor = 12;
3013 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3014 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3015 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3016 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3017 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3021 //Put default value negative
3022 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3023 fCurrentCoefE = 0.0;
3024 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3025 fCurrentCoefE2 = 0.0;
3031 if(fDebugLevel > 1){
3033 if ( !fDebugStreamer ) {
3035 TDirectory *backup = gDirectory;
3036 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3037 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3041 Int_t detector = fCountDet;
3042 Int_t caligroup = idect;
3043 Short_t rowmin = fCalibraMode->GetRowMin(1);
3044 Short_t rowmax = fCalibraMode->GetRowMax(1);
3045 Short_t colmin = fCalibraMode->GetColMin(1);
3046 Short_t colmax = fCalibraMode->GetColMax(1);
3047 Float_t vf = fCurrentCoef[0];
3048 Float_t vs = fCurrentCoef[1];
3049 Float_t vfE = fCurrentCoefE;
3050 Float_t t0f = fCurrentCoef2[0];
3051 Float_t t0s = fCurrentCoef2[1];
3052 Float_t t0E = fCurrentCoefE2;
3056 (* fDebugStreamer) << "FillFillPH"<<
3057 "detector="<<detector<<
3058 "nentries="<<nentries<<
3059 "caligroup="<<caligroup<<
3073 for (Int_t k = 0; k < 2304; k++) {
3074 fCurrentCoefDetector[k] = 0.0;
3075 fCurrentCoefDetector2[k] = 0.0;
3079 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3083 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3084 ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3086 CalculVdriftCoefMean();
3089 //stack 2 and not stack 2
3091 if(GetStack(fCountDet) == 2) factor = 12;
3095 // Fill the fCurrentCoefDetector 2
3096 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3097 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3098 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3099 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3103 // Put the default value
3104 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3105 fCurrentCoefE = 0.0;
3106 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3107 fCurrentCoefE2 = 0.0;
3109 FillFillPH(idect,nentries);
3118 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3119 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3122 // For the case where there are not enough entries in the histograms
3123 // of the calibration group, the value present in the choosen database
3124 // will be put. A negativ sign enables to know that a fit was not possible.
3127 if (fDebugLevel == 1) {
3128 AliInfo("The element has not enough statistic to be fitted");
3130 else if (fNbDet > 0){
3132 Int_t firstdetector = fCountDet;
3133 Int_t lastdetector = fCountDet+fNbDet;
3134 AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3135 ,idect,firstdetector,lastdetector));
3137 // loop over detectors
3138 for(Int_t det = firstdetector; det < lastdetector; det++){
3140 //Set the calibration object again
3144 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3146 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3147 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3148 ,(Int_t) GetStack(fCountDet)
3149 ,(Int_t) GetSector(fCountDet),2);
3150 // Reconstruct row min row max
3151 ReconstructFitRowMinRowMax(idect,2);
3153 // Calcul the coef from the database choosen for the detector
3154 CalculPRFCoefMean();
3156 //stack 2, not stack 2
3158 if(GetStack(fCountDet) == 2) factor = 12;
3161 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3162 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3163 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3164 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3168 //Put default value negative
3169 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3170 fCurrentCoefE = 0.0;
3175 if(fDebugLevel > 1){
3177 if ( !fDebugStreamer ) {
3179 TDirectory *backup = gDirectory;
3180 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3181 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3184 Int_t detector = fCountDet;
3185 Int_t layer = GetLayer(fCountDet);
3186 Int_t caligroup = idect;
3187 Short_t rowmin = fCalibraMode->GetRowMin(2);
3188 Short_t rowmax = fCalibraMode->GetRowMax(2);
3189 Short_t colmin = fCalibraMode->GetColMin(2);
3190 Short_t colmax = fCalibraMode->GetColMax(2);
3191 Float_t widf = fCurrentCoef[0];
3192 Float_t wids = fCurrentCoef[1];
3193 Float_t widfE = fCurrentCoefE;
3195 (* fDebugStreamer) << "FillFillPRF"<<
3196 "detector="<<detector<<
3198 "caligroup="<<caligroup<<
3209 for (Int_t k = 0; k < 2304; k++) {
3210 fCurrentCoefDetector[k] = 0.0;
3214 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3218 AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3219 ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3221 CalculPRFCoefMean();
3223 // stack 2 and not stack 2
3225 if(GetStack(fCountDet) == 2) factor = 12;
3229 // Fill the fCurrentCoefDetector
3230 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3231 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3232 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3236 // Put the default value
3237 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3238 fCurrentCoefE = 0.0;
3246 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3247 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3250 // For the case where there are not enough entries in the histograms
3251 // of the calibration group, the value present in the choosen database
3252 // will be put. A negativ sign enables to know that a fit was not possible.
3255 // Calcul the coef from the database choosen
3256 CalculVdriftLorentzCoef();
3259 if(GetStack(fCountDet) == 2) factor = 1728;
3263 // Fill the fCurrentCoefDetector
3264 for (Int_t k = 0; k < factor; k++) {
3265 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3266 // should be negative
3267 fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3271 //Put default opposite sign
3272 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3273 fCurrentCoefE = 0.0;
3274 fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3275 fCurrentCoefE2 = 0.0;
3277 FillFillLinearFitter();
3282 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3283 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3286 // Fill the coefficients found with the fits or other
3287 // methods from the Fit functions
3290 if (fDebugLevel != 1) {
3292 Int_t firstdetector = fCountDet;
3293 Int_t lastdetector = fCountDet+fNbDet;
3294 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3295 ,idect,firstdetector,lastdetector));
3296 // loop over detectors
3297 for(Int_t det = firstdetector; det < lastdetector; det++){
3299 //Set the calibration object again
3303 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3305 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3306 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3307 ,(Int_t) GetStack(fCountDet)
3308 ,(Int_t) GetSector(fCountDet),0);
3309 // Reconstruct row min row max
3310 ReconstructFitRowMinRowMax(idect,0);
3312 // Calcul the coef from the database choosen for the detector
3313 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3314 else CalculChargeCoefMean(kTRUE);
3316 //stack 2, not stack 2
3318 if(GetStack(fCountDet) == 2) factor = 12;
3321 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3322 Double_t coeftoput = 1.0;
3323 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3324 else coeftoput = fCurrentCoef[0];
3325 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3326 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3327 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3334 if(fDebugLevel > 1){
3336 if ( !fDebugStreamer ) {
3338 TDirectory *backup = gDirectory;
3339 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3340 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3343 Int_t detector = fCountDet;
3344 Int_t caligroup = idect;
3345 Short_t rowmin = fCalibraMode->GetRowMin(0);
3346 Short_t rowmax = fCalibraMode->GetRowMax(0);
3347 Short_t colmin = fCalibraMode->GetColMin(0);
3348 Short_t colmax = fCalibraMode->GetColMax(0);
3349 Float_t gf = fCurrentCoef[0];
3350 Float_t gfs = fCurrentCoef[1];
3351 Float_t gfE = fCurrentCoefE;
3353 (*fDebugStreamer) << "FillFillCH" <<
3354 "detector=" << detector <<
3355 "caligroup=" << caligroup <<
3356 "rowmin=" << rowmin <<
3357 "rowmax=" << rowmax <<
3358 "colmin=" << colmin <<
3359 "colmax=" << colmax <<
3367 for (Int_t k = 0; k < 2304; k++) {
3368 fCurrentCoefDetector[k] = 0.0;
3372 //printf("Check the count now: fCountDet %d\n",fCountDet);
3377 if(GetStack(fCountDet) == 2) factor = 12;
3380 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3381 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3382 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3393 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3394 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3397 // Fill the coefficients found with the fits or other
3398 // methods from the Fit functions
3401 if (fDebugLevel != 1) {
3404 Int_t firstdetector = fCountDet;
3405 Int_t lastdetector = fCountDet+fNbDet;
3406 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3407 ,idect,firstdetector,lastdetector));
3409 // loop over detectors
3410 for(Int_t det = firstdetector; det < lastdetector; det++){
3412 //Set the calibration object again
3416 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3418 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3419 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3420 ,(Int_t) GetStack(fCountDet)
3421 ,(Int_t) GetSector(fCountDet),1);
3422 // Reconstruct row min row max
3423 ReconstructFitRowMinRowMax(idect,1);
3425 // Calcul the coef from the database choosen for the detector
3426 CalculVdriftCoefMean();
3429 //stack 2, not stack 2
3431 if(GetStack(fCountDet) == 2) factor = 12;
3434 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3435 Double_t coeftoput = 1.5;
3436 Double_t coeftoput2 = 0.0;
3438 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3439 else coeftoput = fCurrentCoef[0];
3441 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3442 else coeftoput2 = fCurrentCoef2[0];
3444 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3445 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3446 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3447 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3455 if(fDebugLevel > 1){
3457 if ( !fDebugStreamer ) {
3459 TDirectory *backup = gDirectory;
3460 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3461 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3465 Int_t detector = fCountDet;
3466 Int_t caligroup = idect;
3467 Short_t rowmin = fCalibraMode->GetRowMin(1);
3468 Short_t rowmax = fCalibraMode->GetRowMax(1);
3469 Short_t colmin = fCalibraMode->GetColMin(1);
3470 Short_t colmax = fCalibraMode->GetColMax(1);
3471 Float_t vf = fCurrentCoef[0];
3472 Float_t vs = fCurrentCoef[1];
3473 Float_t vfE = fCurrentCoefE;
3474 Float_t t0f = fCurrentCoef2[0];
3475 Float_t t0s = fCurrentCoef2[1];
3476 Float_t t0E = fCurrentCoefE2;
3480 (* fDebugStreamer) << "FillFillPH"<<
3481 "detector="<<detector<<
3482 "nentries="<<nentries<<
3483 "caligroup="<<caligroup<<
3497 for (Int_t k = 0; k < 2304; k++) {
3498 fCurrentCoefDetector[k] = 0.0;
3499 fCurrentCoefDetector2[k] = 0.0;
3503 //printf("Check the count now: fCountDet %d\n",fCountDet);
3508 if(GetStack(fCountDet) == 2) factor = 12;
3511 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3512 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3513 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3514 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3518 FillFillPH(idect,nentries);
3523 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3524 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3527 // Fill the coefficients found with the fits or other
3528 // methods from the Fit functions
3531 if (fDebugLevel != 1) {
3534 Int_t firstdetector = fCountDet;
3535 Int_t lastdetector = fCountDet+fNbDet;
3536 AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3537 ,idect,firstdetector,lastdetector));
3539 // loop over detectors
3540 for(Int_t det = firstdetector; det < lastdetector; det++){
3542 //Set the calibration object again
3546 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3548 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3549 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3550 ,(Int_t) GetStack(fCountDet)
3551 ,(Int_t) GetSector(fCountDet),2);
3552 // Reconstruct row min row max
3553 ReconstructFitRowMinRowMax(idect,2);
3555 // Calcul the coef from the database choosen for the detector
3556 CalculPRFCoefMean();
3558 //stack 2, not stack 2
3560 if(GetStack(fCountDet) == 2) factor = 12;
3563 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3564 Double_t coeftoput = 1.0;
3565 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3566 else coeftoput = fCurrentCoef[0];
3567 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3568 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3569 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3576 if(fDebugLevel > 1){
3578 if ( !fDebugStreamer ) {
3580 TDirectory *backup = gDirectory;
3581 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3582 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3585 Int_t detector = fCountDet;
3586 Int_t layer = GetLayer(fCountDet);
3587 Int_t caligroup = idect;
3588 Short_t rowmin = fCalibraMode->GetRowMin(2);
3589 Short_t rowmax = fCalibraMode->GetRowMax(2);
3590 Short_t colmin = fCalibraMode->GetColMin(2);
3591 Short_t colmax = fCalibraMode->GetColMax(2);
3592 Float_t widf = fCurrentCoef[0];
3593 Float_t wids = fCurrentCoef[1];
3594 Float_t widfE = fCurrentCoefE;
3596 (* fDebugStreamer) << "FillFillPRF"<<
3597 "detector="<<detector<<
3599 "caligroup="<<caligroup<<
3610 for (Int_t k = 0; k < 2304; k++) {
3611 fCurrentCoefDetector[k] = 0.0;
3615 //printf("Check the count now: fCountDet %d\n",fCountDet);
3620 if(GetStack(fCountDet) == 2) factor = 12;
3623 // Pointer to the branch
3624 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3625 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3626 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3636 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3637 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3640 // Fill the coefficients found with the fits or other
3641 // methods from the Fit functions
3645 if(GetStack(fCountDet) == 2) factor = 1728;
3648 // Pointer to the branch
3649 for (Int_t k = 0; k < factor; k++) {
3650 fCurrentCoefDetector[k] = fCurrentCoef[0];
3651 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3654 FillFillLinearFitter();
3659 //________________________________________________________________________________
3660 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3663 // DebugStream and fVectorFit
3666 // End of one detector
3667 if ((idect == (fCount-1))) {
3670 for (Int_t k = 0; k < 2304; k++) {
3671 fCurrentCoefDetector[k] = 0.0;
3675 if(fDebugLevel > 1){
3677 if ( !fDebugStreamer ) {
3679 TDirectory *backup = gDirectory;
3680 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3681 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3684 Int_t detector = fCountDet;
3685 Int_t caligroup = idect;
3686 Short_t rowmin = fCalibraMode->GetRowMin(0);
3687 Short_t rowmax = fCalibraMode->GetRowMax(0);
3688 Short_t colmin = fCalibraMode->GetColMin(0);
3689 Short_t colmax = fCalibraMode->GetColMax(0);
3690 Float_t gf = fCurrentCoef[0];
3691 Float_t gfs = fCurrentCoef[1];
3692 Float_t gfE = fCurrentCoefE;
3694 (*fDebugStreamer) << "FillFillCH" <<
3695 "detector=" << detector <<
3696 "caligroup=" << caligroup <<
3697 "rowmin=" << rowmin <<
3698 "rowmax=" << rowmax <<
3699 "colmin=" << colmin <<
3700 "colmax=" << colmax <<
3708 //________________________________________________________________________________
3709 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3712 // DebugStream and fVectorFit and fVectorFit2
3715 // End of one detector
3716 if ((idect == (fCount-1))) {
3720 for (Int_t k = 0; k < 2304; k++) {
3721 fCurrentCoefDetector[k] = 0.0;
3722 fCurrentCoefDetector2[k] = 0.0;
3726 if(fDebugLevel > 1){
3728 if ( !fDebugStreamer ) {
3730 TDirectory *backup = gDirectory;
3731 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3732 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3736 Int_t detector = fCountDet;
3737 Int_t caligroup = idect;
3738 Short_t rowmin = fCalibraMode->GetRowMin(1);
3739 Short_t rowmax = fCalibraMode->GetRowMax(1);
3740 Short_t colmin = fCalibraMode->GetColMin(1);
3741 Short_t colmax = fCalibraMode->GetColMax(1);
3742 Float_t vf = fCurrentCoef[0];
3743 Float_t vs = fCurrentCoef[1];
3744 Float_t vfE = fCurrentCoefE;
3745 Float_t t0f = fCurrentCoef2[0];
3746 Float_t t0s = fCurrentCoef2[1];
3747 Float_t t0E = fCurrentCoefE2;
3751 (* fDebugStreamer) << "FillFillPH"<<
3752 "detector="<<detector<<
3753 "nentries="<<nentries<<
3754 "caligroup="<<caligroup<<
3769 //________________________________________________________________________________
3770 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3773 // DebugStream and fVectorFit
3776 // End of one detector
3777 if ((idect == (fCount-1))) {
3780 for (Int_t k = 0; k < 2304; k++) {
3781 fCurrentCoefDetector[k] = 0.0;
3786 if(fDebugLevel > 1){
3788 if ( !fDebugStreamer ) {
3790 TDirectory *backup = gDirectory;
3791 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3792 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3795 Int_t detector = fCountDet;
3796 Int_t layer = GetLayer(fCountDet);
3797 Int_t caligroup = idect;
3798 Short_t rowmin = fCalibraMode->GetRowMin(2);
3799 Short_t rowmax = fCalibraMode->GetRowMax(2);
3800 Short_t colmin = fCalibraMode->GetColMin(2);
3801 Short_t colmax = fCalibraMode->GetColMax(2);
3802 Float_t widf = fCurrentCoef[0];
3803 Float_t wids = fCurrentCoef[1];
3804 Float_t widfE = fCurrentCoefE;
3806 (* fDebugStreamer) << "FillFillPRF"<<
3807 "detector="<<detector<<
3809 "caligroup="<<caligroup<<
3821 //________________________________________________________________________________
3822 void AliTRDCalibraFit::FillFillLinearFitter()
3825 // DebugStream and fVectorFit
3828 // End of one detector
3834 for (Int_t k = 0; k < 2304; k++) {
3835 fCurrentCoefDetector[k] = 0.0;
3836 fCurrentCoefDetector2[k] = 0.0;
3840 if(fDebugLevel > 1){
3842 if ( !fDebugStreamer ) {
3844 TDirectory *backup = gDirectory;
3845 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3846 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3849 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3850 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3851 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3852 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3853 Float_t tiltangle = padplane->GetTiltingAngle();
3854 Int_t detector = fCountDet;
3855 Int_t stack = GetStack(fCountDet);
3856 Int_t layer = GetLayer(fCountDet);
3857 Float_t vf = fCurrentCoef[0];
3858 Float_t vs = fCurrentCoef[1];
3859 Float_t vfE = fCurrentCoefE;
3860 Float_t lorentzangler = fCurrentCoef2[0];
3861 Float_t elorentzangler = fCurrentCoefE2;
3862 Float_t lorentzangles = fCurrentCoef2[1];
3864 (* fDebugStreamer) << "FillFillLinearFitter"<<
3865 "detector="<<detector<<
3870 "tiltangle="<<tiltangle<<
3874 "lorentzangler="<<lorentzangler<<
3875 "Elorentzangler="<<elorentzangler<<
3876 "lorentzangles="<<lorentzangles<<
3882 //____________Calcul Coef Mean_________________________________________________
3884 //_____________________________________________________________________________
3885 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3888 // For the detector Dect calcul the mean time 0
3889 // for the calibration group idect from the choosen database
3892 fCurrentCoef2[1] = 0.0;
3893 if(fDebugLevel != 1){
3894 if(((fCalibraMode->GetNz(1) > 0) ||
3895 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3897 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3898 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3899 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3903 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3909 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3913 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3914 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3915 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3918 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3926 //_____________________________________________________________________________
3927 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3930 // For the detector Dect calcul the mean gain factor
3931 // for the calibration group idect from the choosen database
3934 fCurrentCoef[1] = 0.0;
3935 if(fDebugLevel != 1){
3936 if (((fCalibraMode->GetNz(0) > 0) ||
3937 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
3938 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3939 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3940 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3941 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3944 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3948 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3949 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3954 //_____________________________________________________________________________
3955 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3958 // For the detector Dect calcul the mean sigma of pad response
3959 // function for the calibration group idect from the choosen database
3962 fCurrentCoef[1] = 0.0;
3963 if(fDebugLevel != 1){
3964 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3965 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3966 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3969 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3973 //_____________________________________________________________________________
3974 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3977 // For the detector dect calcul the mean drift velocity for the
3978 // calibration group idect from the choosen database
3981 fCurrentCoef[1] = 0.0;
3982 if(fDebugLevel != 1){
3983 if (((fCalibraMode->GetNz(1) > 0) ||
3984 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
3986 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3987 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3988 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3992 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3997 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4002 //_____________________________________________________________________________
4003 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4006 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4009 fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
4010 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4014 //_____________________________________________________________________________
4015 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4018 // Default width of the PRF if there is no database as reference
4023 //case 0: return 0.515;
4024 //case 1: return 0.502;
4025 //case 2: return 0.491;
4026 //case 3: return 0.481;
4027 //case 4: return 0.471;
4028 //case 5: return 0.463;
4029 //default: return 0.0;
4032 case 0: return 0.538429;
4033 case 1: return 0.524302;
4034 case 2: return 0.511591;
4035 case 3: return 0.500140;
4036 case 4: return 0.489821;
4037 case 5: return 0.480524;
4038 default: return 0.0;
4041 //________________________________________________________________________________
4042 void AliTRDCalibraFit::SetCalROC(Int_t i)
4045 // Set the calib object for fCountDet
4048 Float_t value = 0.0;
4050 //Get the CalDet object
4052 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4054 AliInfo("Could not get calibDB");
4061 fCalROC->~AliTRDCalROC();
4062 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4063 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4067 fCalROC->~AliTRDCalROC();
4068 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4069 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4071 fCalROC2->~AliTRDCalROC();
4072 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4073 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4077 fCalROC->~AliTRDCalROC();
4078 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4079 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4088 if(fCalROC) delete fCalROC;
4089 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4090 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4091 fCalROC->SetValue(k,1.0);
4095 if(fCalROC) delete fCalROC;
4096 if(fCalROC2) delete fCalROC2;
4097 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4098 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4099 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4100 fCalROC->SetValue(k,1.0);
4101 fCalROC2->SetValue(k,0.0);
4105 if(fCalROC) delete fCalROC;
4106 value = GetPRFDefault(GetLayer(fCountDet));
4107 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4108 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4109 fCalROC->SetValue(k,value);
4117 //____________Fit Methods______________________________________________________
4119 //_____________________________________________________________________________
4120 void AliTRDCalibraFit::FitPente(TH1* projPH)
4123 // Slope methode for the drift velocity
4127 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4134 fCurrentCoefE = 0.0;
4135 fCurrentCoefE2 = 0.0;
4136 fCurrentCoef[0] = 0.0;
4137 fCurrentCoef2[0] = 0.0;
4138 TLine *line = new TLine();
4141 TAxis *xpph = projPH->GetXaxis();
4142 Int_t nbins = xpph->GetNbins();
4143 Double_t lowedge = xpph->GetBinLowEdge(1);
4144 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4145 Double_t widbins = (upedge - lowedge) / nbins;
4146 Double_t limit = upedge + 0.5 * widbins;
4149 // Beginning of the signal
4150 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4151 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4152 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4154 binmax = (Int_t) pentea->GetMaximumBin();
4157 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4159 if (binmax >= nbins) {
4162 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4164 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4165 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4166 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4167 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4168 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4170 fPhd[0] = -(l3P1am / (2 * l3P2am));
4173 if((l3P1am != 0.0) && (l3P2am != 0.0)){
4174 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4177 // Amplification region
4180 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4181 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4188 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4190 if (binmax >= nbins) {
4193 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4195 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4196 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4197 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4198 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4199 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4201 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4203 if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
4204 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4207 fCurrentCoefE2 = fCurrentCoefE;
4210 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4211 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4212 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4215 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4218 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4220 if (binmin >= nbins) {
4223 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4228 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4229 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4230 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4231 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4232 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4233 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4235 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4237 if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
4238 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4240 Float_t fPhdt0 = 0.0;
4241 Float_t t0Shift = 0.0;
4244 t0Shift = fT0Shift1;
4248 t0Shift = fT0Shift0;
4251 if ((fPhd[2] > fPhd[0]) &&
4252 (fPhd[2] > fPhd[1]) &&
4253 (fPhd[1] > fPhd[0]) &&
4255 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4256 fNumberFitSuccess++;
4258 if (fPhdt0 >= 0.0) {
4259 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4260 if (fCurrentCoef2[0] < -1.0) {
4261 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4265 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4270 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4271 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4274 if (fDebugLevel == 1) {
4275 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4278 line->SetLineColor(2);
4279 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4280 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4281 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4282 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4283 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4284 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4285 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4286 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4289 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4298 //_____________________________________________________________________________
4299 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4302 // Slope methode but with polynomes de Lagrange
4306 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4309 Double_t *x = new Double_t[5];
4310 Double_t *y = new Double_t[5];
4325 fCurrentCoefE = 0.0;
4326 fCurrentCoefE2 = 1.0;
4327 fCurrentCoef[0] = 0.0;
4328 fCurrentCoef2[0] = 0.0;
4329 TLine *line = new TLine();
4330 TF1 * polynome = 0x0;
4331 TF1 * polynomea = 0x0;
4332 TF1 * polynomeb = 0x0;
4336 TAxis *xpph = projPH->GetXaxis();
4337 Int_t nbins = xpph->GetNbins();
4338 Double_t lowedge = xpph->GetBinLowEdge(1);
4339 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4340 Double_t widbins = (upedge - lowedge) / nbins;
4341 Double_t limit = upedge + 0.5 * widbins;
4346 // Beginning of the signal
4347 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4348 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4349 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4352 binmax = (Int_t) pentea->GetMaximumBin();
4354 Double_t minnn = 0.0;
4355 Double_t maxxx = 0.0;
4357 Int_t kase = nbins-binmax;
4365 minnn = pentea->GetBinCenter(binmax-2);
4366 maxxx = pentea->GetBinCenter(binmax);
4367 x[0] = pentea->GetBinCenter(binmax-2);
4368 x[1] = pentea->GetBinCenter(binmax-1);
4369 x[2] = pentea->GetBinCenter(binmax);
4370 y[0] = pentea->GetBinContent(binmax-2);
4371 y[1] = pentea->GetBinContent(binmax-1);
4372 y[2] = pentea->GetBinContent(binmax);
4373 c = CalculPolynomeLagrange2(x,y);
4374 AliInfo("At the limit for beginning!");
4377 minnn = pentea->GetBinCenter(binmax-2);
4378 maxxx = pentea->GetBinCenter(binmax+1);
4379 x[0] = pentea->GetBinCenter(binmax-2);
4380 x[1] = pentea->GetBinCenter(binmax-1);
4381 x[2] = pentea->GetBinCenter(binmax);
4382 x[3] = pentea->GetBinCenter(binmax+1);
4383 y[0] = pentea->GetBinContent(binmax-2);
4384 y[1] = pentea->GetBinContent(binmax-1);
4385 y[2] = pentea->GetBinContent(binmax);
4386 y[3] = pentea->GetBinContent(binmax+1);
4387 c = CalculPolynomeLagrange3(x,y);
4395 minnn = pentea->GetBinCenter(binmax);
4396 maxxx = pentea->GetBinCenter(binmax+2);
4397 x[0] = pentea->GetBinCenter(binmax);
4398 x[1] = pentea->GetBinCenter(binmax+1);
4399 x[2] = pentea->GetBinCenter(binmax+2);
4400 y[0] = pentea->GetBinContent(binmax);
4401 y[1] = pentea->GetBinContent(binmax+1);
4402 y[2] = pentea->GetBinContent(binmax+2);
4403 c = CalculPolynomeLagrange2(x,y);
4406 minnn = pentea->GetBinCenter(binmax-1);
4407 maxxx = pentea->GetBinCenter(binmax+2);
4408 x[0] = pentea->GetBinCenter(binmax-1);
4409 x[1] = pentea->GetBinCenter(binmax);
4410 x[2] = pentea->GetBinCenter(binmax+1);
4411 x[3] = pentea->GetBinCenter(binmax+2);
4412 y[0] = pentea->GetBinContent(binmax-1);
4413 y[1] = pentea->GetBinContent(binmax);
4414 y[2] = pentea->GetBinContent(binmax+1);
4415 y[3] = pentea->GetBinContent(binmax+2);
4416 c = CalculPolynomeLagrange3(x,y);
4419 minnn = pentea->GetBinCenter(binmax-2);
4420 maxxx = pentea->GetBinCenter(binmax+2);
4421 x[0] = pentea->GetBinCenter(binmax-2);
4422 x[1] = pentea->GetBinCenter(binmax-1);
4423 x[2] = pentea->GetBinCenter(binmax);
4424 x[3] = pentea->GetBinCenter(binmax+1);
4425 x[4] = pentea->GetBinCenter(binmax+2);
4426 y[0] = pentea->GetBinContent(binmax-2);
4427 y[1] = pentea->GetBinContent(binmax-1);
4428 y[2] = pentea->GetBinContent(binmax);
4429 y[3] = pentea->GetBinContent(binmax+1);
4430 y[4] = pentea->GetBinContent(binmax+2);
4431 c = CalculPolynomeLagrange4(x,y);
4439 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4440 polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4442 Double_t step = (maxxx-minnn)/10000;
4444 Double_t maxvalue = 0.0;
4445 Double_t placemaximum = minnn;
4446 for(Int_t o = 0; o < 10000; o++){
4447 if(o == 0) maxvalue = polynomeb->Eval(l);
4448 if(maxvalue < (polynomeb->Eval(l))){
4449 maxvalue = polynomeb->Eval(l);
4454 fPhd[0] = placemaximum;
4457 // Amplification region
4460 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4461 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4467 Double_t minn = 0.0;
4468 Double_t maxx = 0.0;
4480 Int_t kase1 = nbins - binmax;
4482 //Determination of minn and maxx
4483 //case binmax = nbins
4488 minn = projPH->GetBinCenter(binmax-2);
4489 maxx = projPH->GetBinCenter(binmax);
4490 x[0] = projPH->GetBinCenter(binmax-2);
4491 x[1] = projPH->GetBinCenter(binmax-1);
4492 x[2] = projPH->GetBinCenter(binmax);
4493 y[0] = projPH->GetBinContent(binmax-2);
4494 y[1] = projPH->GetBinContent(binmax-1);
4495 y[2] = projPH->GetBinContent(binmax);
4496 c = CalculPolynomeLagrange2(x,y);
4497 //AliInfo("At the limit for the drift!");
4500 minn = projPH->GetBinCenter(binmax-2);
4501 maxx = projPH->GetBinCenter(binmax+1);
4502 x[0] = projPH->GetBinCenter(binmax-2);
4503 x[1] = projPH->GetBinCenter(binmax-1);
4504 x[2] = projPH->GetBinCenter(binmax);
4505 x[3] = projPH->GetBinCenter(binmax+1);
4506 y[0] = projPH->GetBinContent(binmax-2);
4507 y[1] = projPH->GetBinContent(binmax-1);
4508 y[2] = projPH->GetBinContent(binmax);
4509 y[3] = projPH->GetBinContent(binmax+1);
4510 c = CalculPolynomeLagrange3(x,y);
4519 minn = projPH->GetBinCenter(binmax);
4520 maxx = projPH->GetBinCenter(binmax+2);
4521 x[0] = projPH->GetBinCenter(binmax);
4522 x[1] = projPH->GetBinCenter(binmax+1);
4523 x[2] = projPH->GetBinCenter(binmax+2);
4524 y[0] = projPH->GetBinContent(binmax);
4525 y[1] = projPH->GetBinContent(binmax+1);
4526 y[2] = projPH->GetBinContent(binmax+2);
4527 c = CalculPolynomeLagrange2(x,y);
4530 minn = projPH->GetBinCenter(binmax-1);
4531 maxx = projPH->GetBinCenter(binmax+2);
4532 x[0] = projPH->GetBinCenter(binmax-1);
4533 x[1] = projPH->GetBinCenter(binmax);
4534 x[2] = projPH->GetBinCenter(binmax+1);
4535 x[3] = projPH->GetBinCenter(binmax+2);
4536 y[0] = projPH->GetBinContent(binmax-1);
4537 y[1] = projPH->GetBinContent(binmax);
4538 y[2] = projPH->GetBinContent(binmax+1);
4539 y[3] = projPH->GetBinContent(binmax+2);
4540 c = CalculPolynomeLagrange3(x,y);
4543 minn = projPH->GetBinCenter(binmax-2);
4544 maxx = projPH->GetBinCenter(binmax+2);
4545 x[0] = projPH->GetBinCenter(binmax-2);
4546 x[1] = projPH->GetBinCenter(binmax-1);
4547 x[2] = projPH->GetBinCenter(binmax);
4548 x[3] = projPH->GetBinCenter(binmax+1);
4549 x[4] = projPH->GetBinCenter(binmax+2);
4550 y[0] = projPH->GetBinContent(binmax-2);
4551 y[1] = projPH->GetBinContent(binmax-1);
4552 y[2] = projPH->GetBinContent(binmax);
4553 y[3] = projPH->GetBinContent(binmax+1);
4554 y[4] = projPH->GetBinContent(binmax+2);
4555 c = CalculPolynomeLagrange4(x,y);
4562 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4563 polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4565 Double_t step = (maxx-minn)/1000;
4567 Double_t maxvalue = 0.0;
4568 Double_t placemaximum = minn;
4569 for(Int_t o = 0; o < 1000; o++){
4570 if(o == 0) maxvalue = polynomea->Eval(l);
4571 if(maxvalue < (polynomea->Eval(l))){
4572 maxvalue = polynomea->Eval(l);
4577 fPhd[1] = placemaximum;
4581 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4582 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4583 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4586 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4592 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4596 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4597 AliInfo("Too many fluctuations at the end!");
4600 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4601 AliInfo("Too many fluctuations at the end!");
4604 if(pente->GetBinContent(binmin+1)==0){
4605 AliInfo("No entries for the next bin!");
4606 pente->SetBinContent(binmin,0);
4607 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4623 Bool_t case1 = kFALSE;
4624 Bool_t case2 = kFALSE;
4625 Bool_t case4 = kFALSE;
4627 //Determination of min and max
4628 //case binmin <= nbins-3
4630 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4631 min = pente->GetBinCenter(binmin-2);
4632 max = pente->GetBinCenter(binmin+2);
4633 x[0] = pente->GetBinCenter(binmin-2);
4634 x[1] = pente->GetBinCenter(binmin-1);
4635 x[2] = pente->GetBinCenter(binmin);
4636 x[3] = pente->GetBinCenter(binmin+1);
4637 x[4] = pente->GetBinCenter(binmin+2);
4638 y[0] = pente->GetBinContent(binmin-2);
4639 y[1] = pente->GetBinContent(binmin-1);
4640 y[2] = pente->GetBinContent(binmin);
4641 y[3] = pente->GetBinContent(binmin+1);
4642 y[4] = pente->GetBinContent(binmin+2);
4643 //Calcul the polynome de Lagrange
4644 c = CalculPolynomeLagrange4(x,y);
4646 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4647 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4648 //AliInfo("polynome 4 false 1");
4651 if(((binmin+3) <= (nbins-1)) &&
4652 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4653 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4654 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4655 AliInfo("polynome 4 false 2");
4659 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4660 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4661 //AliInfo("polynome 4 case 1");
4664 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4665 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4666 //AliInfo("polynome 4 case 4");
4671 //case binmin = nbins-2
4673 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4675 min = pente->GetBinCenter(binmin-2);
4676 max = pente->GetBinCenter(binmin+1);
4677 x[0] = pente->GetBinCenter(binmin-2);
4678 x[1] = pente->GetBinCenter(binmin-1);
4679 x[2] = pente->GetBinCenter(binmin);
4680 x[3] = pente->GetBinCenter(binmin+1);
4681 y[0] = pente->GetBinContent(binmin-2);
4682 y[1] = pente->GetBinContent(binmin-1);
4683 y[2] = pente->GetBinContent(binmin);
4684 y[3] = pente->GetBinContent(binmin+1);
4685 //Calcul the polynome de Lagrange
4686 c = CalculPolynomeLagrange3(x,y);
4687 //richtung +: nothing
4689 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4690 //AliInfo("polynome 3- case 2");
4695 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4697 min = pente->GetBinCenter(binmin-1);
4698 max = pente->GetBinCenter(binmin+2);
4699 x[0] = pente->GetBinCenter(binmin-1);
4700 x[1] = pente->GetBinCenter(binmin);
4701 x[2] = pente->GetBinCenter(binmin+1);
4702 x[3] = pente->GetBinCenter(binmin+2);
4703 y[0] = pente->GetBinContent(binmin-1);
4704 y[1] = pente->GetBinContent(binmin);
4705 y[2] = pente->GetBinContent(binmin+1);
4706 y[3] = pente->GetBinContent(binmin+2);
4707 //Calcul the polynome de Lagrange
4708 c = CalculPolynomeLagrange3(x,y);
4710 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4711 //AliInfo("polynome 3+ case 2");
4716 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4717 min = pente->GetBinCenter(binmin);
4718 max = pente->GetBinCenter(binmin+2);
4719 x[0] = pente->GetBinCenter(binmin);
4720 x[1] = pente->GetBinCenter(binmin+1);
4721 x[2] = pente->GetBinCenter(binmin+2);
4722 y[0] = pente->GetBinContent(binmin);
4723 y[1] = pente->GetBinContent(binmin+1);
4724 y[2] = pente->GetBinContent(binmin+2);
4725 //Calcul the polynome de Lagrange
4726 c = CalculPolynomeLagrange2(x,y);
4728 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4729 //AliInfo("polynome 2+ false");
4734 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4736 min = pente->GetBinCenter(binmin-1);
4737 max = pente->GetBinCenter(binmin+1);
4738 x[0] = pente->GetBinCenter(binmin-1);
4739 x[1] = pente->GetBinCenter(binmin);
4740 x[2] = pente->GetBinCenter(binmin+1);
4741 y[0] = pente->GetBinContent(binmin-1);
4742 y[1] = pente->GetBinContent(binmin);
4743 y[2] = pente->GetBinContent(binmin+1);
4744 //Calcul the polynome de Lagrange
4745 c = CalculPolynomeLagrange2(x,y);
4746 //richtung +: nothing
4747 //richtung -: nothing
4749 //case binmin = nbins-1
4751 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4752 min = pente->GetBinCenter(binmin-2);
4753 max = pente->GetBinCenter(binmin);
4754 x[0] = pente->GetBinCenter(binmin-2);
4755 x[1] = pente->GetBinCenter(binmin-1);
4756 x[2] = pente->GetBinCenter(binmin);
4757 y[0] = pente->GetBinContent(binmin-2);
4758 y[1] = pente->GetBinContent(binmin-1);
4759 y[2] = pente->GetBinContent(binmin);
4760 //Calcul the polynome de Lagrange
4761 c = CalculPolynomeLagrange2(x,y);
4762 //AliInfo("At the limit for the drift!");
4763 //fluctuation too big!
4764 //richtung +: nothing
4766 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4767 //AliInfo("polynome 2- false ");
4771 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4773 AliInfo("At the limit for the drift and not usable!");
4777 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4779 AliInfo("For the drift...problem!");
4781 //pass but should not happen
4782 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4784 AliInfo("For the drift...problem!");
4788 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4789 polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4790 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4791 Double_t step = (max-min)/1000;
4793 Double_t minvalue = 0.0;
4794 Double_t placeminimum = min;
4795 for(Int_t o = 0; o < 1000; o++){
4796 if(o == 0) minvalue = polynome->Eval(l);
4797 if(minvalue > (polynome->Eval(l))){
4798 minvalue = polynome->Eval(l);
4803 fPhd[2] = placeminimum;
4805 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4806 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4807 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4809 Float_t fPhdt0 = 0.0;
4810 Float_t t0Shift = 0.0;
4813 t0Shift = fT0Shift1;
4817 t0Shift = fT0Shift0;
4820 if ((fPhd[2] > fPhd[0]) &&
4821 (fPhd[2] > fPhd[1]) &&
4822 (fPhd[1] > fPhd[0]) &&
4824 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4825 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4826 else fNumberFitSuccess++;
4827 if (fPhdt0 >= 0.0) {
4828 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4829 if (fCurrentCoef2[0] < -1.0) {
4830 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4834 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4838 //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4839 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4841 if((fPhd[1] > fPhd[0]) &&
4843 if (fPhdt0 >= 0.0) {
4844 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4845 if (fCurrentCoef2[0] < -1.0) {
4846 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4850 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4854 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4855 //printf("Fit failed!\n");
4859 if (fDebugLevel == 1) {
4860 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4863 line->SetLineColor(2);
4864 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4865 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4866 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4867 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4868 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4869 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4870 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4871 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4874 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4879 if(pentea) delete pentea;
4880 if(pente) delete pente;
4881 if(polynome) delete polynome;
4882 if(polynomea) delete polynomea;
4883 if(polynomeb) delete polynomeb;
4887 if(line) delete line;
4892 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4893 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4895 projPH->SetDirectory(0);
4899 //_____________________________________________________________________________
4900 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4903 // Fit methode for the drift velocity
4907 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4910 TAxis *xpph = projPH->GetXaxis();
4911 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4913 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4914 fPH->SetParameter(0,0.469); // Scaling
4915 fPH->SetParameter(1,0.18); // Start
4916 fPH->SetParameter(2,0.0857325); // AR
4917 fPH->SetParameter(3,1.89); // DR
4918 fPH->SetParameter(4,0.08); // QA/QD
4919 fPH->SetParameter(5,0.0); // Baseline
4921 TLine *line = new TLine();
4923 fCurrentCoef[0] = 0.0;
4924 fCurrentCoef2[0] = 0.0;
4925 fCurrentCoefE = 0.0;
4926 fCurrentCoefE2 = 0.0;
4928 if (idect%fFitPHPeriode == 0) {
4930 AliInfo(Form("The detector %d will be fitted",idect));
4931 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4932 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
4933 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
4934 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
4935 fPH->SetParameter(4,0.225); // QA/QD
4936 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4938 if (fDebugLevel != 1) {
4939 projPH->Fit(fPH,"0M","",0.0,upedge);
4942 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4944 projPH->Fit(fPH,"M+","",0.0,upedge);
4946 line->SetLineColor(4);
4947 line->DrawLine(fPH->GetParameter(1)
4949 ,fPH->GetParameter(1)
4950 ,projPH->GetMaximum());
4951 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4953 ,fPH->GetParameter(1)+fPH->GetParameter(2)
4954 ,projPH->GetMaximum());
4955 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4957 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4958 ,projPH->GetMaximum());
4961 if (fPH->GetParameter(3) != 0) {
4962 fNumberFitSuccess++;
4963 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
4964 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4965 fCurrentCoef2[0] = fPH->GetParameter(1);
4966 fCurrentCoefE2 = fPH->GetParError(1);
4969 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4970 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4976 // Put the default value
4977 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4978 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4981 if (fDebugLevel != 1) {
4986 //_____________________________________________________________________________
4987 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
4990 // Fit methode for the sigma of the pad response function
4995 fCurrentCoef[0] = 0.0;
4996 fCurrentCoefE = 0.0;
4998 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5001 fCurrentCoef[0] = -fCurrentCoef[1];
5005 fNumberFitSuccess++;
5006 fCurrentCoef[0] = param[2];
5007 fCurrentCoefE = ret;
5011 //_____________________________________________________________________________
5012 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)
5015 // Fit methode for the sigma of the pad response function
5018 //We should have at least 3 points
5019 if(nBins <=3) return -4.0;
5021 TLinearFitter fitter(3,"pol2");
5022 fitter.StoreData(kFALSE);
5023 fitter.ClearPoints();
5025 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5026 Float_t entries = 0;
5027 Int_t nbbinwithentries = 0;
5028 for (Int_t i=0; i<nBins; i++){
5030 if(arraye[i] > 15) nbbinwithentries++;
5031 //printf("entries for i %d: %f\n",i,arraye[i]);
5033 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5034 //printf("entries %f\n",entries);
5035 //printf("nbbinwithentries %d\n",nbbinwithentries);
5038 Float_t errorm = 0.0;
5039 Float_t errorn = 0.0;
5040 Float_t error = 0.0;
5043 for (Int_t ibin=0;ibin<nBins; ibin++){
5044 Float_t entriesI = arraye[ibin];
5045 Float_t valueI = arraym[ibin];
5046 Double_t xcenter = 0.0;
5048 if ((entriesI>15) && (valueI>0.0)){
5049 xcenter = xMin+(ibin+0.5)*binWidth;
5054 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5055 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5056 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5059 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5060 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5061 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5063 if(error == 0.0) continue;
5064 val = TMath::Log(Float_t(valueI));
5065 fitter.AddPoint(&xcenter,val,error);
5069 if(fDebugLevel > 1){
5071 if ( !fDebugStreamer ) {
5073 TDirectory *backup = gDirectory;
5074 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5075 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5078 Int_t detector = fCountDet;
5079 Int_t layer = GetLayer(fCountDet);
5082 (* fDebugStreamer) << "FitGausMIFill"<<
5083 "detector="<<detector<<
5087 "entriesI="<<entriesI<<
5090 "xcenter="<<xcenter<<
5100 if(npoints <=3) return -4.0;
5105 fitter.GetParameters(par);
5106 chi2 = fitter.GetChisquare()/Float_t(npoints);
5109 if (!param) param = new TVectorD(3);
5110 if(par[2] == 0.0) return -4.0;
5111 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5112 Double_t deltax = (fitter.GetParError(2))/x;
5113 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5116 (*param)[1] = par[1]/(-2.*par[2]);
5117 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5118 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5119 if ( lnparam0>307 ) return -4;
5120 (*param)[0] = TMath::Exp(lnparam0);
5122 if(fDebugLevel > 1){
5124 if ( !fDebugStreamer ) {
5126 TDirectory *backup = gDirectory;
5127 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5128 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5131 Int_t detector = fCountDet;
5132 Int_t layer = GetLayer(fCountDet);
5135 (* fDebugStreamer) << "FitGausMIFit"<<
5136 "detector="<<detector<<
5139 "errorsigma="<<chi2<<
5140 "mean="<<(*param)[1]<<
5141 "sigma="<<(*param)[2]<<
5142 "constant="<<(*param)[0]<<
5147 if((chi2/(*param)[2]) > 0.1){
5149 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5154 if(fDebugLevel == 1){
5155 TString name("PRF");
5156 name += (Int_t)xMin;
5157 name += (Int_t)xMax;
5158 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5161 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5162 for(Int_t k = 0; k < nBins; k++){
5163 histo->SetBinContent(k+1,arraym[k]);
5164 histo->SetBinError(k+1,arrayme[k]);
5167 name += "functionf";
5168 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5169 f1->SetParameter(0, (*param)[0]);
5170 f1->SetParameter(1, (*param)[1]);
5171 f1->SetParameter(2, (*param)[2]);
5179 //_____________________________________________________________________________
5180 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5183 // Fit methode for the sigma of the pad response function
5186 fCurrentCoef[0] = 0.0;
5187 fCurrentCoefE = 0.0;
5189 if (fDebugLevel != 1) {
5190 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5193 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5195 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5198 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5199 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5200 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5202 fNumberFitSuccess++;
5205 //_____________________________________________________________________________
5206 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5209 // Fit methode for the sigma of the pad response function
5211 fCurrentCoef[0] = 0.0;
5212 fCurrentCoefE = 0.0;
5213 if (fDebugLevel == 1) {
5214 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5218 fCurrentCoef[0] = projPRF->GetRMS();
5219 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5221 fNumberFitSuccess++;
5224 //_____________________________________________________________________________
5225 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5228 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5231 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5234 Int_t nbins = (Int_t)(nybins/(2*nbg));
5235 Float_t lowedge = -3.0*nbg;
5236 Float_t upedge = lowedge + 3.0;
5239 Double_t xvalues = -0.2*nbg+0.1;
5241 Int_t total = 2*nbg;
5244 for(Int_t k = 0; k < total; k++){
5245 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5247 y = fCurrentCoef[0]*fCurrentCoef[0];
5248 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5251 if(fDebugLevel > 1){
5253 if ( !fDebugStreamer ) {
5255 TDirectory *backup = gDirectory;
5256 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5257 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5260 Int_t detector = fCountDet;
5261 Int_t layer = GetLayer(fCountDet);
5262 Int_t nbtotal = total;
5264 Float_t low = lowedge;
5265 Float_t up = upedge;
5266 Float_t tnp = xvalues;
5267 Float_t wid = fCurrentCoef[0];
5268 Float_t widfE = fCurrentCoefE;
5270 (* fDebugStreamer) << "FitTnpRange0"<<
5271 "detector="<<detector<<
5273 "nbtotal="<<nbtotal<<
5291 fCurrentCoefE = 0.0;
5292 fCurrentCoef[0] = 0.0;
5294 //printf("npoints\n",npoints);
5297 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5302 linearfitter.Eval();
5303 linearfitter.GetParameters(pars0);
5304 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5305 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5306 Double_t min0 = 0.0;
5307 Double_t ermin0 = 0.0;
5308 //Double_t prfe0 = 0.0;
5309 Double_t prf0 = 0.0;
5310 if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
5311 min0 = -pars0[1]/(2*pars0[2]);
5312 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5313 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5316 prfe0 = linearfitter->GetParError(0)*pointError0
5317 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5318 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5319 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5320 fCurrentCoefE = (Float_t) prfe0;
5322 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5325 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5329 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5332 if(fDebugLevel > 1){
5334 if ( !fDebugStreamer ) {
5336 TDirectory *backup = gDirectory;
5337 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5338 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5341 Int_t detector = fCountDet;
5342 Int_t layer = GetLayer(fCountDet);
5343 Int_t nbtotal = total;
5344 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5345 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5347 (* fDebugStreamer) << "FitTnpRange1"<<
5348 "detector="<<detector<<
5350 "nbtotal="<<nbtotal<<
5354 "npoints="<<npoints<<
5357 "sigmaprf="<<fCurrentCoef[0]<<
5358 "sigprf="<<fCurrentCoef[1]<<
5365 //_____________________________________________________________________________
5366 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5369 // Only mean methode for the gain factor
5372 fCurrentCoef[0] = mean;
5373 fCurrentCoefE = 0.0;
5374 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5375 if (fDebugLevel == 1) {
5376 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5380 CalculChargeCoefMean(kTRUE);
5381 fNumberFitSuccess++;
5383 //_____________________________________________________________________________
5384 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5387 // mean w methode for the gain factor
5391 Int_t nybins = projch->GetNbinsX();
5393 //The weight function
5394 Double_t a = 0.00228515;
5395 Double_t b = -0.00231487;
5396 Double_t c = 0.00044298;
5397 Double_t d = -0.00379239;
5398 Double_t e = 0.00338349;
5408 //A arbitrary error for the moment
5409 fCurrentCoefE = 0.0;
5410 fCurrentCoef[0] = 0.0;
5413 Double_t sumw = 0.0;
5415 Float_t sumAll = (Float_t) nentries;
5416 Int_t sumCurrent = 0;
5417 for(Int_t k = 0; k <nybins; k++){
5418 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5419 if (fraction>0.95) break;
5420 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5421 e*fraction*fraction*fraction*fraction;
5422 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5423 sum += weight*projch->GetBinContent(k+1);
5424 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5425 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5427 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5429 if (fDebugLevel == 1) {
5430 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5434 fNumberFitSuccess++;
5435 CalculChargeCoefMean(kTRUE);
5437 //_____________________________________________________________________________
5438 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5441 // mean w methode for the gain factor
5445 Int_t nybins = projch->GetNbinsX();
5447 //The weight function
5448 Double_t a = 0.00228515;
5449 Double_t b = -0.00231487;
5450 Double_t c = 0.00044298;
5451 Double_t d = -0.00379239;
5452 Double_t e = 0.00338349;
5462 //A arbitrary error for the moment
5463 fCurrentCoefE = 0.0;
5464 fCurrentCoef[0] = 0.0;
5467 Double_t sumw = 0.0;
5469 Int_t sumCurrent = 0;
5470 for(Int_t k = 0; k <nybins; k++){
5471 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5472 if (fraction>0.95) break;
5473 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5474 e*fraction*fraction*fraction*fraction;
5475 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5476 sum += weight*projch->GetBinContent(k+1);
5477 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5478 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5480 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5482 if (fDebugLevel == 1) {
5483 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5487 fNumberFitSuccess++;
5489 //_____________________________________________________________________________
5490 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5493 // Fit methode for the gain factor
5496 fCurrentCoef[0] = 0.0;
5497 fCurrentCoefE = 0.0;
5498 Double_t chisqrl = 0.0;
5499 Double_t chisqrg = 0.0;
5500 Double_t chisqr = 0.0;
5501 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5503 projch->Fit("landau","0",""
5504 ,(Double_t) mean/fBeginFitCharge
5505 ,projch->GetBinCenter(projch->GetNbinsX()));
5506 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
5507 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
5508 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
5509 chisqrl = projch->GetFunction("landau")->GetChisquare();
5511 projch->Fit("gaus","0",""
5512 ,(Double_t) mean/fBeginFitCharge
5513 ,projch->GetBinCenter(projch->GetNbinsX()));
5514 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
5515 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
5516 chisqrg = projch->GetFunction("gaus")->GetChisquare();
5518 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5519 if (fDebugLevel != 1) {
5520 projch->Fit("fLandauGaus","0",""
5521 ,(Double_t) mean/fBeginFitCharge
5522 ,projch->GetBinCenter(projch->GetNbinsX()));
5523 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5526 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5528 projch->Fit("fLandauGaus","+",""
5529 ,(Double_t) mean/fBeginFitCharge
5530 ,projch->GetBinCenter(projch->GetNbinsX()));
5531 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5533 fLandauGaus->Draw("same");
5536 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5537 //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5538 fNumberFitSuccess++;
5539 CalculChargeCoefMean(kTRUE);
5540 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
5541 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
5544 CalculChargeCoefMean(kFALSE);
5545 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5548 if (fDebugLevel != 1) {
5553 //_____________________________________________________________________________
5554 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5557 // Fit methode for the gain factor more time consuming
5561 //Some parameters to initialise
5562 Double_t widthLandau, widthGaus, mPV, integral;
5563 Double_t chisquarel = 0.0;
5564 Double_t chisquareg = 0.0;
5565 projch->Fit("landau","0M+",""
5567 ,projch->GetBinCenter(projch->GetNbinsX()));
5568 widthLandau = projch->GetFunction("landau")->GetParameter(2);
5569 chisquarel = projch->GetFunction("landau")->GetChisquare();
5570 projch->Fit("gaus","0M+",""
5572 ,projch->GetBinCenter(projch->GetNbinsX()));
5573 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
5574 chisquareg = projch->GetFunction("gaus")->GetChisquare();
5576 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5577 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5579 // Setting fit range and start values
5581 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5582 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5583 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
5584 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5585 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5586 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
5587 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
5590 fCurrentCoef[0] = 0.0;
5591 fCurrentCoefE = 0.0;
5595 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5600 Double_t projchPeak;
5601 Double_t projchFWHM;
5602 LanGauPro(fp,projchPeak,projchFWHM);
5604 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5605 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5606 fNumberFitSuccess++;
5607 CalculChargeCoefMean(kTRUE);
5608 fCurrentCoef[0] = fp[1];
5609 fCurrentCoefE = fpe[1];
5610 //chargeCoefE2 = chisqr;
5613 CalculChargeCoefMean(kFALSE);
5614 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5616 if (fDebugLevel == 1) {
5617 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5618 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5621 fitsnr->Draw("same");
5627 //_____________________________________________________________________________
5628 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
5631 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5633 Double_t *c = new Double_t[5];
5634 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5635 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5636 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5641 c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5642 c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5649 //_____________________________________________________________________________
5650 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
5653 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5655 Double_t *c = new Double_t[5];
5656 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5657 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5658 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5659 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5663 c[2] = -(x0*(x[1]+x[2]+x[3])
5664 +x1*(x[0]+x[2]+x[3])
5665 +x2*(x[0]+x[1]+x[3])
5666 +x3*(x[0]+x[1]+x[2]));
5667 c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5668 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5669 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5670 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5672 c[0] = -(x0*x[1]*x[2]*x[3]
5675 +x3*x[0]*x[1]*x[2]);
5683 //_____________________________________________________________________________
5684 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
5687 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5689 Double_t *c = new Double_t[5];
5690 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5691 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5692 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5693 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5694 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5697 c[4] = x0+x1+x2+x3+x4;
5698 c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
5699 +x1*(x[0]+x[2]+x[3]+x[4])
5700 +x2*(x[0]+x[1]+x[3]+x[4])
5701 +x3*(x[0]+x[1]+x[2]+x[4])
5702 +x4*(x[0]+x[1]+x[2]+x[3]));
5703 c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5704 +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])
5705 +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])
5706 +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])
5707 +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]));
5709 c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5710 +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])
5711 +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])
5712 +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])
5713 +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]));
5715 c[0] = (x0*x[1]*x[2]*x[3]*x[4]
5716 +x1*x[0]*x[2]*x[3]*x[4]
5717 +x2*x[0]*x[1]*x[3]*x[4]
5718 +x3*x[0]*x[1]*x[2]*x[4]
5719 +x4*x[0]*x[1]*x[2]*x[3]);
5725 //_____________________________________________________________________________
5726 void AliTRDCalibraFit::NormierungCharge()
5729 // Normalisation of the gain factor resulting for the fits
5732 // Calcul of the mean of choosen method by fFitChargeNDB
5734 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5735 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5737 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5738 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5739 //printf("detector %d coef[0] %f\n",detector,coef[0]);
5740 if (GetStack(detector) == 2) {
5743 if (GetStack(detector) != 2) {
5746 for (Int_t j = 0; j < total; j++) {
5754 fScaleFitFactor = fScaleFitFactor / sum;
5757 fScaleFitFactor = 1.0;
5760 //methode de boeuf mais bon...
5761 Double_t scalefactor = fScaleFitFactor;
5763 if(fDebugLevel > 1){
5765 if ( !fDebugStreamer ) {
5767 TDirectory *backup = gDirectory;
5768 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5769 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5771 (* fDebugStreamer) << "NormierungCharge"<<
5772 "scalefactor="<<scalefactor<<
5776 //_____________________________________________________________________________
5777 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5780 // Rebin of the 1D histo for the gain calibration if needed.
5781 // you have to choose fRebin, divider of fNumberBinCharge
5784 TAxis *xhist = hist->GetXaxis();
5785 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5786 ,xhist->GetBinLowEdge(1)
5787 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5789 AliInfo(Form("fRebin: %d",fRebin));
5791 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5793 for (Int_t ji = i; ji < i+fRebin; ji++) {
5794 sum += hist->GetBinContent(ji);
5797 rehist->SetBinContent(k,sum);
5805 //_____________________________________________________________________________
5806 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5809 // Rebin of the 1D histo for the gain calibration if needed
5810 // you have to choose fRebin divider of fNumberBinCharge
5813 TAxis *xhist = hist->GetXaxis();
5814 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5815 ,xhist->GetBinLowEdge(1)
5816 ,xhist->GetBinUpEdge(xhist->GetNbins()));
5818 AliInfo(Form("fRebin: %d",fRebin));
5820 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5822 for (Int_t ji = i; ji < i+fRebin; ji++) {
5823 sum += hist->GetBinContent(ji);
5826 rehist->SetBinContent(k,sum);
5834 //____________Some basic geometry function_____________________________________
5837 //_____________________________________________________________________________
5838 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5841 // Reconstruct the plane number from the detector number
5844 return ((Int_t) (d % 6));
5848 //_____________________________________________________________________________
5849 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5852 // Reconstruct the stack number from the detector number
5854 const Int_t kNlayer = 6;
5856 return ((Int_t) (d % 30) / kNlayer);
5860 //_____________________________________________________________________________
5861 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5864 // Reconstruct the sector number from the detector number
5868 return ((Int_t) (d / fg));
5873 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5875 //_______________________________________________________________________________
5876 void AliTRDCalibraFit::ResetVectorFit()
5879 // Reset the VectorFits
5882 fVectorFit.SetOwner();
5884 fVectorFit2.SetOwner();
5885 fVectorFit2.Clear();
5889 //____________Private Functions________________________________________________
5892 //_____________________________________________________________________________
5893 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5896 // Function for the fit
5899 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5901 //PARAMETERS FOR FIT PH
5903 //fAsymmGauss->SetParameter(0,0.113755);
5904 //fAsymmGauss->SetParameter(1,0.350706);
5905 //fAsymmGauss->SetParameter(2,0.0604244);
5906 //fAsymmGauss->SetParameter(3,7.65596);
5907 //fAsymmGauss->SetParameter(4,1.00124);
5908 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
5916 Double_t dx = 0.005;
5917 Double_t xs = par[1];
5919 Double_t paras[2] = { 0.0, 0.0 };
5922 if ((xs >= par[1]) &&
5923 (xs < (par[1]+par[2]))) {
5924 //fAsymmGauss->SetParameter(0,par[0]);
5925 //fAsymmGauss->SetParameter(1,xs);
5926 //ss += fAsymmGauss->Eval(xx);
5929 ss += AsymmGauss(&xx,paras);
5931 if ((xs >= (par[1]+par[2])) &&
5932 (xs < (par[1]+par[2]+par[3]))) {
5933 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5934 //fAsymmGauss->SetParameter(1,xs);
5935 //ss += fAsymmGauss->Eval(xx);
5936 paras[0] = par[0]*par[4];
5938 ss += AsymmGauss(&xx,paras);
5947 //_____________________________________________________________________________
5948 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5951 // Function for the fit
5954 //par[0] = normalization
5962 Double_t par1save = par[1];
5963 //Double_t par2save = par[2];
5964 Double_t par2save = 0.0604244;
5965 //Double_t par3save = par[3];
5966 Double_t par3save = 7.65596;
5967 //Double_t par5save = par[5];
5968 Double_t par5save = 0.870597;
5969 Double_t dx = x[0] - par1save;
5971 Double_t sigma2 = par2save*par2save;
5972 Double_t sqrt2 = TMath::Sqrt(2.0);
5973 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5974 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5975 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5976 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5978 //return par[0]*(exp1+par[4]*exp2);
5979 return par[0] * (exp1 + 1.00124 * exp2);
5983 //_____________________________________________________________________________
5984 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5987 // Sum Landau + Gaus with identical mean
5990 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5991 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5992 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5993 Double_t val = valLandau + valGaus;
5999 //_____________________________________________________________________________
6000 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6003 // Function for the fit
6006 // par[0]=Width (scale) parameter of Landau density
6007 // par[1]=Most Probable (MP, location) parameter of Landau density
6008 // par[2]=Total area (integral -inf to inf, normalization constant)
6009 // par[3]=Width (sigma) of convoluted Gaussian function
6011 // In the Landau distribution (represented by the CERNLIB approximation),
6012 // the maximum is located at x=-0.22278298 with the location parameter=0.
6013 // This shift is corrected within this function, so that the actual
6014 // maximum is identical to the MP parameter.
6017 // Numeric constants
6018 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6019 Double_t mpshift = -0.22278298; // Landau maximum location
6021 // Control constants
6022 Double_t np = 100.0; // Number of convolution steps
6023 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6035 // MP shift correction
6036 mpc = par[1] - mpshift * par[0];
6038 // Range of convolution integral
6039 xlow = x[0] - sc * par[3];
6040 xupp = x[0] + sc * par[3];
6042 step = (xupp - xlow) / np;
6044 // Convolution integral of Landau and Gaussian by sum
6045 for (i = 1.0; i <= np/2; i++) {
6047 xx = xlow + (i-.5) * step;
6048 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6049 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6051 xx = xupp - (i-.5) * step;
6052 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6053 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6057 return (par[2] * step * sum * invsq2pi / par[3]);
6060 //_____________________________________________________________________________
6061 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6062 , const Double_t *parlimitslo, const Double_t *parlimitshi
6063 , Double_t *fitparams, Double_t *fiterrors
6064 , Double_t *chiSqr, Int_t *ndf) const
6067 // Function for the fit
6071 Char_t funname[100];
6073 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6078 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6079 ffit->SetParameters(startvalues);
6080 ffit->SetParNames("Width","MP","Area","GSigma");
6082 for (i = 0; i < 4; i++) {
6083 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6086 his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
6088 ffit->GetParameters(fitparams); // Obtain fit parameters
6089 for (i = 0; i < 4; i++) {
6090 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6092 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6093 ndf[0] = ffit->GetNDF(); // Obtain ndf
6095 return (ffit); // Return fit function
6099 //_____________________________________________________________________________
6100 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6103 // Function for the fit
6116 Int_t maxcalls = 10000;
6118 // Search for maximum
6119 p = params[1] - 0.1 * params[0];
6120 step = 0.05 * params[0];
6124 while ((l != lold) && (i < maxcalls)) {
6128 l = LanGauFun(&x,params);
6130 step = -step / 10.0;
6135 if (i == maxcalls) {
6141 // Search for right x location of fy
6142 p = maxx + params[0];
6148 while ( (l != lold) && (i < maxcalls) ) {
6153 l = TMath::Abs(LanGauFun(&x,params) - fy);
6167 // Search for left x location of fy
6169 p = maxx - 0.5 * params[0];
6175 while ((l != lold) && (i < maxcalls)) {
6179 l = TMath::Abs(LanGauFun(&x,params) - fy);
6181 step = -step / 10.0;
6186 if (i == maxcalls) {
6195 //_____________________________________________________________________________
6196 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6199 // Gaus with identical mean
6202 Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];