1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
25 // if they have enough entries, otherwise the (default) value of the choosen database
26 // will be put. For the relative gain calibration the resulted factors will be globally
27 // normalized to the gain factors of the choosen database. It unables to precise
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
32 // _fDebug = 2: a comparaison of the coefficients found and the default values
33 // in the choosen database.
34 // fCoef , histogram of the coefs as function of the calibration group number
35 // fDelta , histogram of the relative difference of the coef with the default
36 // value in the database as function of the calibration group number
37 // fError , dirstribution of this relative difference
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 // pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 // also the comparaison histograms of the 1 for this detector
45 // R. Bailhache (R.Bailhache@gsi.de)
47 //////////////////////////////////////////////////////////////////////////////////////
52 #include <TProfile2D.h>
54 #include <TGraphErrors.h>
55 #include <TObjArray.h>
61 #include <TDirectory.h>
62 #include <TTreeStream.h>
63 #include <TLinearFitter.h>
70 #include "AliMathBase.h"
72 #include "AliTRDCalibraFit.h"
73 #include "AliTRDCalibraMode.h"
74 #include "AliTRDCalibraVector.h"
75 #include "AliTRDCalibraVdriftLinearFit.h"
76 #include "AliTRDCalibraExbAltFit.h"
77 #include "AliTRDcalibDB.h"
78 #include "AliTRDgeometry.h"
79 #include "AliTRDpadPlane.h"
80 #include "AliTRDgeometry.h"
81 #include "AliTRDCommonParam.h"
82 #include "./Cal/AliTRDCalROC.h"
83 #include "./Cal/AliTRDCalPad.h"
84 #include "./Cal/AliTRDCalDet.h"
87 ClassImp(AliTRDCalibraFit)
89 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
90 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
92 //_____________singleton implementation_________________________________________________
93 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
96 // Singleton implementation
99 if (fgTerminated != kFALSE) {
103 if (fgInstance == 0) {
104 fgInstance = new AliTRDCalibraFit();
110 //______________________________________________________________________________________
111 void AliTRDCalibraFit::Terminate()
114 // Singleton implementation
115 // Deletes the instance of this class
118 fgTerminated = kTRUE;
120 if (fgInstance != 0) {
126 //______________________________________________________________________________________
127 AliTRDCalibraFit::AliTRDCalibraFit()
130 ,fNumberOfBinsExpected(0)
132 ,fBeginFitCharge(3.5)
133 ,fOutliersFitChargeLow(0.03)
134 ,fOutliersFitChargeHigh(0.80)
136 ,fTakeTheMaxPH(kTRUE)
145 ,fNumberFitSuccess(0)
152 ,fCalibraMode(new AliTRDCalibraMode())
157 ,fScaleFitFactor(0.0)
166 ,fCalDetVdriftUsed(0x0)
168 ,fCurrentCoefDetector(0x0)
169 ,fCurrentCoefDetector2(0x0)
174 // Default constructor
177 fGeo = new AliTRDgeometry();
179 // Current variables initialised
180 for (Int_t k = 0; k < 2; k++) {
181 fCurrentCoef[k] = 0.0;
182 fCurrentCoef2[k] = 0.0;
184 for (Int_t i = 0; i < 3; i++) {
190 //______________________________________________________________________________________
191 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
194 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
196 ,fBeginFitCharge(c.fBeginFitCharge)
197 ,fOutliersFitChargeLow(c.fOutliersFitChargeLow)
198 ,fOutliersFitChargeHigh(c.fOutliersFitChargeHigh)
199 ,fFitPHPeriode(c.fFitPHPeriode)
200 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
201 ,fT0Shift0(c.fT0Shift0)
202 ,fT0Shift1(c.fT0Shift1)
203 ,fRangeFitPRF(c.fRangeFitPRF)
205 ,fMinEntries(c.fMinEntries)
207 ,fScaleGain(c.fScaleGain)
208 ,fNumberFit(c.fNumberFit)
209 ,fNumberFitSuccess(c.fNumberFitSuccess)
210 ,fNumberEnt(c.fNumberEnt)
211 ,fStatisticMean(c.fStatisticMean)
213 ,fDebugLevel(c.fDebugLevel)
214 ,fFitVoir(c.fFitVoir)
215 ,fMagneticField(c.fMagneticField)
217 ,fCurrentCoefE(c.fCurrentCoefE)
218 ,fCurrentCoefE2(c.fCurrentCoefE2)
221 ,fScaleFitFactor(c.fScaleFitFactor)
222 ,fEntriesCurrent(c.fEntriesCurrent)
223 ,fCountDet(c.fCountDet)
230 ,fCalDetVdriftUsed(0x0)
232 ,fCurrentCoefDetector(0x0)
233 ,fCurrentCoefDetector2(0x0)
241 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
243 //Current variables initialised
244 for (Int_t k = 0; k < 2; k++) {
245 fCurrentCoef[k] = 0.0;
246 fCurrentCoef2[k] = 0.0;
248 for (Int_t i = 0; i < 3; i++) {
252 if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet);
253 if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
255 if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
256 if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
258 if(c.fCalDetVdriftUsed) fCalDetVdriftUsed = new AliTRDCalDet(*c.fCalDetVdriftUsed);
259 if(c.fCalDetExBUsed) fCalDetExBUsed = new AliTRDCalDet(*c.fCalDetExBUsed);
261 fVectorFit.SetName(c.fVectorFit.GetName());
262 for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
263 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
264 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
266 if (GetStack(detector) == 2) {
272 Float_t *coef = new Float_t[ntotal];
273 for (Int_t i = 0; i < ntotal; i++) {
274 coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
276 fitInfo->SetCoef(coef);
277 fitInfo->SetDetector(detector);
278 fVectorFit.Add((TObject *) fitInfo);
280 fVectorFit.SetName(c.fVectorFit.GetName());
281 for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
282 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
283 Int_t detector = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
285 if (GetStack(detector) == 2) {
291 Float_t *coef = new Float_t[ntotal];
292 for (Int_t i = 0; i < ntotal; i++) {
293 coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
295 fitInfo->SetCoef(coef);
296 fitInfo->SetDetector(detector);
297 fVectorFit2.Add((TObject *) fitInfo);
302 fGeo = new AliTRDgeometry();
305 //____________________________________________________________________________________
306 AliTRDCalibraFit::~AliTRDCalibraFit()
309 // AliTRDCalibraFit destructor
311 if ( fDebugStreamer ) delete fDebugStreamer;
312 if ( fCalDet ) delete fCalDet;
313 if ( fCalDet2 ) delete fCalDet2;
314 if ( fCalROC ) delete fCalROC;
315 if ( fCalROC2 ) delete fCalROC2;
316 if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed;
317 if ( fCalDetExBUsed) delete fCalDetExBUsed;
318 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
319 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
321 fVectorFit2.Delete();
327 //_____________________________________________________________________________
328 void AliTRDCalibraFit::Destroy()
340 //_____________________________________________________________________________
341 void AliTRDCalibraFit::DestroyDebugStreamer()
344 // Delete DebugStreamer
347 if ( fDebugStreamer ) delete fDebugStreamer;
348 fDebugStreamer = 0x0;
351 //____________Functions fit Online CH2d________________________________________
352 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
355 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
356 // calibration group normalized the resulted coefficients (to 1 normally)
359 // Set the calibration mode
360 //const char *name = ch->GetTitle();
361 TString name = ch->GetTitle();
362 if(!SetModeCalibration(name,0)) return kFALSE;
364 // Number of Ybins (detectors or groups of pads)
365 Int_t nbins = ch->GetNbinsX();// charge
366 Int_t nybins = ch->GetNbinsY();// groups number
367 if (!InitFit(nybins,0)) {
373 fStatisticMean = 0.0;
375 fNumberFitSuccess = 0;
377 // Init fCountDet and fCount
378 InitfCountDetAndfCount(0);
379 // Beginning of the loop betwwen dect1 and dect2
380 for (Int_t idect = fDect1; idect < fDect2; idect++) {
381 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
382 UpdatefCountDetAndfCount(idect,0);
383 ReconstructFitRowMinRowMax(idect, 0);
385 TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
386 projch->SetDirectory(0);
387 // Number of entries for this calibration group
388 Double_t nentries = 0.0;
390 for (Int_t k = 0; k < nbins; k++) {
391 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
392 nentries += ch->GetBinContent(binnb);
393 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
394 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
396 projch->SetEntries(nentries);
397 //printf("The number of entries for the group %d is %f\n",idect,nentries);
402 // Rebin and statistic stuff
404 projch = ReBin((TH1I *) projch);
406 // This detector has not enough statistics or was off
407 if (nentries <= fMinEntries) {
408 NotEnoughStatisticCH(idect);
409 if (fDebugLevel != 1) {
414 // Statistics of the group fitted
415 fStatisticMean += nentries;
420 case 0: FitMeanW((TH1 *) projch, nentries); break;
421 case 1: FitMean((TH1 *) projch, nentries, mean); break;
422 case 2: FitLandau((TH1 *) projch, mean, nentries); break;
423 case 3: FitCH((TH1 *) projch, mean, nentries); break;
424 case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
425 case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
426 default: return kFALSE;
429 FillInfosFitCH(idect);
431 if (fDebugLevel != 1) {
436 if (fDebugLevel != 1) {
440 if (fNumberFit > 0) {
441 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit), fNumberFitSuccess));
442 fStatisticMean = fStatisticMean / fNumberFit;
445 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
447 delete fDebugStreamer;
448 fDebugStreamer = 0x0;
452 //____________Functions fit Online CH2d________________________________________
453 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
456 // Reconstruct a 1D histo from the vectorCH for each calibration group,
457 // fit the histo, normalized the resulted coefficients (to 1 normally)
460 // Set the calibraMode
461 //const char *name = calvect->GetNameCH();
462 TString name = calvect->GetNameCH();
463 if(!SetModeCalibration(name,0)) return kFALSE;
465 // Number of Xbins (detectors or groups of pads)
466 if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
472 fStatisticMean = 0.0;
474 fNumberFitSuccess = 0;
476 // Init fCountDet and fCount
477 InitfCountDetAndfCount(0);
478 // Beginning of the loop between dect1 and dect2
479 for (Int_t idect = fDect1; idect < fDect2; idect++) {
480 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
481 UpdatefCountDetAndfCount(idect,0);
482 ReconstructFitRowMinRowMax(idect,0);
484 Double_t nentries = 0.0;
486 if(!calvect->GetCHEntries(fCountDet)) {
487 NotEnoughStatisticCH(idect);
493 TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
494 projch->SetDirectory(0);
495 for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
496 nentries += projch->GetBinContent(k+1);
497 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
503 //printf("The number of entries for the group %d is %f\n",idect,nentries);
506 projch = ReBin((TH1F *) projch);
508 // This detector has not enough statistics or was not found in VectorCH
509 if (nentries <= fMinEntries) {
510 NotEnoughStatisticCH(idect);
513 // Statistic of the histos fitted
514 fStatisticMean += nentries;
519 case 0: FitMeanW((TH1 *) projch, nentries); break;
520 case 1: FitMean((TH1 *) projch, nentries, mean); break;
521 case 2: FitLandau((TH1 *) projch, mean, nentries); break;
522 case 3: FitCH((TH1 *) projch, mean, nentries); break;
523 case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
524 case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
525 default: return kFALSE;
528 FillInfosFitCH(idect);
531 if (fDebugLevel != 1) {
535 if (fNumberFit > 0) {
536 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit), fNumberFitSuccess));
537 fStatisticMean = fStatisticMean / fNumberFit;
540 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
542 delete fDebugStreamer;
543 fDebugStreamer = 0x0;
546 //____________Functions fit Online CH2d________________________________________
547 Double_t AliTRDCalibraFit::AnalyseCHAllTogether(const TH2I *ch)
550 // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
551 // calibration group normalized the resulted coefficients (to 1 normally)
553 Int_t nbins = ch->GetNbinsX();// charge
554 Int_t nybins = ch->GetNbinsY();// groups number
556 TH1I *projch = (TH1I *) ch->ProjectionX("projch",1,nybins+1,(Option_t *)"e");
557 projch->SetDirectory(0);
558 // Number of entries for this calibration group
559 Double_t nentries = 0.0;
561 for (Int_t k = 0; k < nbins; k++) {
562 nentries += projch->GetBinContent(k+1);
563 mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
564 projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
566 projch->SetEntries(nentries);
567 //printf("The number of entries for the group %d is %f\n",idect,nentries);
571 // This detector has not enough statistics or was off
572 if (nentries <= fMinEntries) {
574 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
580 case 0: FitMeanW((TH1 *) projch, nentries); break;
581 case 1: FitMean((TH1 *) projch, nentries, mean); break;
582 case 2: FitLandau((TH1 *) projch, mean, nentries); break;
583 case 3: FitCH((TH1 *) projch, mean, nentries); break;
584 case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
585 case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
586 default: return -100.0;
588 delete fDebugStreamer;
589 fDebugStreamer = 0x0;
591 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
595 //________________functions fit Online PH2d____________________________________
596 Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph)
599 // Take the 1D profiles (average pulse height), projections of the 2D PH
600 // on the Xaxis, for each calibration group
601 // Reconstruct a drift velocity
602 // A first calibration of T0 is also made using the same method
605 // Number of Xbins (detectors or groups of pads)
606 Int_t nbins = ph->GetNbinsX();// time
607 Int_t nybins = ph->GetNbinsY();// calibration group
610 TH1D *projph = (TH1D *) ph->ProjectionX("projph",1,nybins+1,(Option_t *) "e");
611 projph->SetDirectory(0);
612 // Number of entries for this calibration group
613 Double_t nentries = 0;
614 for(Int_t idect = 0; idect < nybins; idect++){
615 for (Int_t k = 0; k < nbins; k++) {
616 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
617 nentries += ph->GetBinEntries(binnb);
620 //printf("AnalysePHAllTogether:: the number of entries is %f\n",nentries);
621 // This detector has not enough statistics or was off
622 if (nentries <= fMinEntries) {
623 AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
624 if (fDebugLevel != 1) {
630 //printf("Method\n");
633 case 0: FitLagrangePoly((TH1 *) projph); break;
634 case 1: FitPente((TH1 *) projph); break;
635 case 2: FitPH((TH1 *) projph,0); break;
636 default: return -100.0;
639 if (fDebugLevel != 1) {
642 delete fDebugStreamer;
643 fDebugStreamer = 0x0;
645 if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
649 //____________Functions fit Online PH2d________________________________________
650 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
653 // Reconstruct the average pulse height from the vectorPH for each
655 // Reconstruct a drift velocity
656 // A first calibration of T0 is also made using the same method (slope method)
659 // Set the calibration mode
660 //const char *name = calvect->GetNamePH();
661 TString name = calvect->GetNamePH();
662 if(!SetModeCalibration(name,1)) return kFALSE;
664 // Number of Xbins (detectors or groups of pads)
665 if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
671 fStatisticMean = 0.0;
673 fNumberFitSuccess = 0;
675 // Init fCountDet and fCount
676 InitfCountDetAndfCount(1);
677 // Beginning of the loop
678 for (Int_t idect = fDect1; idect < fDect2; idect++) {
679 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
680 UpdatefCountDetAndfCount(idect,1);
681 ReconstructFitRowMinRowMax(idect,1);
684 if(!calvect->GetPHEntries(fCountDet)) {
685 NotEnoughStatisticPH(idect,fEntriesCurrent);
690 TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
691 projph->SetDirectory(0);
692 if(fEntriesCurrent > 0) fNumberEnt++;
693 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
694 // This detector has not enough statistics or was off
695 if (fEntriesCurrent <= fMinEntries) {
696 //printf("Not enough stat!\n");
697 NotEnoughStatisticPH(idect,fEntriesCurrent);
700 // Statistic of the histos fitted
702 fStatisticMean += fEntriesCurrent;
703 // Calcul of "real" coef
704 CalculVdriftCoefMean();
709 case 0: FitLagrangePoly((TH1 *) projph); break;
710 case 1: FitPente((TH1 *) projph); break;
711 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
712 default: return kFALSE;
714 // Fill the tree if end of a detector or only the pointer to the branch!!!
715 FillInfosFitPH(idect,fEntriesCurrent);
719 if (fNumberFit > 0) {
720 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
721 fStatisticMean = fStatisticMean / fNumberFit;
724 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
726 delete fDebugStreamer;
727 fDebugStreamer = 0x0;
730 //________________functions fit Online PH2d____________________________________
731 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
734 // Take the 1D profiles (average pulse height), projections of the 2D PH
735 // on the Xaxis, for each calibration group
736 // Reconstruct a drift velocity
737 // A first calibration of T0 is also made using the same method
740 // Set the calibration mode
741 //const char *name = ph->GetTitle();
742 TString name = ph->GetTitle();
743 if(!SetModeCalibration(name,1)) return kFALSE;
745 //printf("Mode calibration set\n");
747 // Number of Xbins (detectors or groups of pads)
748 Int_t nbins = ph->GetNbinsX();// time
749 Int_t nybins = ph->GetNbinsY();// calibration group
750 if (!InitFit(nybins,1)) {
754 //printf("Init fit\n");
760 //printf("Init fit PH\n");
762 fStatisticMean = 0.0;
764 fNumberFitSuccess = 0;
766 // Init fCountDet and fCount
767 InitfCountDetAndfCount(1);
768 //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
770 // Beginning of the loop
771 for (Int_t idect = fDect1; idect < fDect2; idect++) {
772 //printf("idect = %d\n",idect);
773 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
774 UpdatefCountDetAndfCount(idect,1);
775 ReconstructFitRowMinRowMax(idect,1);
777 TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
778 projph->SetDirectory(0);
779 // Number of entries for this calibration group
780 Double_t nentries = 0;
781 for (Int_t k = 0; k < nbins; k++) {
782 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
783 nentries += ph->GetBinEntries(binnb);
788 //printf("The number of entries for the group %d is %f\n",idect,nentries);
789 // This detector has not enough statistics or was off
790 if (nentries <= fMinEntries) {
791 //printf("Not enough statistic!\n");
792 NotEnoughStatisticPH(idect,nentries);
793 if (fDebugLevel != 1) {
798 // Statistics of the histos fitted
800 fStatisticMean += nentries;
801 // Calcul of "real" coef
802 CalculVdriftCoefMean();
805 //printf("Method\n");
808 case 0: FitLagrangePoly((TH1 *) projph); break;
809 case 1: FitPente((TH1 *) projph); break;
810 case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
811 default: return kFALSE;
813 // Fill the tree if end of a detector or only the pointer to the branch!!!
814 FillInfosFitPH(idect,nentries);
816 if (fDebugLevel != 1) {
821 if (fNumberFit > 0) {
822 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
823 fStatisticMean = fStatisticMean / fNumberFit;
826 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
828 delete fDebugStreamer;
829 fDebugStreamer = 0x0;
832 //____________Functions fit Online PRF2d_______________________________________
833 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
836 // Take the 1D profiles (pad response function), projections of the 2D PRF
837 // on the Xaxis, for each calibration group
838 // Fit with a gaussian to reconstruct the sigma of the pad response function
841 // Set the calibration mode
842 //const char *name = prf->GetTitle();
843 TString name = prf->GetTitle();
844 if(!SetModeCalibration(name,2)) return kFALSE;
846 // Number of Ybins (detectors or groups of pads)
847 Int_t nybins = prf->GetNbinsY();// calibration groups
848 Int_t nbins = prf->GetNbinsX();// bins
849 Int_t nbg = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
850 if((nbg > 0) || (nbg == -1)) return kFALSE;
851 if (!InitFit(nybins,2)) {
857 fStatisticMean = 0.0;
859 fNumberFitSuccess = 0;
861 // Init fCountDet and fCount
862 InitfCountDetAndfCount(2);
863 // Beginning of the loop
864 for (Int_t idect = fDect1; idect < fDect2; idect++) {
865 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
866 UpdatefCountDetAndfCount(idect,2);
867 ReconstructFitRowMinRowMax(idect,2);
869 TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
870 projprf->SetDirectory(0);
871 // Number of entries for this calibration group
872 Double_t nentries = 0;
873 for (Int_t k = 0; k < nbins; k++) {
874 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
875 nentries += prf->GetBinEntries(binnb);
877 if(nentries > 0) fNumberEnt++;
878 // This detector has not enough statistics or was off
879 if (nentries <= fMinEntries) {
880 NotEnoughStatisticPRF(idect);
881 if (fDebugLevel != 1) {
886 // Statistics of the histos fitted
888 fStatisticMean += nentries;
889 // Calcul of "real" coef
894 case 0: FitPRF((TH1 *) projprf); break;
895 case 1: RmsPRF((TH1 *) projprf); break;
896 default: return kFALSE;
898 // Fill the tree if end of a detector or only the pointer to the branch!!!
899 FillInfosFitPRF(idect);
901 if (fDebugLevel != 1) {
906 if (fNumberFit > 0) {
907 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
908 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
909 AliInfo(Form("There is a mean statistic of: %f over these fitted histograms and %d successfulled fits"
910 ,(Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
911 fStatisticMean = fStatisticMean / fNumberFit;
914 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
916 delete fDebugStreamer;
917 fDebugStreamer = 0x0;
920 //____________Functions fit Online PRF2d_______________________________________
921 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
924 // Take the 1D profiles (pad response function), projections of the 2D PRF
925 // on the Xaxis, for each calibration group
926 // Fit with a gaussian to reconstruct the sigma of the pad response function
929 // Set the calibration mode
930 //const char *name = prf->GetTitle();
931 TString name = prf->GetTitle();
932 if(!SetModeCalibration(name,2)) return kFALSE;
934 // Number of Ybins (detectors or groups of pads)
935 const TAxis *xprf = prf->GetXaxis();
936 const TAxis *yprf = prf->GetYaxis();
937 Int_t nybins = yprf->GetNbins();// calibration groups
938 Int_t nbins = xprf->GetNbins();// bins
939 Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
940 Float_t upedge = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
941 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
942 if(nbg == -1) return kFALSE;
943 if(nbg > 0) fMethod = 1;
945 if (!InitFit(nybins,2)) {
951 fStatisticMean = 0.0;
953 fNumberFitSuccess = 0;
955 // Init fCountDet and fCount
956 InitfCountDetAndfCount(2);
957 // Beginning of the loop
958 for (Int_t idect = fDect1; idect < fDect2; idect++) {
959 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
960 UpdatefCountDetAndfCount(idect,2);
961 ReconstructFitRowMinRowMax(idect,2);
962 // Build the array of entries and sum
963 TArrayD arraye = TArrayD(nbins);
964 TArrayD arraym = TArrayD(nbins);
965 TArrayD arrayme = TArrayD(nbins);
966 Double_t nentries = 0;
967 //printf("nbins %d\n",nbins);
968 for (Int_t k = 0; k < nbins; k++) {
969 Int_t binnb = (nbins+2)*(idect+1)+(k+1);
970 Double_t entries = (Double_t)prf->GetBinEntries(binnb);
971 Double_t mean = (Double_t)prf->GetBinContent(binnb);
972 Double_t error = (Double_t)prf->GetBinError(binnb);
973 //printf("for %d we have %f\n",k,entries);
975 arraye.AddAt(entries,k);
976 arraym.AddAt(mean,k);
977 arrayme.AddAt(error,k);
979 if(nentries > 0) fNumberEnt++;
980 //printf("The number of entries for the group %d is %f\n",idect,nentries);
981 // This detector has not enough statistics or was off
982 if (nentries <= fMinEntries) {
983 NotEnoughStatisticPRF(idect);
986 // Statistics of the histos fitted
988 fStatisticMean += nentries;
989 // Calcul of "real" coef
994 case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
995 case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
996 default: return kFALSE;
998 // Fill the tree if end of a detector or only the pointer to the branch!!!
999 FillInfosFitPRF(idect);
1002 if (fNumberFit > 0) {
1003 AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
1004 AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
1005 AliInfo(Form("There is a mean statistic of: %f over these fitted histograms and %d successfulled fits"
1006 ,(Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
1007 fStatisticMean = fStatisticMean / fNumberFit;
1010 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1012 delete fDebugStreamer;
1013 fDebugStreamer = 0x0;
1016 //____________Functions fit Online PRF2d_______________________________________
1017 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
1020 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1021 // each calibration group
1022 // Fit with a gaussian to reconstruct the sigma of the pad response function
1025 // Set the calibra mode
1026 //const char *name = calvect->GetNamePRF();
1027 TString name = calvect->GetNamePRF();
1028 if(!SetModeCalibration(name,2)) return kFALSE;
1029 //printf("test0 %s\n",name);
1031 // Number of Xbins (detectors or groups of pads)
1032 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1033 //printf("test1\n");
1036 if (!InitFitPRF()) {
1037 ///printf("test2\n");
1040 fStatisticMean = 0.0;
1042 fNumberFitSuccess = 0;
1044 // Init fCountDet and fCount
1045 InitfCountDetAndfCount(2);
1046 // Beginning of the loop
1047 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1048 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
1049 UpdatefCountDetAndfCount(idect,2);
1050 ReconstructFitRowMinRowMax(idect,2);
1052 fEntriesCurrent = 0;
1053 if(!calvect->GetPRFEntries(fCountDet)) {
1054 NotEnoughStatisticPRF(idect);
1057 TString tname("PRF");
1059 TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
1060 projprf->SetDirectory(0);
1061 if(fEntriesCurrent > 0) fNumberEnt++;
1062 // This detector has not enough statistics or was off
1063 if (fEntriesCurrent <= fMinEntries) {
1064 NotEnoughStatisticPRF(idect);
1067 // Statistic of the histos fitted
1069 fStatisticMean += fEntriesCurrent;
1070 // Calcul of "real" coef
1071 CalculPRFCoefMean();
1075 case 1: FitPRF((TH1 *) projprf); break;
1076 case 2: RmsPRF((TH1 *) projprf); break;
1077 default: return kFALSE;
1079 // Fill the tree if end of a detector or only the pointer to the branch!!!
1080 FillInfosFitPRF(idect);
1083 if (fNumberFit > 0) {
1084 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
1087 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1089 delete fDebugStreamer;
1090 fDebugStreamer = 0x0;
1093 //____________Functions fit Online PRF2d_______________________________________
1094 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
1097 // Reconstruct the 1D histo (pad response function) from the vectorPRD for
1098 // each calibration group
1099 // Fit with a gaussian to reconstruct the sigma of the pad response function
1102 // Set the calibra mode
1103 //const char *name = calvect->GetNamePRF();
1104 TString name = calvect->GetNamePRF();
1105 if(!SetModeCalibration(name,2)) return kFALSE;
1106 //printf("test0 %s\n",name);
1107 Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1108 //printf("test1 %d\n",nbg);
1109 if(nbg == -1) return kFALSE;
1110 if(nbg > 0) fMethod = 1;
1112 // Number of Xbins (detectors or groups of pads)
1113 if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
1114 //printf("test2\n");
1117 if (!InitFitPRF()) {
1118 //printf("test3\n");
1121 fStatisticMean = 0.0;
1123 fNumberFitSuccess = 0;
1127 Double_t *arrayx = 0;
1128 Double_t *arraye = 0;
1129 Double_t *arraym = 0;
1130 Double_t *arrayme = 0;
1131 Float_t lowedge = 0.0;
1132 Float_t upedge = 0.0;
1133 // Init fCountDet and fCount
1134 InitfCountDetAndfCount(2);
1135 // Beginning of the loop
1136 for (Int_t idect = fDect1; idect < fDect2; idect++) {
1137 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1138 UpdatefCountDetAndfCount(idect,2);
1139 ReconstructFitRowMinRowMax(idect,2);
1141 fEntriesCurrent = 0;
1142 if(!calvect->GetPRFEntries(fCountDet)) {
1143 NotEnoughStatisticPRF(idect);
1146 TString tname("PRF");
1148 TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1149 nbins = projprftree->GetN();
1150 arrayx = (Double_t *)projprftree->GetX();
1151 arraye = (Double_t *)projprftree->GetEX();
1152 arraym = (Double_t *)projprftree->GetY();
1153 arrayme = (Double_t *)projprftree->GetEY();
1154 Float_t step = arrayx[1]-arrayx[0];
1155 lowedge = arrayx[0] - step/2.0;
1156 upedge = arrayx[(nbins-1)] + step/2.0;
1157 //printf("nbins est %d\n",nbins);
1158 for(Int_t k = 0; k < nbins; k++){
1159 fEntriesCurrent += (Int_t)arraye[k];
1160 //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1161 if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1163 if(fEntriesCurrent > 0) fNumberEnt++;
1164 //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1165 // This detector has not enough statistics or was off
1166 if (fEntriesCurrent <= fMinEntries) {
1167 NotEnoughStatisticPRF(idect);
1170 // Statistic of the histos fitted
1172 fStatisticMean += fEntriesCurrent;
1173 // Calcul of "real" coef
1174 CalculPRFCoefMean();
1178 case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1179 case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1180 default: return kFALSE;
1182 // Fill the tree if end of a detector or only the pointer to the branch!!!
1183 FillInfosFitPRF(idect);
1186 if (fNumberFit > 0) {
1187 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
1190 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1192 delete fDebugStreamer;
1193 fDebugStreamer = 0x0;
1196 //____________Functions fit Online CH2d________________________________________
1197 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1200 // The linear method
1203 fStatisticMean = 0.0;
1205 fNumberFitSuccess = 0;
1207 if(!InitFitLinearFitter()) return kFALSE;
1210 for(Int_t idet = 0; idet < 540; idet++){
1213 //printf("detector number %d\n",idet);
1218 Double_t entriesCurrent = 0;
1220 Bool_t here = calivdli->GetParam(idet,¶m);
1221 Bool_t heree = calivdli->GetError(idet,&error);
1222 //printf("here %d and heree %d\n",here, heree);
1224 entriesCurrent = error[2];
1227 //printf("Number of entries %d\n",fEntriesCurrent);
1228 // Nothing found or not enough statistic
1229 if((!heree) || (!here) || (entriesCurrent <= fMinEntries)) {
1230 NotEnoughStatisticLinearFitter();
1237 fStatisticMean += entriesCurrent;
1240 if((-(param[1])) <= 0.000001) {
1241 NotEnoughStatisticLinearFitter();
1245 // CalculDatabaseVdriftandTan
1246 CalculVdriftLorentzCoef();
1247 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFit detector %d, vdrift %f and %f and exB %f and %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCurrentCoef[1],fCalDetExBUsed->GetValue(idet),fCurrentCoef2[1]);
1250 fNumberFitSuccess ++;
1252 // Put the fCurrentCoef
1253 fCurrentCoef[0] = -param[1];
1254 // here the database must be the one of the reconstruction for the lorentz angle....
1255 fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1256 fCurrentCoefE = error[1];
1257 fCurrentCoefE2 = error[0];
1258 if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
1259 fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1263 FillInfosFitLinearFitter();
1268 if (fNumberFit > 0) {
1269 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
1272 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1274 delete fDebugStreamer;
1275 fDebugStreamer = 0x0;
1279 //______________________________________________________________________________________
1280 Bool_t AliTRDCalibraFit::AnalyseExbAltFit(AliTRDCalibraExbAltFit *calivdli)
1283 // The linear method
1286 fStatisticMean = 0.0;
1288 fNumberFitSuccess = 0;
1290 if(!InitFitExbAlt()) return kFALSE;
1293 for(Int_t idet = 0; idet < 540; idet++){
1296 //printf("detector number %d\n",idet);
1301 Double_t entriesCurrent = 0;
1303 Bool_t here = calivdli->GetParam(idet,¶m);
1304 Bool_t heree = calivdli->GetError(idet,&error);
1305 //printf("here %d and heree %d\n",here, heree);
1307 entriesCurrent = error[2];
1310 //printf("Number of entries %d\n",fEntriesCurrent);
1311 // Nothing found or not enough statistic
1312 if((!heree) || (!here) || (entriesCurrent <= fMinEntries)) {
1313 NotEnoughStatisticExbAlt();
1320 fStatisticMean += entriesCurrent;
1323 fNumberFitSuccess ++;
1325 // Put the fCurrentCoef
1326 if(TMath::Abs(param[2])>0.0001){
1327 fCurrentCoef2[0] = -param[1]/2/param[2];
1328 fCurrentCoefE2 = 0;//error[1];
1330 fCurrentCoef2[0] = 100;
1331 fCurrentCoefE2 = 0;//error[1];
1335 FillInfosFitExbAlt();
1339 if (fNumberFit > 0) {
1340 AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
1343 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1345 delete fDebugStreamer;
1346 fDebugStreamer = 0x0;
1350 //____________Functions fit Online CH2d________________________________________
1351 void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall)
1354 // The linear method
1357 // Get the mean vdrift and exb used
1358 Double_t meanvdriftused = 0.0;
1359 Double_t meanexbused = 0.0;
1360 Double_t counterdet = 0.0;
1361 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) {
1362 vdriftoverall = -100.0;
1369 TH2S *linearfitterhisto = 0x0;
1371 for(Int_t idet = 0; idet < 540; idet++){
1373 TH2S * u = calivdli->GetLinearFitterHistoForce(idet);
1374 Double_t detectorentries = u->Integral();
1375 meanvdriftused += fCalDetVdriftUsed->GetValue(idet)*detectorentries;
1376 meanexbused += fCalDetExBUsed->GetValue(idet)*detectorentries;
1377 counterdet += detectorentries;
1379 //printf("detectorentries %f\n",detectorentries);
1381 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether detector %d, vdrift %f and exB %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCalDetExBUsed->GetValue(idet));
1383 if(idet == 0) linearfitterhisto = u;
1384 else linearfitterhisto->Add(u);
1387 if(counterdet > 0.0){
1388 meanvdriftused = meanvdriftused/counterdet;
1389 meanexbused = meanexbused/counterdet;
1392 vdriftoverall = -100.0;
1398 //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether MEAN vdrift %f and exB %f\n",meanvdriftused,meanexbused);
1402 Double_t entries = 0;
1403 TAxis *xaxis = linearfitterhisto->GetXaxis();
1404 TAxis *yaxis = linearfitterhisto->GetYaxis();
1405 TLinearFitter linearfitter = TLinearFitter(2,"pol1");
1407 Double_t integral = linearfitterhisto->Integral();
1408 //printf("Integral is %f\n",integral);
1409 Bool_t securitybreaking = kFALSE;
1410 if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
1411 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
1412 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
1413 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
1414 Double_t x = xaxis->GetBinCenter(ibinx+1);
1415 Double_t y = yaxis->GetBinCenter(ibiny+1);
1417 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
1418 if(!securitybreaking){
1419 linearfitter.AddPoint(&x,y);
1420 entries = entries+1.;
1423 if(entries< 1198.0){
1424 linearfitter.AddPoint(&x,y);
1425 entries = entries + 1.;
1434 //printf("AnalyseLinearFittersAllTogether::Find %d entries\n",entries);
1435 //printf("Minstats %d\n",fMinEntries);
1439 // Eval the linear fitter
1440 if(entries > fMinEntries){
1441 TVectorD par = TVectorD(2);
1443 if((linearfitter.EvalRobust(0.8)==0)) {
1444 //printf("Take the param\n");
1445 linearfitter.GetParameters(par);
1448 //printf("Finish\n");
1449 // Put the fCurrentCoef
1450 fCurrentCoef[0] = -par[1];
1451 // here the database must be the one of the reconstruction for the lorentz angle....
1452 if(fCurrentCoef[0] > 0.00001) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
1453 else fCurrentCoef2[0] = 100.0;
1458 fCurrentCoef[0] = -100.0;
1459 fCurrentCoef2[0] = 100.0;
1467 fCurrentCoef[0] = -100.0;
1468 fCurrentCoef2[0] = 100.0;
1472 vdriftoverall = fCurrentCoef[0];
1473 exboverall = fCurrentCoef2[0];
1476 delete linearfitterhisto;
1477 delete fDebugStreamer;
1478 fDebugStreamer = 0x0;
1481 //____________Functions for seeing if the pad is really okey___________________
1482 //_____________________________________________________________________________
1483 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
1486 // Get numberofgroupsprf
1490 const Char_t *pattern0 = "Ngp0";
1491 const Char_t *pattern1 = "Ngp1";
1492 const Char_t *pattern2 = "Ngp2";
1493 const Char_t *pattern3 = "Ngp3";
1494 const Char_t *pattern4 = "Ngp4";
1495 const Char_t *pattern5 = "Ngp5";
1496 const Char_t *pattern6 = "Ngp6";
1499 if (strstr(nametitle.Data(),pattern0)) {
1502 if (strstr(nametitle.Data(),pattern1)) {
1505 if (strstr(nametitle.Data(),pattern2)) {
1508 if (strstr(nametitle.Data(),pattern3)) {
1511 if (strstr(nametitle.Data(),pattern4)) {
1514 if (strstr(nametitle.Data(),pattern5)) {
1517 if (strstr(nametitle.Data(),pattern6)){
1524 //_____________________________________________________________________________
1525 Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
1528 // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1529 // corresponding to the given name
1532 if(!SetNzFromTObject(name,i)) return kFALSE;
1533 if(!SetNrphiFromTObject(name,i)) return kFALSE;
1538 //_____________________________________________________________________________
1539 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
1542 // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1543 // corresponding to the given TObject
1547 const Char_t *patternrphi0 = "Nrphi0";
1548 const Char_t *patternrphi1 = "Nrphi1";
1549 const Char_t *patternrphi2 = "Nrphi2";
1550 const Char_t *patternrphi3 = "Nrphi3";
1551 const Char_t *patternrphi4 = "Nrphi4";
1552 const Char_t *patternrphi5 = "Nrphi5";
1553 const Char_t *patternrphi6 = "Nrphi6";
1556 const Char_t *patternrphi10 = "Nrphi10";
1557 const Char_t *patternrphi100 = "Nrphi100";
1558 const Char_t *patternz10 = "Nz10";
1559 const Char_t *patternz100 = "Nz100";
1562 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1563 fCalibraMode->SetAllTogether(i);
1565 if (fDebugLevel > 1) {
1566 AliInfo(Form("fNbDet %d and 100",fNbDet));
1570 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1571 fCalibraMode->SetPerSuperModule(i);
1573 if (fDebugLevel > 1) {
1574 AliInfo(Form("fNDet %d and 100",fNbDet));
1579 if (strstr(name.Data(),patternrphi0)) {
1580 fCalibraMode->SetNrphi(i ,0);
1581 if (fDebugLevel > 1) {
1582 AliInfo(Form("fNbDet %d and 0",fNbDet));
1586 if (strstr(name.Data(),patternrphi1)) {
1587 fCalibraMode->SetNrphi(i, 1);
1588 if (fDebugLevel > 1) {
1589 AliInfo(Form("fNbDet %d and 1",fNbDet));
1593 if (strstr(name.Data(),patternrphi2)) {
1594 fCalibraMode->SetNrphi(i, 2);
1595 if (fDebugLevel > 1) {
1596 AliInfo(Form("fNbDet %d and 2",fNbDet));
1600 if (strstr(name.Data(),patternrphi3)) {
1601 fCalibraMode->SetNrphi(i, 3);
1602 if (fDebugLevel > 1) {
1603 AliInfo(Form("fNbDet %d and 3",fNbDet));
1607 if (strstr(name.Data(),patternrphi4)) {
1608 fCalibraMode->SetNrphi(i, 4);
1609 if (fDebugLevel > 1) {
1610 AliInfo(Form("fNbDet %d and 4",fNbDet));
1614 if (strstr(name.Data(),patternrphi5)) {
1615 fCalibraMode->SetNrphi(i, 5);
1616 if (fDebugLevel > 1) {
1617 AliInfo(Form("fNbDet %d and 5",fNbDet));
1621 if (strstr(name.Data(),patternrphi6)) {
1622 fCalibraMode->SetNrphi(i, 6);
1623 if (fDebugLevel > 1) {
1624 AliInfo(Form("fNbDet %d and 6",fNbDet));
1629 if (fDebugLevel > 1) {
1630 AliInfo(Form("fNbDet %d and rest",fNbDet));
1632 fCalibraMode->SetNrphi(i ,0);
1636 //_____________________________________________________________________________
1637 Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
1640 // Set fNz[i] of the AliTRDCalibraFit::Instance()
1641 // corresponding to the given TObject
1645 const Char_t *patternz0 = "Nz0";
1646 const Char_t *patternz1 = "Nz1";
1647 const Char_t *patternz2 = "Nz2";
1648 const Char_t *patternz3 = "Nz3";
1649 const Char_t *patternz4 = "Nz4";
1651 const Char_t *patternrphi10 = "Nrphi10";
1652 const Char_t *patternrphi100 = "Nrphi100";
1653 const Char_t *patternz10 = "Nz10";
1654 const Char_t *patternz100 = "Nz100";
1656 if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
1657 fCalibraMode->SetAllTogether(i);
1659 if (fDebugLevel > 1) {
1660 AliInfo(Form("fNbDet %d and 100",fNbDet));
1664 if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
1665 fCalibraMode->SetPerSuperModule(i);
1667 if (fDebugLevel > 1) {
1668 AliInfo(Form("fNbDet %d and 10",fNbDet));
1672 if (strstr(name.Data(),patternz0)) {
1673 fCalibraMode->SetNz(i, 0);
1674 if (fDebugLevel > 1) {
1675 AliInfo(Form("fNbDet %d and 0",fNbDet));
1679 if (strstr(name.Data(),patternz1)) {
1680 fCalibraMode->SetNz(i ,1);
1681 if (fDebugLevel > 1) {
1682 AliInfo(Form("fNbDet %d and 1",fNbDet));
1686 if (strstr(name.Data(),patternz2)) {
1687 fCalibraMode->SetNz(i ,2);
1688 if (fDebugLevel > 1) {
1689 AliInfo(Form("fNbDet %d and 2",fNbDet));
1693 if (strstr(name.Data(),patternz3)) {
1694 fCalibraMode->SetNz(i ,3);
1695 if (fDebugLevel > 1) {
1696 AliInfo(Form("fNbDet %d and 3",fNbDet));
1700 if (strstr(name.Data(),patternz4)) {
1701 fCalibraMode->SetNz(i ,4);
1702 if (fDebugLevel > 1) {
1703 AliInfo(Form("fNbDet %d and 4",fNbDet));
1708 if (fDebugLevel > 1) {
1709 AliInfo(Form("fNbDet %d and rest",fNbDet));
1711 fCalibraMode->SetNz(i ,0);
1714 //______________________________________________________________________
1715 void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
1717 // Remove the results too far from the mean value and rms
1718 // type: 0 gain, 1 vdrift
1722 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1724 AliInfo("The Vector Fit is not complete!");
1727 Int_t detector = -1;
1729 //Int_t sector = -1;
1730 Float_t value = 0.0;
1732 /////////////////////////////////
1733 // Calculate the mean values
1734 ////////////////////////////////
1736 ////////////////////////
1737 Double_t meanAll = 0.0;
1738 Double_t rmsAll = 0.0;
1743 for (Int_t k = 0; k < loop; k++) {
1744 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1746 //sector = GetSector(detector);
1748 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1750 rmsAll += value*value;
1756 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1757 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1758 for (Int_t row = 0; row < rowMax; row++) {
1759 for (Int_t col = 0; col < colMax; col++) {
1760 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1762 rmsAll += value*value;
1772 meanAll = meanAll/countAll;
1773 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1775 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1776 /////////////////////////////////////////////////
1778 ////////////////////////////////////////////////
1779 Double_t defaultvalue = -1.0;
1780 if(type==1) defaultvalue = -1.5;
1781 for (Int_t k = 0; k < loop; k++) {
1782 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1784 //sector = GetSector(detector);
1785 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1786 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1787 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1789 // remove the results too far away
1790 for (Int_t row = 0; row < rowMax; row++) {
1791 for (Int_t col = 0; col < colMax; col++) {
1792 value = coef[(Int_t)(col*rowMax+row)];
1793 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1794 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1800 //______________________________________________________________________
1801 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1803 // Remove the results too far from the mean and rms
1807 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1809 AliInfo("The Vector Fit is not complete!");
1812 Int_t detector = -1;
1814 //Int_t sector = -1;
1815 Float_t value = 0.0;
1817 /////////////////////////////////
1818 // Calculate the mean values
1819 ////////////////////////////////
1821 ////////////////////////
1822 Double_t meanAll = 0.0;
1823 Double_t rmsAll = 0.0;
1828 for (Int_t k = 0; k < loop; k++) {
1829 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1831 //sector = GetSector(detector);
1833 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1836 rmsAll += value*value;
1841 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1842 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1843 for (Int_t row = 0; row < rowMax; row++) {
1844 for (Int_t col = 0; col < colMax; col++) {
1845 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1847 rmsAll += value*value;
1856 meanAll = meanAll/countAll;
1857 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1859 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1860 /////////////////////////////////////////////////
1862 ////////////////////////////////////////////////
1863 for (Int_t k = 0; k < loop; k++) {
1864 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1866 //sector = GetSector(detector);
1867 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1868 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1869 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1871 // remove the results too far away
1872 for (Int_t row = 0; row < rowMax; row++) {
1873 for (Int_t col = 0; col < colMax; col++) {
1874 value = coef[(Int_t)(col*rowMax+row)];
1875 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) {
1876 //printf("value outlier %f\n",value);
1877 coef[(Int_t)(col*rowMax+row)] = 100.0;
1883 //______________________________________________________________________
1884 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1886 // ofwhat is equaled to 0: mean value of all passing detectors
1887 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1890 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1892 AliInfo("The Vector Fit is not complete!");
1895 Int_t detector = -1;
1897 Float_t value = 0.0;
1899 /////////////////////////////////
1900 // Calculate the mean values
1901 ////////////////////////////////
1903 ////////////////////////
1904 Double_t meanAll = 0.0;
1905 Double_t meanSupermodule[18];
1906 Double_t meanDetector[540];
1907 Double_t rmsAll = 0.0;
1908 Double_t rmsSupermodule[18];
1909 Double_t rmsDetector[540];
1911 Int_t countSupermodule[18];
1912 Int_t countDetector[540];
1913 for(Int_t sm = 0; sm < 18; sm++){
1914 rmsSupermodule[sm] = 0.0;
1915 meanSupermodule[sm] = 0.0;
1916 countSupermodule[sm] = 0;
1918 for(Int_t det = 0; det < 540; det++){
1919 rmsDetector[det] = 0.0;
1920 meanDetector[det] = 0.0;
1921 countDetector[det] = 0;
1926 for (Int_t k = 0; k < loop; k++) {
1927 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1928 sector = GetSector(detector);
1930 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1932 rmsDetector[detector] += value*value;
1933 meanDetector[detector] += value;
1934 countDetector[detector]++;
1935 rmsSupermodule[sector] += value*value;
1936 meanSupermodule[sector] += value;
1937 countSupermodule[sector]++;
1938 rmsAll += value*value;
1944 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1945 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1946 for (Int_t row = 0; row < rowMax; row++) {
1947 for (Int_t col = 0; col < colMax; col++) {
1948 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1950 rmsDetector[detector] += value*value;
1951 meanDetector[detector] += value;
1952 countDetector[detector]++;
1953 rmsSupermodule[sector] += value*value;
1954 meanSupermodule[sector] += value;
1955 countSupermodule[sector]++;
1956 rmsAll += value*value;
1966 meanAll = meanAll/countAll;
1967 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1969 for(Int_t sm = 0; sm < 18; sm++){
1970 if(countSupermodule[sm] > 0) {
1971 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1972 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1975 for(Int_t det = 0; det < 540; det++){
1976 if(countDetector[det] > 0) {
1977 meanDetector[det] = meanDetector[det]/countDetector[det];
1978 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1981 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1982 ///////////////////////////////////////////////
1983 // Put the mean value for the no-fitted
1984 /////////////////////////////////////////////
1985 for (Int_t k = 0; k < loop; k++) {
1986 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1987 sector = GetSector(detector);
1988 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1989 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1990 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1992 for (Int_t row = 0; row < rowMax; row++) {
1993 for (Int_t col = 0; col < colMax; col++) {
1994 value = coef[(Int_t)(col*rowMax+row)];
1996 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1998 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1999 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
2000 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
2004 if(fDebugLevel > 1){
2006 if ( !fDebugStreamer ) {
2008 TDirectory *backup = gDirectory;
2009 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2010 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2013 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2015 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
2016 "detector="<<detector<<
2028 //______________________________________________________________________
2029 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
2031 // ofwhat is equaled to 0: mean value of all passing detectors
2032 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
2035 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
2037 AliInfo("The Vector Fit is not complete!");
2040 Int_t detector = -1;
2042 Float_t value = 0.0;
2044 /////////////////////////////////
2045 // Calculate the mean values
2046 ////////////////////////////////
2048 ////////////////////////
2049 Double_t meanAll = 0.0;
2050 Double_t rmsAll = 0.0;
2051 Double_t meanSupermodule[18];
2052 Double_t rmsSupermodule[18];
2053 Double_t meanDetector[540];
2054 Double_t rmsDetector[540];
2056 Int_t countSupermodule[18];
2057 Int_t countDetector[540];
2058 for(Int_t sm = 0; sm < 18; sm++){
2059 rmsSupermodule[sm] = 0.0;
2060 meanSupermodule[sm] = 0.0;
2061 countSupermodule[sm] = 0;
2063 for(Int_t det = 0; det < 540; det++){
2064 rmsDetector[det] = 0.0;
2065 meanDetector[det] = 0.0;
2066 countDetector[det] = 0;
2070 for (Int_t k = 0; k < loop; k++) {
2071 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2072 sector = GetSector(detector);
2074 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
2076 rmsDetector[detector] += value*value;
2077 meanDetector[detector] += value;
2078 countDetector[detector]++;
2079 rmsSupermodule[sector] += value*value;
2080 meanSupermodule[sector] += value;
2081 countSupermodule[sector]++;
2083 rmsAll += value*value;
2088 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2089 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2090 for (Int_t row = 0; row < rowMax; row++) {
2091 for (Int_t col = 0; col < colMax; col++) {
2092 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2094 rmsDetector[detector] += value*value;
2095 meanDetector[detector] += value;
2096 countDetector[detector]++;
2097 rmsSupermodule[sector] += value*value;
2098 meanSupermodule[sector] += value;
2099 countSupermodule[sector]++;
2100 rmsAll += value*value;
2110 meanAll = meanAll/countAll;
2111 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
2113 for(Int_t sm = 0; sm < 18; sm++){
2114 if(countSupermodule[sm] > 0) {
2115 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
2116 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
2119 for(Int_t det = 0; det < 540; det++){
2120 if(countDetector[det] > 0) {
2121 meanDetector[det] = meanDetector[det]/countDetector[det];
2122 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
2125 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
2126 ////////////////////////////////////////////
2127 // Put the mean value for the no-fitted
2128 /////////////////////////////////////////////
2129 for (Int_t k = 0; k < loop; k++) {
2130 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2131 sector = GetSector(detector);
2132 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2133 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2134 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
2136 for (Int_t row = 0; row < rowMax; row++) {
2137 for (Int_t col = 0; col < colMax; col++) {
2138 value = coef[(Int_t)(col*rowMax+row)];
2140 if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2142 if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
2143 else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
2144 else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2148 if(fDebugLevel > 1){
2150 if ( !fDebugStreamer ) {
2152 TDirectory *backup = gDirectory;
2153 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2154 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2157 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2159 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
2160 "detector="<<detector<<
2173 //_____________________________________________________________________________
2174 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
2177 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2178 // It takes the mean value of the coefficients per detector
2179 // This object has to be written in the database
2182 // Create the DetObject
2183 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
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;
2191 for (Int_t k = 0; k < loop; k++) {
2192 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2195 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
2199 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2200 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2201 for (Int_t row = 0; row < rowMax; row++) {
2202 for (Int_t col = 0; col < colMax; col++) {
2203 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2204 mean += TMath::Abs(value);
2208 if(count > 0) mean = mean/count;
2210 object->SetValue(detector,mean);
2215 //_____________________________________________________________________________
2216 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
2219 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2220 // It takes the mean value of the coefficients per detector
2221 // This object has to be written in the database
2224 // Create the DetObject
2225 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2227 fScaleGain = scaleFitFactor;
2229 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2230 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2231 Int_t detector = -1;
2232 Float_t value = 0.0;
2234 for (Int_t k = 0; k < loop; k++) {
2235 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2238 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2239 if(!meanOtherBefore){
2240 if(value > 0) value = value*scaleFitFactor;
2242 else value = value*scaleFitFactor;
2243 mean = TMath::Abs(value);
2247 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2248 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2249 for (Int_t row = 0; row < rowMax; row++) {
2250 for (Int_t col = 0; col < colMax; col++) {
2251 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2252 if(!meanOtherBefore) {
2253 if(value > 0) value = value*scaleFitFactor;
2255 else value = value*scaleFitFactor;
2256 mean += TMath::Abs(value);
2260 if(count > 0) mean = mean/count;
2262 if(mean < 0.1) mean = 0.1;
2263 object->SetValue(detector,mean);
2268 //_____________________________________________________________________________
2269 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
2272 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2273 // It takes the min value of the coefficients per detector
2274 // This object has to be written in the database
2277 // Create the DetObject
2278 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2280 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2281 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2282 Int_t detector = -1;
2283 Float_t value = 0.0;
2285 for (Int_t k = 0; k < loop; k++) {
2286 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2287 Float_t min = 100.0;
2289 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2290 //printf("Create det object %f for %d\n",value,k);
2292 if(value > 70.0) value = value-100.0;
2297 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2298 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2299 for (Int_t row = 0; row < rowMax; row++) {
2300 for (Int_t col = 0; col < colMax; col++) {
2301 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2303 if(value > 70.0) value = value-100.0;
2305 if(min > value) min = value;
2309 object->SetValue(detector,min);
2315 //_____________________________________________________________________________
2316 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2319 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2320 // It takes the min value of the coefficients per detector
2321 // This object has to be written in the database
2324 // Create the DetObject
2325 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2328 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2329 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2330 Int_t detector = -1;
2331 Float_t value = 0.0;
2333 for (Int_t k = 0; k < loop; k++) {
2334 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2336 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2337 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2338 Float_t min = 100.0;
2339 for (Int_t row = 0; row < rowMax; row++) {
2340 for (Int_t col = 0; col < colMax; col++) {
2341 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2342 mean += -TMath::Abs(value);
2346 if(count > 0) mean = mean/count;
2348 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2349 if(value > 70.0) value = value-100.0;
2350 object->SetValue(detector,value);
2356 //_____________________________________________________________________________
2357 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit)
2360 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2361 // It takes the min value of the coefficients per detector
2362 // This object has to be written in the database
2365 // Create the DetObject
2366 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2369 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2370 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2371 Int_t detector = -1;
2372 Float_t value = 0.0;
2374 for (Int_t k = 0; k < loop; k++) {
2375 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2377 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2378 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2379 Float_t min = 100.0;
2380 for (Int_t row = 0; row < rowMax; row++) {
2381 for (Int_t col = 0; col < colMax; col++) {
2382 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2383 mean += -TMath::Abs(value);
2387 if(count > 0) mean = mean/count;
2389 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2390 //if(value > 70.0) value = value-100.0;
2391 object->SetValue(detector,value);
2397 //_____________________________________________________________________________
2398 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2401 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2402 // You need first to create the object for the detectors,
2403 // where the mean value is put.
2404 // This object has to be written in the database
2407 // Create the DetObject
2408 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2411 for(Int_t k = 0; k < 540; k++){
2412 AliTRDCalROC *calROC = object->GetCalROC(k);
2413 Int_t nchannels = calROC->GetNchannels();
2414 for(Int_t ch = 0; ch < nchannels; ch++){
2415 calROC->SetValue(ch,1.0);
2421 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2422 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2423 Int_t detector = -1;
2424 Float_t value = 0.0;
2426 for (Int_t k = 0; k < loop; k++) {
2427 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2428 AliTRDCalROC *calROC = object->GetCalROC(detector);
2429 Float_t mean = detobject->GetValue(detector);
2430 if(TMath::Abs(mean) <= 0.0000000001) continue;
2431 Int_t rowMax = calROC->GetNrows();
2432 Int_t colMax = calROC->GetNcols();
2433 for (Int_t row = 0; row < rowMax; row++) {
2434 for (Int_t col = 0; col < colMax; col++) {
2435 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2436 if(value > 0) value = value*scaleFitFactor;
2437 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2445 //_____________________________________________________________________________
2446 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2449 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2450 // You need first to create the object for the detectors,
2451 // where the mean value is put.
2452 // This object has to be written in the database
2455 // Create the DetObject
2456 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2459 for(Int_t k = 0; k < 540; k++){
2460 AliTRDCalROC *calROC = object->GetCalROC(k);
2461 Int_t nchannels = calROC->GetNchannels();
2462 for(Int_t ch = 0; ch < nchannels; ch++){
2463 calROC->SetValue(ch,1.0);
2469 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2470 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2471 Int_t detector = -1;
2472 Float_t value = 0.0;
2474 for (Int_t k = 0; k < loop; k++) {
2475 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2476 AliTRDCalROC *calROC = object->GetCalROC(detector);
2477 Float_t mean = detobject->GetValue(detector);
2478 if(mean == 0) continue;
2479 Int_t rowMax = calROC->GetNrows();
2480 Int_t colMax = calROC->GetNcols();
2481 for (Int_t row = 0; row < rowMax; row++) {
2482 for (Int_t col = 0; col < colMax; col++) {
2483 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2484 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2492 //_____________________________________________________________________________
2493 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2496 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2497 // You need first to create the object for the detectors,
2498 // where the mean value is put.
2499 // This object has to be written in the database
2502 // Create the DetObject
2503 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2506 for(Int_t k = 0; k < 540; k++){
2507 AliTRDCalROC *calROC = object->GetCalROC(k);
2508 Int_t nchannels = calROC->GetNchannels();
2509 for(Int_t ch = 0; ch < nchannels; ch++){
2510 calROC->SetValue(ch,0.0);
2516 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2517 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2518 Int_t detector = -1;
2519 Float_t value = 0.0;
2521 for (Int_t k = 0; k < loop; k++) {
2522 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2523 AliTRDCalROC *calROC = object->GetCalROC(detector);
2524 Float_t min = detobject->GetValue(detector);
2525 Int_t rowMax = calROC->GetNrows();
2526 Int_t colMax = calROC->GetNcols();
2527 for (Int_t row = 0; row < rowMax; row++) {
2528 for (Int_t col = 0; col < colMax; col++) {
2529 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2531 if(value > 70.0) value = value - 100.0;
2533 calROC->SetValue(col,row,value-min);
2541 //_____________________________________________________________________________
2542 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2545 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2546 // This object has to be written in the database
2549 // Create the DetObject
2550 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2552 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2553 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2554 Int_t detector = -1;
2555 Float_t value = 0.0;
2557 for (Int_t k = 0; k < loop; k++) {
2558 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2559 AliTRDCalROC *calROC = object->GetCalROC(detector);
2560 Int_t rowMax = calROC->GetNrows();
2561 Int_t colMax = calROC->GetNcols();
2562 for (Int_t row = 0; row < rowMax; row++) {
2563 for (Int_t col = 0; col < colMax; col++) {
2564 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2565 calROC->SetValue(col,row,TMath::Abs(value));
2573 //_____________________________________________________________________________
2574 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2577 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2578 // 0 successful fit 1 not successful fit
2579 // mean is the mean value over the successful fit
2580 // do not use it for t0: no meaning
2583 // Create the CalObject
2584 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2588 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2590 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2591 for(Int_t k = 0; k < 540; k++){
2592 object->SetValue(k,1.0);
2595 Int_t detector = -1;
2596 Float_t value = 0.0;
2598 for (Int_t k = 0; k < loop; k++) {
2599 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2600 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2601 if(value <= 0) object->SetValue(detector,1.0);
2603 object->SetValue(detector,0.0);
2608 if(count > 0) mean /= count;
2611 //_____________________________________________________________________________
2612 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2615 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2616 // 0 not successful fit 1 successful fit
2617 // mean mean value over the successful fit
2620 // Create the CalObject
2621 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2625 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2627 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2628 for(Int_t k = 0; k < 540; k++){
2629 AliTRDCalROC *calROC = object->GetCalROC(k);
2630 Int_t nchannels = calROC->GetNchannels();
2631 for(Int_t ch = 0; ch < nchannels; ch++){
2632 calROC->SetValue(ch,1.0);
2636 Int_t detector = -1;
2637 Float_t value = 0.0;
2639 for (Int_t k = 0; k < loop; k++) {
2640 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2641 AliTRDCalROC *calROC = object->GetCalROC(detector);
2642 Int_t nchannels = calROC->GetNchannels();
2643 for (Int_t ch = 0; ch < nchannels; ch++) {
2644 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2645 if(value <= 0) calROC->SetValue(ch,1.0);
2647 calROC->SetValue(ch,0.0);
2653 if(count > 0) mean /= count;
2656 //_____________________________________________________________________________
2657 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2660 // Set FitPH if 1 then each detector will be fitted
2663 if (periodeFitPH > 0) {
2664 fFitPHPeriode = periodeFitPH;
2667 AliInfo("periodeFitPH must be higher than 0!");
2671 //_____________________________________________________________________________
2672 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2675 // The fit of the deposited charge distribution begins at
2676 // histo->Mean()/beginFitCharge
2677 // You can here set beginFitCharge
2680 if (beginFitCharge > 0) {
2681 fBeginFitCharge = beginFitCharge;
2684 AliInfo("beginFitCharge must be strict positif!");
2689 //_____________________________________________________________________________
2690 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2693 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2694 // You can here set t0Shift0
2698 fT0Shift0 = t0Shift;
2701 AliInfo("t0Shift0 must be strict positif!");
2706 //_____________________________________________________________________________
2707 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2710 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2711 // You can here set t0Shift1
2715 fT0Shift1 = t0Shift;
2718 AliInfo("t0Shift must be strict positif!");
2723 //_____________________________________________________________________________
2724 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2727 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2728 // You can here set rangeFitPRF
2731 if ((rangeFitPRF > 0) &&
2732 (rangeFitPRF <= 1.5)) {
2733 fRangeFitPRF = rangeFitPRF;
2736 AliInfo("rangeFitPRF must be between 0 and 1.0");
2741 //_____________________________________________________________________________
2742 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2745 // Minimum entries for fitting
2748 if (minEntries > 0) {
2749 fMinEntries = minEntries;
2752 AliInfo("fMinEntries must be >= 0.");
2757 //_____________________________________________________________________________
2758 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2761 // Rebin with rebin time less bins the Ch histo
2762 // You can set here rebin that should divide the number of bins of CH histo
2767 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2770 AliInfo("You have to choose a positiv value!");
2774 //_____________________________________________________________________________
2775 Bool_t AliTRDCalibraFit::FillVectorFit()
2778 // For the Fit functions fill the vector Fit
2781 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2784 if (GetStack(fCountDet) == 2) {
2791 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2792 Float_t *coef = new Float_t[ntotal];
2793 for (Int_t i = 0; i < ntotal; i++) {
2794 coef[i] = fCurrentCoefDetector[i];
2797 Int_t detector = fCountDet;
2799 fitInfo->SetCoef(coef);
2800 fitInfo->SetDetector(detector);
2801 fVectorFit.Add((TObject *) fitInfo);
2806 //_____________________________________________________________________________
2807 Bool_t AliTRDCalibraFit::FillVectorFit2()
2810 // For the Fit functions fill the vector Fit
2813 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2816 if (GetStack(fCountDet) == 2) {
2823 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2824 Float_t *coef = new Float_t[ntotal];
2825 for (Int_t i = 0; i < ntotal; i++) {
2826 coef[i] = fCurrentCoefDetector2[i];
2829 Int_t detector = fCountDet;
2831 fitInfo->SetCoef(coef);
2832 fitInfo->SetDetector(detector);
2833 fVectorFit2.Add((TObject *) fitInfo);
2838 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2839 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2842 // Init the number of expected bins and fDect1[i] fDect2[i]
2845 gStyle->SetPalette(1);
2846 gStyle->SetOptStat(1111);
2847 gStyle->SetPadBorderMode(0);
2848 gStyle->SetCanvasColor(10);
2849 gStyle->SetPadLeftMargin(0.13);
2850 gStyle->SetPadRightMargin(0.01);
2852 // Mode groups of pads: the total number of bins!
2853 CalculNumberOfBinsExpected(i);
2855 // Quick verification that we have the good pad calibration mode!
2856 if (fNumberOfBinsExpected != nbins) {
2857 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2861 // Security for fDebug 3 and 4
2862 if ((fDebugLevel >= 3) &&
2866 AliInfo("This detector doesn't exit!");
2870 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2871 CalculDect1Dect2(i);
2876 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2877 Bool_t AliTRDCalibraFit::InitFitCH()
2880 // Init the fVectorFitCH for normalisation
2881 // Init the histo for debugging
2886 fScaleFitFactor = 0.0;
2887 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2888 fCurrentCoefDetector = new Float_t[2304];
2889 for (Int_t k = 0; k < 2304; k++) {
2890 fCurrentCoefDetector[k] = 0.0;
2892 fVectorFit.SetName("gainfactorscoefficients");
2894 // fDebug == 0 nothing
2895 // fDebug == 1 and fFitVoir no histo
2896 if (fDebugLevel == 1) {
2897 if(!CheckFitVoir()) return kFALSE;
2899 //Get the CalDet object
2901 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2903 AliInfo("Could not get calibDB");
2906 if(fCalDet) delete fCalDet;
2907 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2910 Float_t devalue = 1.0;
2911 if(fCalDet) delete fCalDet;
2912 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2913 for(Int_t k = 0; k < 540; k++){
2914 fCalDet->SetValue(k,devalue);
2920 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2921 Bool_t AliTRDCalibraFit::InitFitPH()
2924 // Init the arrays of results
2925 // Init the histos for debugging
2929 fVectorFit.SetName("driftvelocitycoefficients");
2930 fVectorFit2.SetName("t0coefficients");
2932 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2933 fCurrentCoefDetector = new Float_t[2304];
2934 for (Int_t k = 0; k < 2304; k++) {
2935 fCurrentCoefDetector[k] = 0.0;
2937 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
2938 fCurrentCoefDetector2 = new Float_t[2304];
2939 for (Int_t k = 0; k < 2304; k++) {
2940 fCurrentCoefDetector2[k] = 0.0;
2943 //fDebug == 0 nothing
2944 // fDebug == 1 and fFitVoir no histo
2945 if (fDebugLevel == 1) {
2946 if(!CheckFitVoir()) return kFALSE;
2948 //Get the CalDet object
2950 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2952 AliInfo("Could not get calibDB");
2955 if(fCalDet) delete fCalDet;
2956 if(fCalDet2) delete fCalDet2;
2957 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2958 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2961 Float_t devalue = 1.5;
2962 Float_t devalue2 = 0.0;
2963 if(fCalDet) delete fCalDet;
2964 if(fCalDet2) delete fCalDet2;
2965 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2966 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2967 for(Int_t k = 0; k < 540; k++){
2968 fCalDet->SetValue(k,devalue);
2969 fCalDet2->SetValue(k,devalue2);
2974 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2975 Bool_t AliTRDCalibraFit::InitFitPRF()
2978 // Init the calibration mode (Nz, Nrphi), the histograms for
2979 // debugging the fit methods if fDebug > 0,
2983 fVectorFit.SetName("prfwidthcoefficients");
2985 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2986 fCurrentCoefDetector = new Float_t[2304];
2987 for (Int_t k = 0; k < 2304; k++) {
2988 fCurrentCoefDetector[k] = 0.0;
2991 // fDebug == 0 nothing
2992 // fDebug == 1 and fFitVoir no histo
2993 if (fDebugLevel == 1) {
2994 if(!CheckFitVoir()) return kFALSE;
2998 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2999 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
3002 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3007 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
3008 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3009 fCurrentCoefDetector = new Float_t[2304];
3010 fCurrentCoefDetector2 = new Float_t[2304];
3011 for (Int_t k = 0; k < 2304; k++) {
3012 fCurrentCoefDetector[k] = 0.0;
3013 fCurrentCoefDetector2[k] = 0.0;
3016 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) return kFALSE;
3020 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3021 Bool_t AliTRDCalibraFit::InitFitExbAlt()
3024 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3029 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3030 fCurrentCoefDetector2 = new Float_t[2304];
3031 for (Int_t k = 0; k < 2304; k++) {
3032 fCurrentCoefDetector2[k] = 0.0;
3037 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3038 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
3041 // Init the current detector where we are fCountDet and the
3042 // next fCount for the functions Fit...
3045 // Loop on the Xbins of ch!!
3046 fCountDet = -1; // Current detector
3047 fCount = 0; // To find the next detector
3050 if (fDebugLevel >= 3) {
3051 // Set countdet to the detector
3052 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3053 // Set counter to write at the end of the detector
3055 // Get the right calib objects
3058 if(fDebugLevel == 1) {
3060 fCalibraMode->CalculXBins(fCountDet,i);
3061 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
3062 while(fCalibraMode->GetXbins(i) <=fFitVoir){
3064 fCalibraMode->CalculXBins(fCountDet,i);
3065 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
3071 fCount = fCalibraMode->GetXbins(i);
3073 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3074 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3075 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3076 ,(Int_t) GetStack(fCountDet)
3077 ,(Int_t) GetSector(fCountDet),i);
3080 //_______________________________________________________________________________
3081 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
3084 // Calculate the number of bins expected (calibration groups)
3087 fNumberOfBinsExpected = 0;
3089 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
3090 fNumberOfBinsExpected = 1;
3094 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
3095 fNumberOfBinsExpected = 18;
3099 fCalibraMode->ModePadCalibration(2,i);
3100 fCalibraMode->ModePadFragmentation(0,2,0,i);
3101 fCalibraMode->SetDetChamb2(i);
3102 if (fDebugLevel > 1) {
3103 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
3105 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
3106 fCalibraMode->ModePadCalibration(0,i);
3107 fCalibraMode->ModePadFragmentation(0,0,0,i);
3108 fCalibraMode->SetDetChamb0(i);
3109 if (fDebugLevel > 1) {
3110 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
3112 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
3115 //_______________________________________________________________________________
3116 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
3119 // Calculate the range of fits
3124 if (fDebugLevel == 1) {
3128 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
3130 fDect2 = fNumberOfBinsExpected;
3132 if (fDebugLevel >= 3) {
3133 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3134 fCalibraMode->CalculXBins(fCountDet,i);
3135 fDect1 = fCalibraMode->GetXbins(i);
3136 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3137 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3138 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3139 ,(Int_t) GetStack(fCountDet)
3140 ,(Int_t) GetSector(fCountDet),i);
3141 // Set for the next detector
3142 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3145 //_______________________________________________________________________________
3146 Bool_t AliTRDCalibraFit::CheckFitVoir()
3149 // Check if fFitVoir is in the range
3152 if (fFitVoir < fNumberOfBinsExpected) {
3153 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
3156 AliInfo("fFitVoir is out of range of the histo!");
3161 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3162 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
3165 // See if we are in a new detector and update the
3166 // variables fNfragZ and fNfragRphi if yes
3167 // Will never happen for only one detector (3 and 4)
3168 // Doesn't matter for 2
3170 if (fCount == idect) {
3171 // On en est au detector (or first detector in the group)
3173 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
3174 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3175 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3176 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3177 ,(Int_t) GetStack(fCountDet)
3178 ,(Int_t) GetSector(fCountDet),i);
3179 // Set for the next detector
3180 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3185 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3186 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
3189 // Reconstruct the min pad row, max pad row, min pad col and
3190 // max pad col of the calibration group for the Fit functions
3191 // idect is the calibration group inside the detector
3193 if (fDebugLevel != 1) {
3194 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
3196 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
3197 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
3199 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3200 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
3203 // For the case where there are not enough entries in the histograms
3204 // of the calibration group, the value present in the choosen database
3205 // will be put. A negativ sign enables to know that a fit was not possible.
3208 if (fDebugLevel == 1) {
3209 AliInfo("The element has not enough statistic to be fitted");
3211 else if (fNbDet > 0){
3212 Int_t firstdetector = fCountDet;
3213 Int_t lastdetector = fCountDet+fNbDet;
3214 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3215 // loop over detectors
3216 for(Int_t det = firstdetector; det < lastdetector; det++){
3218 //Set the calibration object again
3222 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3224 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3225 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3226 ,(Int_t) GetStack(fCountDet)
3227 ,(Int_t) GetSector(fCountDet),0);
3228 // Reconstruct row min row max
3229 ReconstructFitRowMinRowMax(idect,0);
3231 // Calcul the coef from the database choosen for the detector
3232 CalculChargeCoefMean(kFALSE);
3234 //stack 2, not stack 2
3236 if(GetStack(fCountDet) == 2) factor = 12;
3239 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3240 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3241 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3242 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3246 //Put default value negative
3247 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3248 fCurrentCoefE = 0.0;
3253 if(fDebugLevel > 1){
3255 if ( !fDebugStreamer ) {
3257 TDirectory *backup = gDirectory;
3258 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3259 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3262 Int_t detector = fCountDet;
3263 Int_t caligroup = idect;
3264 Short_t rowmin = fCalibraMode->GetRowMin(0);
3265 Short_t rowmax = fCalibraMode->GetRowMax(0);
3266 Short_t colmin = fCalibraMode->GetColMin(0);
3267 Short_t colmax = fCalibraMode->GetColMax(0);
3268 Float_t gf = fCurrentCoef[0];
3269 Float_t gfs = fCurrentCoef[1];
3270 Float_t gfE = fCurrentCoefE;
3272 (*fDebugStreamer) << "FillFillCH" <<
3273 "detector=" << detector <<
3274 "caligroup=" << caligroup <<
3275 "rowmin=" << rowmin <<
3276 "rowmax=" << rowmax <<
3277 "colmin=" << colmin <<
3278 "colmax=" << colmax <<
3286 for (Int_t k = 0; k < 2304; k++) {
3287 fCurrentCoefDetector[k] = 0.0;
3291 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3295 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
3297 // Calcul the coef from the database choosen
3298 CalculChargeCoefMean(kFALSE);
3300 //stack 2, not stack 2
3302 if(GetStack(fCountDet) == 2) factor = 12;
3305 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3306 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3307 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3308 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3312 //Put default value negative
3313 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3314 fCurrentCoefE = 0.0;
3323 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3324 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
3327 // For the case where there are not enough entries in the histograms
3328 // of the calibration group, the value present in the choosen database
3329 // will be put. A negativ sign enables to know that a fit was not possible.
3331 if (fDebugLevel == 1) {
3332 AliInfo("The element has not enough statistic to be fitted");
3334 else if (fNbDet > 0) {
3336 Int_t firstdetector = fCountDet;
3337 Int_t lastdetector = fCountDet+fNbDet;
3338 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3339 // loop over detectors
3340 for(Int_t det = firstdetector; det < lastdetector; det++){
3342 //Set the calibration object again
3346 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3348 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3349 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3350 ,(Int_t) GetStack(fCountDet)
3351 ,(Int_t) GetSector(fCountDet),1);
3352 // Reconstruct row min row max
3353 ReconstructFitRowMinRowMax(idect,1);
3355 // Calcul the coef from the database choosen for the detector
3356 CalculVdriftCoefMean();
3359 //stack 2, not stack 2
3361 if(GetStack(fCountDet) == 2) factor = 12;
3364 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3365 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3366 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3367 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3368 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3372 //Put default value negative
3373 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3374 fCurrentCoefE = 0.0;
3375 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3376 fCurrentCoefE2 = 0.0;
3382 if(fDebugLevel > 1){
3384 if ( !fDebugStreamer ) {
3386 TDirectory *backup = gDirectory;
3387 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3388 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3392 Int_t detector = fCountDet;
3393 Int_t caligroup = idect;
3394 Short_t rowmin = fCalibraMode->GetRowMin(1);
3395 Short_t rowmax = fCalibraMode->GetRowMax(1);
3396 Short_t colmin = fCalibraMode->GetColMin(1);
3397 Short_t colmax = fCalibraMode->GetColMax(1);
3398 Float_t vf = fCurrentCoef[0];
3399 Float_t vs = fCurrentCoef[1];
3400 Float_t vfE = fCurrentCoefE;
3401 Float_t t0f = fCurrentCoef2[0];
3402 Float_t t0s = fCurrentCoef2[1];
3403 Float_t t0E = fCurrentCoefE2;
3407 (* fDebugStreamer) << "FillFillPH"<<
3408 "detector="<<detector<<
3409 "nentries="<<nentries<<
3410 "caligroup="<<caligroup<<
3424 for (Int_t k = 0; k < 2304; k++) {
3425 fCurrentCoefDetector[k] = 0.0;
3426 fCurrentCoefDetector2[k] = 0.0;
3430 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3434 //AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3436 CalculVdriftCoefMean();
3439 //stack 2 and not stack 2
3441 if(GetStack(fCountDet) == 2) factor = 12;
3445 // Fill the fCurrentCoefDetector 2
3446 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3447 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3448 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3449 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3453 // Put the default value
3454 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3455 fCurrentCoefE = 0.0;
3456 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3457 fCurrentCoefE2 = 0.0;
3459 FillFillPH(idect,nentries);
3468 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3469 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3472 // For the case where there are not enough entries in the histograms
3473 // of the calibration group, the value present in the choosen database
3474 // will be put. A negativ sign enables to know that a fit was not possible.
3477 if (fDebugLevel == 1) {
3478 AliInfo("The element has not enough statistic to be fitted");
3480 else if (fNbDet > 0){
3482 Int_t firstdetector = fCountDet;
3483 Int_t lastdetector = fCountDet+fNbDet;
3484 // AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3486 // loop over detectors
3487 for(Int_t det = firstdetector; det < lastdetector; det++){
3489 //Set the calibration object again
3493 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3495 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3496 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3497 ,(Int_t) GetStack(fCountDet)
3498 ,(Int_t) GetSector(fCountDet),2);
3499 // Reconstruct row min row max
3500 ReconstructFitRowMinRowMax(idect,2);
3502 // Calcul the coef from the database choosen for the detector
3503 CalculPRFCoefMean();
3505 //stack 2, not stack 2
3507 if(GetStack(fCountDet) == 2) factor = 12;
3510 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3511 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3512 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3513 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3517 //Put default value negative
3518 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3519 fCurrentCoefE = 0.0;
3524 if(fDebugLevel > 1){
3526 if ( !fDebugStreamer ) {
3528 TDirectory *backup = gDirectory;
3529 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3530 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3533 Int_t detector = fCountDet;
3534 Int_t layer = GetLayer(fCountDet);
3535 Int_t caligroup = idect;
3536 Short_t rowmin = fCalibraMode->GetRowMin(2);
3537 Short_t rowmax = fCalibraMode->GetRowMax(2);
3538 Short_t colmin = fCalibraMode->GetColMin(2);
3539 Short_t colmax = fCalibraMode->GetColMax(2);
3540 Float_t widf = fCurrentCoef[0];
3541 Float_t wids = fCurrentCoef[1];
3542 Float_t widfE = fCurrentCoefE;
3544 (* fDebugStreamer) << "FillFillPRF"<<
3545 "detector="<<detector<<
3547 "caligroup="<<caligroup<<
3558 for (Int_t k = 0; k < 2304; k++) {
3559 fCurrentCoefDetector[k] = 0.0;
3563 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3567 // AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3569 CalculPRFCoefMean();
3571 // stack 2 and not stack 2
3573 if(GetStack(fCountDet) == 2) factor = 12;
3577 // Fill the fCurrentCoefDetector
3578 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3579 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3580 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3584 // Put the default value
3585 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3586 fCurrentCoefE = 0.0;
3594 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3595 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3598 // For the case where there are not enough entries in the histograms
3599 // of the calibration group, the value present in the choosen database
3600 // will be put. A negativ sign enables to know that a fit was not possible.
3603 // Calcul the coef from the database choosen
3604 CalculVdriftLorentzCoef();
3607 if(GetStack(fCountDet) == 2) factor = 1728;
3611 // Fill the fCurrentCoefDetector
3612 for (Int_t k = 0; k < factor; k++) {
3613 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3614 // should be negative
3615 fCurrentCoefDetector2[k] = fCurrentCoef2[1]+100.0;
3619 //Put default opposite sign only for vdrift
3620 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3621 fCurrentCoefE = 0.0;
3622 fCurrentCoef2[0] = fCurrentCoef2[1]+100.0;
3623 fCurrentCoefE2 = 0.0;
3625 FillFillLinearFitter();
3630 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3631 Bool_t AliTRDCalibraFit::NotEnoughStatisticExbAlt()
3634 // For the case where there are not enough entries in the histograms
3635 // of the calibration group, the value present in the choosen database
3636 // will be put. A negativ sign enables to know that a fit was not possible.
3640 if(GetStack(fCountDet) == 2) factor = 1728;
3644 // Fill the fCurrentCoefDetector
3645 for (Int_t k = 0; k < factor; k++) {
3646 fCurrentCoefDetector2[k] = 100.0;
3649 fCurrentCoef2[0] = 100.0;
3650 fCurrentCoefE2 = 0.0;
3657 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3658 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3661 // Fill the coefficients found with the fits or other
3662 // methods from the Fit functions
3665 if (fDebugLevel != 1) {
3667 Int_t firstdetector = fCountDet;
3668 Int_t lastdetector = fCountDet+fNbDet;
3669 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3670 // loop over detectors
3671 for(Int_t det = firstdetector; det < lastdetector; det++){
3673 //Set the calibration object again
3677 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3679 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3680 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3681 ,(Int_t) GetStack(fCountDet)
3682 ,(Int_t) GetSector(fCountDet),0);
3683 // Reconstruct row min row max
3684 ReconstructFitRowMinRowMax(idect,0);
3686 // Calcul the coef from the database choosen for the detector
3687 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3688 else CalculChargeCoefMean(kTRUE);
3690 //stack 2, not stack 2
3692 if(GetStack(fCountDet) == 2) factor = 12;
3695 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3696 Double_t coeftoput = 1.0;
3697 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3698 else coeftoput = fCurrentCoef[0];
3699 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3700 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3701 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3708 if(fDebugLevel > 1){
3710 if ( !fDebugStreamer ) {
3712 TDirectory *backup = gDirectory;
3713 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3714 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3717 Int_t detector = fCountDet;
3718 Int_t caligroup = idect;
3719 Short_t rowmin = fCalibraMode->GetRowMin(0);
3720 Short_t rowmax = fCalibraMode->GetRowMax(0);
3721 Short_t colmin = fCalibraMode->GetColMin(0);
3722 Short_t colmax = fCalibraMode->GetColMax(0);
3723 Float_t gf = fCurrentCoef[0];
3724 Float_t gfs = fCurrentCoef[1];
3725 Float_t gfE = fCurrentCoefE;
3727 (*fDebugStreamer) << "FillFillCH" <<
3728 "detector=" << detector <<
3729 "caligroup=" << caligroup <<
3730 "rowmin=" << rowmin <<
3731 "rowmax=" << rowmax <<
3732 "colmin=" << colmin <<
3733 "colmax=" << colmax <<
3741 for (Int_t k = 0; k < 2304; k++) {
3742 fCurrentCoefDetector[k] = 0.0;
3746 //printf("Check the count now: fCountDet %d\n",fCountDet);
3751 if(GetStack(fCountDet) == 2) factor = 12;
3754 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3755 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3756 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3767 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3768 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3771 // Fill the coefficients found with the fits or other
3772 // methods from the Fit functions
3775 if (fDebugLevel != 1) {
3778 Int_t firstdetector = fCountDet;
3779 Int_t lastdetector = fCountDet+fNbDet;
3780 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3782 // loop over detectors
3783 for(Int_t det = firstdetector; det < lastdetector; det++){
3785 //Set the calibration object again
3789 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3791 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3792 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3793 ,(Int_t) GetStack(fCountDet)
3794 ,(Int_t) GetSector(fCountDet),1);
3795 // Reconstruct row min row max
3796 ReconstructFitRowMinRowMax(idect,1);
3798 // Calcul the coef from the database choosen for the detector
3799 CalculVdriftCoefMean();
3802 //stack 2, not stack 2
3804 if(GetStack(fCountDet) == 2) factor = 12;
3807 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3808 Double_t coeftoput = 1.5;
3809 Double_t coeftoput2 = 0.0;
3811 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3812 else coeftoput = fCurrentCoef[0];
3814 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3815 else coeftoput2 = fCurrentCoef2[0];
3817 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3818 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3819 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3820 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3828 if(fDebugLevel > 1){
3830 if ( !fDebugStreamer ) {
3832 TDirectory *backup = gDirectory;
3833 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3834 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3838 Int_t detector = fCountDet;
3839 Int_t caligroup = idect;
3840 Short_t rowmin = fCalibraMode->GetRowMin(1);
3841 Short_t rowmax = fCalibraMode->GetRowMax(1);
3842 Short_t colmin = fCalibraMode->GetColMin(1);
3843 Short_t colmax = fCalibraMode->GetColMax(1);
3844 Float_t vf = fCurrentCoef[0];
3845 Float_t vs = fCurrentCoef[1];
3846 Float_t vfE = fCurrentCoefE;
3847 Float_t t0f = fCurrentCoef2[0];
3848 Float_t t0s = fCurrentCoef2[1];
3849 Float_t t0E = fCurrentCoefE2;
3853 (* fDebugStreamer) << "FillFillPH"<<
3854 "detector="<<detector<<
3855 "nentries="<<nentries<<
3856 "caligroup="<<caligroup<<
3870 for (Int_t k = 0; k < 2304; k++) {
3871 fCurrentCoefDetector[k] = 0.0;
3872 fCurrentCoefDetector2[k] = 0.0;
3876 //printf("Check the count now: fCountDet %d\n",fCountDet);
3881 if(GetStack(fCountDet) == 2) factor = 12;
3884 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3885 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3886 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3887 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3891 FillFillPH(idect,nentries);
3896 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3897 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3900 // Fill the coefficients found with the fits or other
3901 // methods from the Fit functions
3904 if (fDebugLevel != 1) {
3907 Int_t firstdetector = fCountDet;
3908 Int_t lastdetector = fCountDet+fNbDet;
3909 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3911 // loop over detectors
3912 for(Int_t det = firstdetector; det < lastdetector; det++){
3914 //Set the calibration object again
3918 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3920 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3921 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3922 ,(Int_t) GetStack(fCountDet)
3923 ,(Int_t) GetSector(fCountDet),2);
3924 // Reconstruct row min row max
3925 ReconstructFitRowMinRowMax(idect,2);
3927 // Calcul the coef from the database choosen for the detector
3928 CalculPRFCoefMean();
3930 //stack 2, not stack 2
3932 if(GetStack(fCountDet) == 2) factor = 12;
3935 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3936 Double_t coeftoput = 1.0;
3937 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3938 else coeftoput = fCurrentCoef[0];
3939 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3940 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3941 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3948 if(fDebugLevel > 1){
3950 if ( !fDebugStreamer ) {
3952 TDirectory *backup = gDirectory;
3953 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3954 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3957 Int_t detector = fCountDet;
3958 Int_t layer = GetLayer(fCountDet);
3959 Int_t caligroup = idect;
3960 Short_t rowmin = fCalibraMode->GetRowMin(2);
3961 Short_t rowmax = fCalibraMode->GetRowMax(2);
3962 Short_t colmin = fCalibraMode->GetColMin(2);
3963 Short_t colmax = fCalibraMode->GetColMax(2);
3964 Float_t widf = fCurrentCoef[0];
3965 Float_t wids = fCurrentCoef[1];
3966 Float_t widfE = fCurrentCoefE;
3968 (* fDebugStreamer) << "FillFillPRF"<<
3969 "detector="<<detector<<
3971 "caligroup="<<caligroup<<
3982 for (Int_t k = 0; k < 2304; k++) {
3983 fCurrentCoefDetector[k] = 0.0;
3987 //printf("Check the count now: fCountDet %d\n",fCountDet);
3992 if(GetStack(fCountDet) == 2) factor = 12;
3995 // Pointer to the branch
3996 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3997 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3998 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
4008 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4009 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
4012 // Fill the coefficients found with the fits or other
4013 // methods from the Fit functions
4017 if(GetStack(fCountDet) == 2) factor = 1728;
4020 // Pointer to the branch
4021 for (Int_t k = 0; k < factor; k++) {
4022 fCurrentCoefDetector[k] = fCurrentCoef[0];
4023 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4026 FillFillLinearFitter();
4031 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4032 Bool_t AliTRDCalibraFit::FillInfosFitExbAlt()
4035 // Fill the coefficients found with the fits or other
4036 // methods from the Fit functions
4040 if(GetStack(fCountDet) == 2) factor = 1728;
4043 // Pointer to the branch
4044 for (Int_t k = 0; k < factor; k++) {
4045 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4053 //________________________________________________________________________________
4054 void AliTRDCalibraFit::FillFillCH(Int_t idect)
4057 // DebugStream and fVectorFit
4060 // End of one detector
4061 if (idect == (fCount-1)) {
4064 for (Int_t k = 0; k < 2304; k++) {
4065 fCurrentCoefDetector[k] = 0.0;
4069 if(fDebugLevel > 1){
4071 if ( !fDebugStreamer ) {
4073 TDirectory *backup = gDirectory;
4074 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
4075 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4078 Int_t detector = fCountDet;
4079 Int_t caligroup = idect;
4080 Short_t rowmin = fCalibraMode->GetRowMin(0);
4081 Short_t rowmax = fCalibraMode->GetRowMax(0);
4082 Short_t colmin = fCalibraMode->GetColMin(0);
4083 Short_t colmax = fCalibraMode->GetColMax(0);
4084 Float_t gf = fCurrentCoef[0];
4085 Float_t gfs = fCurrentCoef[1];
4086 Float_t gfE = fCurrentCoefE;
4088 (*fDebugStreamer) << "FillFillCH" <<
4089 "detector=" << detector <<
4090 "caligroup=" << caligroup <<
4091 "rowmin=" << rowmin <<
4092 "rowmax=" << rowmax <<
4093 "colmin=" << colmin <<
4094 "colmax=" << colmax <<
4102 //________________________________________________________________________________
4103 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
4106 // DebugStream and fVectorFit and fVectorFit2
4109 // End of one detector
4110 if (idect == (fCount-1)) {
4114 for (Int_t k = 0; k < 2304; k++) {
4115 fCurrentCoefDetector[k] = 0.0;
4116 fCurrentCoefDetector2[k] = 0.0;
4120 if(fDebugLevel > 1){
4122 if ( !fDebugStreamer ) {
4124 TDirectory *backup = gDirectory;
4125 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
4126 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4130 Int_t detector = fCountDet;
4131 Int_t caligroup = idect;
4132 Short_t rowmin = fCalibraMode->GetRowMin(1);
4133 Short_t rowmax = fCalibraMode->GetRowMax(1);
4134 Short_t colmin = fCalibraMode->GetColMin(1);
4135 Short_t colmax = fCalibraMode->GetColMax(1);
4136 Float_t vf = fCurrentCoef[0];
4137 Float_t vs = fCurrentCoef[1];
4138 Float_t vfE = fCurrentCoefE;
4139 Float_t t0f = fCurrentCoef2[0];
4140 Float_t t0s = fCurrentCoef2[1];
4141 Float_t t0E = fCurrentCoefE2;
4145 (* fDebugStreamer) << "FillFillPH"<<
4146 "detector="<<detector<<
4147 "nentries="<<nentries<<
4148 "caligroup="<<caligroup<<
4163 //________________________________________________________________________________
4164 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
4167 // DebugStream and fVectorFit
4170 // End of one detector
4171 if (idect == (fCount-1)) {
4174 for (Int_t k = 0; k < 2304; k++) {
4175 fCurrentCoefDetector[k] = 0.0;
4180 if(fDebugLevel > 1){
4182 if ( !fDebugStreamer ) {
4184 TDirectory *backup = gDirectory;
4185 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4186 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4189 Int_t detector = fCountDet;
4190 Int_t layer = GetLayer(fCountDet);
4191 Int_t caligroup = idect;
4192 Short_t rowmin = fCalibraMode->GetRowMin(2);
4193 Short_t rowmax = fCalibraMode->GetRowMax(2);
4194 Short_t colmin = fCalibraMode->GetColMin(2);
4195 Short_t colmax = fCalibraMode->GetColMax(2);
4196 Float_t widf = fCurrentCoef[0];
4197 Float_t wids = fCurrentCoef[1];
4198 Float_t widfE = fCurrentCoefE;
4200 (* fDebugStreamer) << "FillFillPRF"<<
4201 "detector="<<detector<<
4203 "caligroup="<<caligroup<<
4215 //________________________________________________________________________________
4216 void AliTRDCalibraFit::FillFillLinearFitter()
4219 // DebugStream and fVectorFit
4222 // End of one detector
4228 for (Int_t k = 0; k < 2304; k++) {
4229 fCurrentCoefDetector[k] = 0.0;
4230 fCurrentCoefDetector2[k] = 0.0;
4234 if(fDebugLevel > 1){
4236 if ( !fDebugStreamer ) {
4238 TDirectory *backup = gDirectory;
4239 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
4240 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4243 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4244 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4245 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4246 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4247 Float_t tiltangle = padplane->GetTiltingAngle();
4248 Int_t detector = fCountDet;
4249 Int_t stack = GetStack(fCountDet);
4250 Int_t layer = GetLayer(fCountDet);
4251 Float_t vf = fCurrentCoef[0];
4252 Float_t vs = fCurrentCoef[1];
4253 Float_t vfE = fCurrentCoefE;
4254 Float_t lorentzangler = fCurrentCoef2[0];
4255 Float_t elorentzangler = fCurrentCoefE2;
4256 Float_t lorentzangles = fCurrentCoef2[1];
4258 (* fDebugStreamer) << "FillFillLinearFitter"<<
4259 "detector="<<detector<<
4264 "tiltangle="<<tiltangle<<
4268 "lorentzangler="<<lorentzangler<<
4269 "Elorentzangler="<<elorentzangler<<
4270 "lorentzangles="<<lorentzangles<<
4275 //________________________________________________________________________________
4276 void AliTRDCalibraFit::FillFillExbAlt()
4279 // DebugStream and fVectorFit
4282 // End of one detector
4287 for (Int_t k = 0; k < 2304; k++) {
4288 fCurrentCoefDetector2[k] = 0.0;
4292 if(fDebugLevel > 1){
4294 if ( !fDebugStreamer ) {
4296 TDirectory *backup = gDirectory;
4297 fDebugStreamer = new TTreeSRedirector("TRDDebugFitExbAlt.root");
4298 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4301 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4302 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4303 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4304 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4305 Float_t tiltangle = padplane->GetTiltingAngle();
4306 Int_t detector = fCountDet;
4307 Int_t stack = GetStack(fCountDet);
4308 Int_t layer = GetLayer(fCountDet);
4309 Float_t vf = fCurrentCoef2[0];
4310 Float_t vfE = fCurrentCoefE2;
4312 (* fDebugStreamer) << "FillFillLinearFitter"<<
4313 "detector="<<detector<<
4318 "tiltangle="<<tiltangle<<
4326 //____________Calcul Coef Mean_________________________________________________
4328 //_____________________________________________________________________________
4329 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
4332 // For the detector Dect calcul the mean time 0
4333 // for the calibration group idect from the choosen database
4336 fCurrentCoef2[1] = 0.0;
4337 if(fDebugLevel != 1){
4338 if(((fCalibraMode->GetNz(1) > 0) ||
4339 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4341 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4342 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4343 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4347 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4353 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4357 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
4358 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
4359 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4362 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
4370 //_____________________________________________________________________________
4371 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
4374 // For the detector Dect calcul the mean gain factor
4375 // for the calibration group idect from the choosen database
4378 fCurrentCoef[1] = 0.0;
4379 if(fDebugLevel != 1){
4380 if (((fCalibraMode->GetNz(0) > 0) ||
4381 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
4382 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4383 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4384 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4385 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4388 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4392 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4393 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
4398 //_____________________________________________________________________________
4399 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
4402 // For the detector Dect calcul the mean sigma of pad response
4403 // function for the calibration group idect from the choosen database
4406 fCurrentCoef[1] = 0.0;
4407 if(fDebugLevel != 1){
4408 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4409 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4410 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
4413 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
4417 //_____________________________________________________________________________
4418 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
4421 // For the detector dect calcul the mean drift velocity for the
4422 // calibration group idect from the choosen database
4425 fCurrentCoef[1] = 0.0;
4426 if(fDebugLevel != 1){
4427 if (((fCalibraMode->GetNz(1) > 0) ||
4428 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4430 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4431 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4432 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4436 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4441 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4446 //_____________________________________________________________________________
4447 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4450 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4453 fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet);
4454 fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet);
4458 //_____________________________________________________________________________
4459 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4462 // Default width of the PRF if there is no database as reference
4467 //case 0: return 0.515;
4468 //case 1: return 0.502;
4469 //case 2: return 0.491;
4470 //case 3: return 0.481;
4471 //case 4: return 0.471;
4472 //case 5: return 0.463;
4473 //default: return 0.0;
4476 case 0: return 0.538429;
4477 case 1: return 0.524302;
4478 case 2: return 0.511591;
4479 case 3: return 0.500140;
4480 case 4: return 0.489821;
4481 case 5: return 0.480524;
4482 default: return 0.0;
4485 //________________________________________________________________________________
4486 void AliTRDCalibraFit::SetCalROC(Int_t i)
4489 // Set the calib object for fCountDet
4492 Float_t value = 0.0;
4494 //Get the CalDet object
4496 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4498 AliInfo("Could not get calibDB");
4505 fCalROC->~AliTRDCalROC();
4506 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4507 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4511 fCalROC->~AliTRDCalROC();
4512 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4513 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4515 fCalROC2->~AliTRDCalROC();
4516 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4517 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4521 fCalROC->~AliTRDCalROC();
4522 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4523 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4532 if(fCalROC) delete fCalROC;
4533 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4534 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4535 fCalROC->SetValue(k,1.0);
4539 if(fCalROC) delete fCalROC;
4540 if(fCalROC2) delete fCalROC2;
4541 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4542 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4543 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4544 fCalROC->SetValue(k,1.0);
4545 fCalROC2->SetValue(k,0.0);
4549 if(fCalROC) delete fCalROC;
4550 value = GetPRFDefault(GetLayer(fCountDet));
4551 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4552 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4553 fCalROC->SetValue(k,value);
4561 //____________Fit Methods______________________________________________________
4563 //_____________________________________________________________________________
4564 void AliTRDCalibraFit::FitPente(TH1* projPH)
4567 // Slope methode for the drift velocity
4571 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4578 fCurrentCoefE = 0.0;
4579 fCurrentCoefE2 = 0.0;
4580 fCurrentCoef[0] = 0.0;
4581 fCurrentCoef2[0] = 0.0;
4582 TLine *line = new TLine();
4585 TAxis *xpph = projPH->GetXaxis();
4586 Int_t nbins = xpph->GetNbins();
4587 Double_t lowedge = xpph->GetBinLowEdge(1);
4588 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4589 Double_t widbins = (upedge - lowedge) / nbins;
4590 Double_t limit = upedge + 0.5 * widbins;
4593 // Beginning of the signal
4594 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4595 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4596 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4598 binmax = (Int_t) pentea->GetMaximumBin();
4601 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4603 if (binmax >= nbins) {
4606 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4608 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4609 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4610 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4611 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4612 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4613 if (TMath::Abs(l3P2am) > 0.00000001) {
4614 fPhd[0] = -(l3P1am / (2 * l3P2am));
4617 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4618 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4621 // Amplification region
4624 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4625 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4632 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4634 if (binmax >= nbins) {
4637 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4639 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4640 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4641 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4642 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4643 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4644 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4645 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4647 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4648 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4651 fCurrentCoefE2 = fCurrentCoefE;
4654 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4655 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4656 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4659 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4662 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4664 if (binmin >= nbins) {
4667 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4672 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4673 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4674 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4675 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4676 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4677 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4678 if (TMath::Abs(l3P2dr) > 0.00000001) {
4679 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4681 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4682 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4684 Float_t fPhdt0 = 0.0;
4685 Float_t t0Shift = 0.0;
4688 t0Shift = fT0Shift1;
4692 t0Shift = fT0Shift0;
4695 if ((fPhd[2] > fPhd[0]) &&
4696 (fPhd[2] > fPhd[1]) &&
4697 (fPhd[1] > fPhd[0]) &&
4699 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4700 fNumberFitSuccess++;
4702 if (fPhdt0 >= 0.0) {
4703 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4704 if (fCurrentCoef2[0] < -3.0) {
4705 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4709 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4714 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4715 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4718 if (fDebugLevel == 1) {
4719 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4722 line->SetLineColor(2);
4723 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4724 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4725 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4726 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4727 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4728 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4729 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4730 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4733 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4742 //_____________________________________________________________________________
4743 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4746 // Slope methode but with polynomes de Lagrange
4750 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4753 //Double_t *x = new Double_t[5];
4754 //Double_t *y = new Double_t[5];
4771 fCurrentCoefE = 0.0;
4772 fCurrentCoefE2 = 1.0;
4773 fCurrentCoef[0] = 0.0;
4774 fCurrentCoef2[0] = 0.0;
4775 TLine *line = new TLine();
4776 TF1 * polynome = 0x0;
4777 TF1 * polynomea = 0x0;
4778 TF1 * polynomeb = 0x0;
4786 TAxis *xpph = projPH->GetXaxis();
4787 Int_t nbins = xpph->GetNbins();
4788 Double_t lowedge = xpph->GetBinLowEdge(1);
4789 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4790 Double_t widbins = (upedge - lowedge) / nbins;
4791 Double_t limit = upedge + 0.5 * widbins;
4796 ///////////////////////////////
4797 // Beginning of the signal
4798 //////////////////////////////
4799 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4800 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4801 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4804 binmax = (Int_t) pentea->GetMaximumBin();
4806 Double_t minnn = 0.0;
4807 Double_t maxxx = 0.0;
4809 Int_t kase = nbins-binmax;
4817 minnn = pentea->GetBinCenter(binmax-2);
4818 maxxx = pentea->GetBinCenter(binmax);
4819 x[0] = pentea->GetBinCenter(binmax-2);
4820 x[1] = pentea->GetBinCenter(binmax-1);
4821 x[2] = pentea->GetBinCenter(binmax);
4822 y[0] = pentea->GetBinContent(binmax-2);
4823 y[1] = pentea->GetBinContent(binmax-1);
4824 y[2] = pentea->GetBinContent(binmax);
4825 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4826 //AliInfo("At the limit for beginning!");
4829 minnn = pentea->GetBinCenter(binmax-2);
4830 maxxx = pentea->GetBinCenter(binmax+1);
4831 x[0] = pentea->GetBinCenter(binmax-2);
4832 x[1] = pentea->GetBinCenter(binmax-1);
4833 x[2] = pentea->GetBinCenter(binmax);
4834 x[3] = pentea->GetBinCenter(binmax+1);
4835 y[0] = pentea->GetBinContent(binmax-2);
4836 y[1] = pentea->GetBinContent(binmax-1);
4837 y[2] = pentea->GetBinContent(binmax);
4838 y[3] = pentea->GetBinContent(binmax+1);
4839 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4847 minnn = pentea->GetBinCenter(binmax);
4848 maxxx = pentea->GetBinCenter(binmax+2);
4849 x[0] = pentea->GetBinCenter(binmax);
4850 x[1] = pentea->GetBinCenter(binmax+1);
4851 x[2] = pentea->GetBinCenter(binmax+2);
4852 y[0] = pentea->GetBinContent(binmax);
4853 y[1] = pentea->GetBinContent(binmax+1);
4854 y[2] = pentea->GetBinContent(binmax+2);
4855 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4858 minnn = pentea->GetBinCenter(binmax-1);
4859 maxxx = pentea->GetBinCenter(binmax+2);
4860 x[0] = pentea->GetBinCenter(binmax-1);
4861 x[1] = pentea->GetBinCenter(binmax);
4862 x[2] = pentea->GetBinCenter(binmax+1);
4863 x[3] = pentea->GetBinCenter(binmax+2);
4864 y[0] = pentea->GetBinContent(binmax-1);
4865 y[1] = pentea->GetBinContent(binmax);
4866 y[2] = pentea->GetBinContent(binmax+1);
4867 y[3] = pentea->GetBinContent(binmax+2);
4868 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4871 minnn = pentea->GetBinCenter(binmax-2);
4872 maxxx = pentea->GetBinCenter(binmax+2);
4873 x[0] = pentea->GetBinCenter(binmax-2);
4874 x[1] = pentea->GetBinCenter(binmax-1);
4875 x[2] = pentea->GetBinCenter(binmax);
4876 x[3] = pentea->GetBinCenter(binmax+1);
4877 x[4] = pentea->GetBinCenter(binmax+2);
4878 y[0] = pentea->GetBinContent(binmax-2);
4879 y[1] = pentea->GetBinContent(binmax-1);
4880 y[2] = pentea->GetBinContent(binmax);
4881 y[3] = pentea->GetBinContent(binmax+1);
4882 y[4] = pentea->GetBinContent(binmax+2);
4883 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4891 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4892 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4894 Double_t step = (maxxx-minnn)/10000;
4896 Double_t maxvalue = 0.0;
4897 Double_t placemaximum = minnn;
4898 for(Int_t o = 0; o < 10000; o++){
4899 if(o == 0) maxvalue = polynomeb->Eval(l);
4900 if(maxvalue < (polynomeb->Eval(l))){
4901 maxvalue = polynomeb->Eval(l);
4906 fPhd[0] = placemaximum;
4909 /////////////////////////////
4910 // Amplification region
4911 /////////////////////////////
4914 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4915 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4921 Double_t minn = 0.0;
4922 Double_t maxx = 0.0;
4934 Int_t kase1 = nbins - binmax;
4936 //Determination of minn and maxx
4937 //case binmax = nbins
4942 minn = projPH->GetBinCenter(binmax-2);
4943 maxx = projPH->GetBinCenter(binmax);
4944 x[0] = projPH->GetBinCenter(binmax-2);
4945 x[1] = projPH->GetBinCenter(binmax-1);
4946 x[2] = projPH->GetBinCenter(binmax);
4947 y[0] = projPH->GetBinContent(binmax-2);
4948 y[1] = projPH->GetBinContent(binmax-1);
4949 y[2] = projPH->GetBinContent(binmax);
4950 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4951 //AliInfo("At the limit for the drift!");
4954 minn = projPH->GetBinCenter(binmax-2);
4955 maxx = projPH->GetBinCenter(binmax+1);
4956 x[0] = projPH->GetBinCenter(binmax-2);
4957 x[1] = projPH->GetBinCenter(binmax-1);
4958 x[2] = projPH->GetBinCenter(binmax);
4959 x[3] = projPH->GetBinCenter(binmax+1);
4960 y[0] = projPH->GetBinContent(binmax-2);
4961 y[1] = projPH->GetBinContent(binmax-1);
4962 y[2] = projPH->GetBinContent(binmax);
4963 y[3] = projPH->GetBinContent(binmax+1);
4964 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4973 minn = projPH->GetBinCenter(binmax);
4974 maxx = projPH->GetBinCenter(binmax+2);
4975 x[0] = projPH->GetBinCenter(binmax);
4976 x[1] = projPH->GetBinCenter(binmax+1);
4977 x[2] = projPH->GetBinCenter(binmax+2);
4978 y[0] = projPH->GetBinContent(binmax);
4979 y[1] = projPH->GetBinContent(binmax+1);
4980 y[2] = projPH->GetBinContent(binmax+2);
4981 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4984 minn = projPH->GetBinCenter(binmax-1);
4985 maxx = projPH->GetBinCenter(binmax+2);
4986 x[0] = projPH->GetBinCenter(binmax-1);
4987 x[1] = projPH->GetBinCenter(binmax);
4988 x[2] = projPH->GetBinCenter(binmax+1);
4989 x[3] = projPH->GetBinCenter(binmax+2);
4990 y[0] = projPH->GetBinContent(binmax-1);
4991 y[1] = projPH->GetBinContent(binmax);
4992 y[2] = projPH->GetBinContent(binmax+1);
4993 y[3] = projPH->GetBinContent(binmax+2);
4994 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4997 minn = projPH->GetBinCenter(binmax-2);
4998 maxx = projPH->GetBinCenter(binmax+2);
4999 x[0] = projPH->GetBinCenter(binmax-2);
5000 x[1] = projPH->GetBinCenter(binmax-1);
5001 x[2] = projPH->GetBinCenter(binmax);
5002 x[3] = projPH->GetBinCenter(binmax+1);
5003 x[4] = projPH->GetBinCenter(binmax+2);
5004 y[0] = projPH->GetBinContent(binmax-2);
5005 y[1] = projPH->GetBinContent(binmax-1);
5006 y[2] = projPH->GetBinContent(binmax);
5007 y[3] = projPH->GetBinContent(binmax+1);
5008 y[4] = projPH->GetBinContent(binmax+2);
5009 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5016 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
5017 polynomea->SetParameters(c0,c1,c2,c3,c4);
5019 Double_t step = (maxx-minn)/1000;
5021 Double_t maxvalue = 0.0;
5022 Double_t placemaximum = minn;
5023 for(Int_t o = 0; o < 1000; o++){
5024 if(o == 0) maxvalue = polynomea->Eval(l);
5025 if(maxvalue < (polynomea->Eval(l))){
5026 maxvalue = polynomea->Eval(l);
5031 fPhd[1] = placemaximum;
5037 Bool_t putd = kTRUE;
5038 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
5039 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5040 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5043 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5049 //AliInfo("Put the binmax from 1 to 2 to enable the fit");
5053 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
5054 //AliInfo("Too many fluctuations at the end!");
5057 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
5058 //AliInfo("Too many fluctuations at the end!");
5061 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
5062 //AliInfo("No entries for the next bin!");
5063 pente->SetBinContent(binmin,0);
5064 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5080 Bool_t case1 = kFALSE;
5081 Bool_t case2 = kFALSE;
5082 Bool_t case4 = kFALSE;
5084 //Determination of min and max
5085 //case binmin <= nbins-3
5087 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5088 min = pente->GetBinCenter(binmin-2);
5089 max = pente->GetBinCenter(binmin+2);
5090 x[0] = pente->GetBinCenter(binmin-2);
5091 x[1] = pente->GetBinCenter(binmin-1);
5092 x[2] = pente->GetBinCenter(binmin);
5093 x[3] = pente->GetBinCenter(binmin+1);
5094 x[4] = pente->GetBinCenter(binmin+2);
5095 y[0] = pente->GetBinContent(binmin-2);
5096 y[1] = pente->GetBinContent(binmin-1);
5097 y[2] = pente->GetBinContent(binmin);
5098 y[3] = pente->GetBinContent(binmin+1);
5099 y[4] = pente->GetBinContent(binmin+2);
5100 //Calcul the polynome de Lagrange
5101 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5103 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5104 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5105 //AliInfo("polynome 4 false 1");
5108 if(((binmin+3) <= (nbins-1)) &&
5109 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
5110 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
5111 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
5112 //AliInfo("polynome 4 false 2");
5116 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5117 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
5118 //AliInfo("polynome 4 case 1");
5121 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
5122 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5123 //AliInfo("polynome 4 case 4");
5128 //case binmin = nbins-2
5130 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5132 min = pente->GetBinCenter(binmin-2);
5133 max = pente->GetBinCenter(binmin+1);
5134 x[0] = pente->GetBinCenter(binmin-2);
5135 x[1] = pente->GetBinCenter(binmin-1);
5136 x[2] = pente->GetBinCenter(binmin);
5137 x[3] = pente->GetBinCenter(binmin+1);
5138 y[0] = pente->GetBinContent(binmin-2);
5139 y[1] = pente->GetBinContent(binmin-1);
5140 y[2] = pente->GetBinContent(binmin);
5141 y[3] = pente->GetBinContent(binmin+1);
5142 //Calcul the polynome de Lagrange
5143 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5144 //richtung +: nothing
5146 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5147 //AliInfo("polynome 3- case 2");
5152 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5154 min = pente->GetBinCenter(binmin-1);
5155 max = pente->GetBinCenter(binmin+2);
5156 x[0] = pente->GetBinCenter(binmin-1);
5157 x[1] = pente->GetBinCenter(binmin);
5158 x[2] = pente->GetBinCenter(binmin+1);
5159 x[3] = pente->GetBinCenter(binmin+2);
5160 y[0] = pente->GetBinContent(binmin-1);
5161 y[1] = pente->GetBinContent(binmin);
5162 y[2] = pente->GetBinContent(binmin+1);
5163 y[3] = pente->GetBinContent(binmin+2);
5164 //Calcul the polynome de Lagrange
5165 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5167 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5168 //AliInfo("polynome 3+ case 2");
5173 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
5174 min = pente->GetBinCenter(binmin);
5175 max = pente->GetBinCenter(binmin+2);
5176 x[0] = pente->GetBinCenter(binmin);
5177 x[1] = pente->GetBinCenter(binmin+1);
5178 x[2] = pente->GetBinCenter(binmin+2);
5179 y[0] = pente->GetBinContent(binmin);
5180 y[1] = pente->GetBinContent(binmin+1);
5181 y[2] = pente->GetBinContent(binmin+2);
5182 //Calcul the polynome de Lagrange
5183 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5185 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5186 //AliInfo("polynome 2+ false");
5191 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5193 min = pente->GetBinCenter(binmin-1);
5194 max = pente->GetBinCenter(binmin+1);
5195 x[0] = pente->GetBinCenter(binmin-1);
5196 x[1] = pente->GetBinCenter(binmin);
5197 x[2] = pente->GetBinCenter(binmin+1);
5198 y[0] = pente->GetBinContent(binmin-1);
5199 y[1] = pente->GetBinContent(binmin);
5200 y[2] = pente->GetBinContent(binmin+1);
5201 //Calcul the polynome de Lagrange
5202 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5203 //richtung +: nothing
5204 //richtung -: nothing
5206 //case binmin = nbins-1
5208 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5209 min = pente->GetBinCenter(binmin-2);
5210 max = pente->GetBinCenter(binmin);
5211 x[0] = pente->GetBinCenter(binmin-2);
5212 x[1] = pente->GetBinCenter(binmin-1);
5213 x[2] = pente->GetBinCenter(binmin);
5214 y[0] = pente->GetBinContent(binmin-2);
5215 y[1] = pente->GetBinContent(binmin-1);
5216 y[2] = pente->GetBinContent(binmin);
5217 //Calcul the polynome de Lagrange
5218 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5219 //AliInfo("At the limit for the drift!");
5220 //fluctuation too big!
5221 //richtung +: nothing
5223 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5224 //AliInfo("polynome 2- false ");
5228 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
5230 //AliInfo("At the limit for the drift and not usable!");
5234 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5236 //AliInfo("For the drift...problem!");
5238 //pass but should not happen
5239 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
5241 //AliInfo("For the drift...problem!");
5245 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5246 polynome->SetParameters(c0,c1,c2,c3,c4);
5247 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5248 Double_t step = (max-min)/1000;
5250 Double_t minvalue = 0.0;
5251 Double_t placeminimum = min;
5252 for(Int_t o = 0; o < 1000; o++){
5253 if(o == 0) minvalue = polynome->Eval(l);
5254 if(minvalue > (polynome->Eval(l))){
5255 minvalue = polynome->Eval(l);
5260 fPhd[2] = placeminimum;
5262 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
5263 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
5264 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
5266 Float_t fPhdt0 = 0.0;
5267 Float_t t0Shift = 0.0;
5270 t0Shift = fT0Shift1;
5274 t0Shift = fT0Shift0;
5277 if ((fPhd[2] > fPhd[0]) &&
5278 (fPhd[2] > fPhd[1]) &&
5279 (fPhd[1] > fPhd[0]) &&
5282 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5283 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5284 else fNumberFitSuccess++;
5285 if (fPhdt0 >= 0.0) {
5286 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5287 //printf("Value of timeoffset %f\n",fCurrentCoef2[0]);
5288 if (fCurrentCoef2[0] < -3.0) {
5289 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5293 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5297 ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
5298 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5300 if((fPhd[1] > fPhd[0]) &&
5302 if (fPhdt0 >= 0.0) {
5303 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5304 if (fCurrentCoef2[0] < -3.0) {
5305 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5307 else fNumberFitSuccess++;
5310 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5314 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5315 //printf("Fit failed!\n");
5319 if (fDebugLevel == 1) {
5320 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5323 if(polynomea) polynomea->Draw("same");
5324 line->SetLineColor(2);
5325 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5326 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5327 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5328 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5329 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5330 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5331 AliInfo(Form("Vdrift (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
5332 AliInfo(Form("Timeoffset (with only the drift region(default)): %f",(Float_t) fCurrentCoef2[0]));
5333 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5336 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5343 if(polynome) delete polynome;
5344 if(polynomea) delete polynomea;
5345 if(polynomeb) delete polynomeb;
5346 //if(x) delete [] x;
5347 //if(y) delete [] y;
5348 if(line) delete line;
5353 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5354 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5355 //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]);
5356 projPH->SetDirectory(0);
5360 //_____________________________________________________________________________
5361 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5364 // Fit methode for the drift velocity
5368 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5371 TAxis *xpph = projPH->GetXaxis();
5372 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5374 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5375 fPH->SetParameter(0,0.469); // Scaling
5376 fPH->SetParameter(1,0.18); // Start
5377 fPH->SetParameter(2,0.0857325); // AR
5378 fPH->SetParameter(3,1.89); // DR
5379 fPH->SetParameter(4,0.08); // QA/QD
5380 fPH->SetParameter(5,0.0); // Baseline
5382 TLine *line = new TLine();
5384 fCurrentCoef[0] = 0.0;
5385 fCurrentCoef2[0] = 0.0;
5386 fCurrentCoefE = 0.0;
5387 fCurrentCoefE2 = 0.0;
5389 if (idect%fFitPHPeriode == 0) {
5391 AliInfo(Form("The detector %d will be fitted",idect));
5392 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5393 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5394 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5395 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5396 fPH->SetParameter(4,0.225); // QA/QD
5397 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5399 if (fDebugLevel != 1) {
5400 projPH->Fit(fPH,"0M","",0.0,upedge);
5403 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5405 projPH->Fit(fPH,"M+","",0.0,upedge);
5407 line->SetLineColor(4);
5408 line->DrawLine(fPH->GetParameter(1)
5410 ,fPH->GetParameter(1)
5411 ,projPH->GetMaximum());
5412 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5414 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5415 ,projPH->GetMaximum());
5416 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5418 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5419 ,projPH->GetMaximum());
5422 if (fPH->GetParameter(3) != 0) {
5423 fNumberFitSuccess++;
5424 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
5425 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
5426 fCurrentCoef2[0] = fPH->GetParameter(1);
5427 fCurrentCoefE2 = fPH->GetParError(1);
5430 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5431 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5437 // Put the default value
5438 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5439 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5442 if (fDebugLevel != 1) {
5447 //_____________________________________________________________________________
5448 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5451 // Fit methode for the sigma of the pad response function
5456 fCurrentCoef[0] = 0.0;
5457 fCurrentCoefE = 0.0;
5459 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5461 if(TMath::Abs(ret+4) <= 0.000000001){
5462 fCurrentCoef[0] = -fCurrentCoef[1];
5466 fNumberFitSuccess++;
5467 fCurrentCoef[0] = param[2];
5468 fCurrentCoefE = ret;
5472 //_____________________________________________________________________________
5473 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)
5476 // Fit methode for the sigma of the pad response function
5479 //We should have at least 3 points
5480 if(nBins <=3) return -4.0;
5482 TLinearFitter fitter(3,"pol2");
5483 fitter.StoreData(kFALSE);
5484 fitter.ClearPoints();
5486 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5487 Float_t entries = 0;
5488 Int_t nbbinwithentries = 0;
5489 for (Int_t i=0; i<nBins; i++){
5491 if(arraye[i] > 15) nbbinwithentries++;
5492 //printf("entries for i %d: %f\n",i,arraye[i]);
5494 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5495 //printf("entries %f\n",entries);
5496 //printf("nbbinwithentries %d\n",nbbinwithentries);
5499 Float_t errorm = 0.0;
5500 Float_t errorn = 0.0;
5501 Float_t error = 0.0;
5504 for (Int_t ibin=0;ibin<nBins; ibin++){
5505 Float_t entriesI = arraye[ibin];
5506 Float_t valueI = arraym[ibin];
5507 Double_t xcenter = 0.0;
5509 if ((entriesI>15) && (valueI>0.0)){
5510 xcenter = xMin+(ibin+0.5)*binWidth;
5515 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5516 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5517 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5520 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5521 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5522 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5524 if(TMath::Abs(error) < 0.000000001) continue;
5525 val = TMath::Log(Float_t(valueI));
5526 fitter.AddPoint(&xcenter,val,error);
5530 if(fDebugLevel > 1){
5532 if ( !fDebugStreamer ) {
5534 TDirectory *backup = gDirectory;
5535 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5536 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5539 Int_t detector = fCountDet;
5540 Int_t layer = GetLayer(fCountDet);
5543 (* fDebugStreamer) << "FitGausMIFill"<<
5544 "detector="<<detector<<
5548 "entriesI="<<entriesI<<
5551 "xcenter="<<xcenter<<
5561 if(npoints <=3) return -4.0;
5566 fitter.GetParameters(par);
5567 chi2 = fitter.GetChisquare()/Float_t(npoints);
5570 if (!param) param = new TVectorD(3);
5571 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5572 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5573 Double_t deltax = (fitter.GetParError(2))/x;
5574 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5577 (*param)[1] = par[1]/(-2.*par[2]);
5578 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5579 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5580 if ( lnparam0>307 ) return -4;
5581 (*param)[0] = TMath::Exp(lnparam0);
5583 if(fDebugLevel > 1){
5585 if ( !fDebugStreamer ) {
5587 TDirectory *backup = gDirectory;
5588 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5589 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5592 Int_t detector = fCountDet;
5593 Int_t layer = GetLayer(fCountDet);
5596 (* fDebugStreamer) << "FitGausMIFit"<<
5597 "detector="<<detector<<
5600 "errorsigma="<<chi2<<
5601 "mean="<<(*param)[1]<<
5602 "sigma="<<(*param)[2]<<
5603 "constant="<<(*param)[0]<<
5608 if((chi2/(*param)[2]) > 0.1){
5610 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5615 if(fDebugLevel == 1){
5616 TString name("PRF");
5617 name += (Int_t)xMin;
5618 name += (Int_t)xMax;
5619 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5622 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5623 for(Int_t k = 0; k < nBins; k++){
5624 histo->SetBinContent(k+1,arraym[k]);
5625 histo->SetBinError(k+1,arrayme[k]);
5628 name += "functionf";
5629 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5630 f1->SetParameter(0, (*param)[0]);
5631 f1->SetParameter(1, (*param)[1]);
5632 f1->SetParameter(2, (*param)[2]);
5640 //_____________________________________________________________________________
5641 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5644 // Fit methode for the sigma of the pad response function
5647 fCurrentCoef[0] = 0.0;
5648 fCurrentCoefE = 0.0;
5650 if (fDebugLevel != 1) {
5651 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5654 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5656 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5659 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5660 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5661 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5663 fNumberFitSuccess++;
5666 //_____________________________________________________________________________
5667 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5670 // Fit methode for the sigma of the pad response function
5672 fCurrentCoef[0] = 0.0;
5673 fCurrentCoefE = 0.0;
5674 if (fDebugLevel == 1) {
5675 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5679 fCurrentCoef[0] = projPRF->GetRMS();
5680 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5682 fNumberFitSuccess++;
5685 //_____________________________________________________________________________
5686 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5689 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5692 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5695 Int_t nbins = (Int_t)(nybins/(2*nbg));
5696 Float_t lowedge = -3.0*nbg;
5697 Float_t upedge = lowedge + 3.0;
5700 Double_t xvalues = -0.2*nbg+0.1;
5702 Int_t total = 2*nbg;
5705 for(Int_t k = 0; k < total; k++){
5706 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5708 y = fCurrentCoef[0]*fCurrentCoef[0];
5709 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5712 if(fDebugLevel > 1){
5714 if ( !fDebugStreamer ) {
5716 TDirectory *backup = gDirectory;
5717 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5718 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5721 Int_t detector = fCountDet;
5722 Int_t layer = GetLayer(fCountDet);
5723 Int_t nbtotal = total;
5725 Float_t low = lowedge;
5726 Float_t up = upedge;
5727 Float_t tnp = xvalues;
5728 Float_t wid = fCurrentCoef[0];
5729 Float_t widfE = fCurrentCoefE;
5731 (* fDebugStreamer) << "FitTnpRange0"<<
5732 "detector="<<detector<<
5734 "nbtotal="<<nbtotal<<
5752 fCurrentCoefE = 0.0;
5753 fCurrentCoef[0] = 0.0;
5755 //printf("npoints\n",npoints);
5758 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5763 linearfitter.Eval();
5764 linearfitter.GetParameters(pars0);
5766 //Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5767 //Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5768 Double_t min0 = 0.0;
5770 //Double_t ermin0 = 0.0;
5771 //Double_t prfe0 = 0.0;
5772 Double_t prf0 = 0.0;
5773 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5774 min0 = -pars0[1]/(2*pars0[2]);
5776 //ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5777 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5780 prfe0 = linearfitter->GetParError(0)*pointError0
5781 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5782 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5783 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5784 fCurrentCoefE = (Float_t) prfe0;
5786 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5789 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5793 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5796 if(fDebugLevel > 1){
5798 if ( !fDebugStreamer ) {
5800 TDirectory *backup = gDirectory;
5801 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5802 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5805 Int_t detector = fCountDet;
5806 Int_t layer = GetLayer(fCountDet);
5807 Int_t nbtotal = total;
5808 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5809 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5811 (* fDebugStreamer) << "FitTnpRange1"<<
5812 "detector="<<detector<<
5814 "nbtotal="<<nbtotal<<
5818 "npoints="<<npoints<<
5821 "sigmaprf="<<fCurrentCoef[0]<<
5822 "sigprf="<<fCurrentCoef[1]<<
5829 //_____________________________________________________________________________
5830 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5833 // Only mean methode for the gain factor
5836 fCurrentCoef[0] = mean;
5837 fCurrentCoefE = 0.0;
5838 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5839 if (fDebugLevel == 1) {
5840 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5844 CalculChargeCoefMean(kTRUE);
5845 fNumberFitSuccess++;
5847 //_____________________________________________________________________________
5848 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5851 // mean w methode for the gain factor
5855 Int_t nybins = projch->GetNbinsX();
5857 //The weight function
5858 Double_t a = 0.00228515;
5859 Double_t b = -0.00231487;
5860 Double_t c = 0.00044298;
5861 Double_t d = -0.00379239;
5862 Double_t e = 0.00338349;
5872 //A arbitrary error for the moment
5873 fCurrentCoefE = 0.0;
5874 fCurrentCoef[0] = 0.0;
5877 Double_t sumw = 0.0;
5879 Float_t sumAll = (Float_t) nentries;
5880 Int_t sumCurrent = 0;
5881 for(Int_t k = 0; k <nybins; k++){
5882 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5883 if (fraction<fOutliersFitChargeLow) {
5884 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5885 //printf("Take only after bin %d\n",k);
5888 if (fraction>fOutliersFitChargeHigh) {
5889 //printf("Break by the bin %d\n",k);
5892 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5893 e*fraction*fraction*fraction*fraction;
5894 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5895 sum += weight*projch->GetBinContent(k+1);
5896 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5897 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5899 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5901 if (fDebugLevel == 1) {
5902 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5905 TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
5908 fNumberFitSuccess++;
5909 CalculChargeCoefMean(kTRUE);
5911 //_____________________________________________________________________________
5912 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5915 // mean w methode for the gain factor
5919 Int_t nybins = projch->GetNbinsX();
5921 //The weight function
5922 Double_t a = 0.00228515;
5923 Double_t b = -0.00231487;
5924 Double_t c = 0.00044298;
5925 Double_t d = -0.00379239;
5926 Double_t e = 0.00338349;
5936 //A arbitrary error for the moment
5937 fCurrentCoefE = 0.0;
5938 fCurrentCoef[0] = 0.0;
5941 Double_t sumw = 0.0;
5943 Int_t sumCurrent = 0;
5944 for(Int_t k = 0; k <nybins; k++){
5945 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5946 if (fraction>0.95) break;
5947 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5948 e*fraction*fraction*fraction*fraction;
5949 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5950 sum += weight*projch->GetBinContent(k+1);
5951 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5952 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5954 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5956 if (fDebugLevel == 1) {
5957 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5960 TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
5963 fNumberFitSuccess++;
5965 //_____________________________________________________________________________
5966 void AliTRDCalibraFit::FitLandau(TH1 *projch, Double_t mean, Double_t nentries)
5969 // Fit methode for the gain factor
5973 //Calcul Range of the fit
5974 Double_t lastvalue = 0.0;
5975 Float_t sumAll = (Float_t) nentries;
5976 Int_t sumCurrent = 0;
5977 //printf("There are %d bins\n",nybins);
5978 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
5979 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5980 if (fraction>fOutliersFitChargeHigh) {
5981 lastvalue = projch->GetBinCenter(k+1);
5982 //printf("Break by %f\n",lastvalue);
5985 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5989 fCurrentCoef[0] = 0.0;
5990 fCurrentCoefE = 0.0;
5992 //Double_t chisqrl = 0.0;
5994 projch->Fit("landau","WWQ+",""
5995 ,(Double_t) mean/fBeginFitCharge
5998 //chisqrl = projch->GetFunction("landau")->GetChisquare();
6000 if (fDebugLevel == 1) {
6001 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
6004 TLine *line = new TLine( projch->GetFunction("landau")->GetParameter(1),0.0,projch->GetFunction("landau")->GetParameter(1),20000.0);
6008 if ((projch->GetFunction("landau")->GetParameter(1) > 0) && (projch->GetFunction("landau")->GetParError(1) < (0.05*projch->GetFunction("landau")->GetParameter(1)))) {
6009 fNumberFitSuccess++;
6010 CalculChargeCoefMean(kTRUE);
6011 fCurrentCoef[0] = projch->GetFunction("landau")->GetParameter(1);
6012 fCurrentCoefE = projch->GetFunction("landau")->GetParError(1);
6015 CalculChargeCoefMean(kFALSE);
6016 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6022 //_____________________________________________________________________________
6023 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean, Double_t nentries)
6026 // Fit methode for the gain factor
6029 //Calcul Range of the fit
6030 Double_t lastvalue = 0.0;
6031 Float_t sumAll = (Float_t) nentries;
6032 Int_t sumCurrent = 0;
6033 //printf("There are %d bins\n",nybins);
6034 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6035 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6036 if (fraction>fOutliersFitChargeHigh) {
6037 lastvalue = projch->GetBinCenter(k+1);
6038 //printf("Break by %f\n",lastvalue);
6041 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6045 fCurrentCoef[0] = 0.0;
6046 fCurrentCoefE = 0.0;
6047 Double_t chisqrl = 0.0;
6048 Double_t chisqrg = 0.0;
6049 Double_t chisqr = 0.0;
6050 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,(Double_t) mean/fBeginFitCharge,lastvalue,5);
6052 projch->Fit("landau","WWQ0",""
6053 ,(Double_t) mean/fBeginFitCharge
6055 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
6056 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
6057 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
6058 chisqrl = projch->GetFunction("landau")->GetChisquare();
6060 projch->Fit("gaus","WWQ0",""
6061 ,(Double_t) mean/fBeginFitCharge
6063 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
6064 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
6065 chisqrg = projch->GetFunction("gaus")->GetChisquare();
6067 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
6068 if (fDebugLevel != 1) {
6069 projch->Fit("fLandauGaus","WWQ0",""
6070 ,(Double_t) mean/fBeginFitCharge
6072 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6075 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
6077 projch->Fit("fLandauGaus","WWQ+",""
6078 ,(Double_t) mean/fBeginFitCharge
6080 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6082 fLandauGaus->Draw("same");
6083 TLine *line = new TLine(projch->GetFunction("fLandauGaus")->GetParameter(1),0.0,projch->GetFunction("fLandauGaus")->GetParameter(1),20000.0);
6087 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
6088 fNumberFitSuccess++;
6089 CalculChargeCoefMean(kTRUE);
6090 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
6091 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
6094 CalculChargeCoefMean(kFALSE);
6095 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6098 if (fDebugLevel != 1) {
6103 //_____________________________________________________________________________
6104 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean, Double_t nentries)
6107 // Fit methode for the gain factor more time consuming
6110 //Calcul Range of the fit
6111 Double_t lastvalue = 0.0;
6112 Float_t sumAll = (Float_t) nentries;
6113 Int_t sumCurrent = 0;
6114 //printf("There are %d bins\n",nybins);
6115 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6116 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6117 if (fraction>fOutliersFitChargeHigh) {
6118 lastvalue = projch->GetBinCenter(k+1);
6119 //printf("Break by %f\n",lastvalue);
6122 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6126 //Some parameters to initialise
6127 Double_t widthLandau, widthGaus, mPV, integral;
6128 Double_t chisquarel = 0.0;
6129 Double_t chisquareg = 0.0;
6130 projch->Fit("landau","WWQ0M+",""
6131 ,(Double_t) mean/fBeginFitCharge
6133 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6134 chisquarel = projch->GetFunction("landau")->GetChisquare();
6135 projch->Fit("gaus","WWQ0M+",""
6136 ,(Double_t) mean/fBeginFitCharge
6138 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6139 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6141 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6142 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
6144 // Setting fit range and start values
6146 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
6147 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
6148 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
6149 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
6150 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
6151 fr[0] = mean/fBeginFitCharge;
6153 fCurrentCoef[0] = 0.0;
6154 fCurrentCoefE = 0.0;
6158 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
6163 //Double_t projchPeak;
6164 //Double_t projchFWHM;
6165 //LanGauPro(fp,projchPeak,projchFWHM);
6167 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6168 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6169 fNumberFitSuccess++;
6170 CalculChargeCoefMean(kTRUE);
6171 fCurrentCoef[0] = fp[1];
6172 fCurrentCoefE = fpe[1];
6173 //chargeCoefE2 = chisqr;
6176 CalculChargeCoefMean(kFALSE);
6177 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6179 if (fDebugLevel == 1) {
6180 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
6181 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6184 fitsnr->Draw("same");
6185 TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
6192 //_____________________________________________________________________________
6193 void AliTRDCalibraFit::FitBisCHEx(TH1* projch, Double_t mean, Double_t nentries)
6196 // Fit methode for the gain factor more time consuming
6199 //Calcul Range of the fit
6200 Double_t lastvalue = 0.0;
6201 Float_t sumAll = (Float_t) nentries;
6202 Int_t sumCurrent = 0;
6203 //printf("There are %d bins\n",nybins);
6204 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6205 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6206 if (fraction>fOutliersFitChargeHigh) {
6207 lastvalue = projch->GetBinCenter(k+1);
6208 //printf("Break by %f\n",lastvalue);
6211 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6216 //Some parameters to initialise
6217 Double_t widthLandau, widthGaus, mPV, integral;
6218 Double_t chisquarel = 0.0;
6219 Double_t chisquareg = 0.0;
6220 projch->Fit("landau","WWQM+",""
6221 ,(Double_t) mean/fBeginFitCharge
6223 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6224 chisquarel = projch->GetFunction("landau")->GetChisquare();
6225 projch->Fit("gaus","WWQM+",""
6226 ,(Double_t) mean/fBeginFitCharge
6228 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6229 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6231 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6232 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
6234 // Setting fit range and start values
6236 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
6237 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
6238 Double_t sv[5] = { widthLandau, mPV, integral, widthGaus, 0.0};
6239 Double_t pllo[5] = { 0.001, 0.001, projch->Integral()/3, 0.001, 0.0};
6240 Double_t plhi[5] = { 300.0, 300.0, 30*projch->Integral(), 300.0, 2.0};
6241 Double_t fp[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
6242 Double_t fpe[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
6244 //fr[0] = 0.3 * mean;
6245 //fr[1] = 3.0 * mean;
6247 fr[0] = mean/fBeginFitCharge;
6250 fCurrentCoef[0] = 0.0;
6251 fCurrentCoefE = 0.0;
6253 Double_t chisqr = 100.0;
6258 if((mPV > 0.0) && (projch->GetFunction("gaus")->GetParameter(1) > 0.0)) {
6259 fitsnr = LanGauFitEx(projch,&fr[0],&sv[0]
6265 //Double_t projchPeak;
6266 //Double_t projchFWHM;
6267 //LanGauProEx(fp,projchPeak,projchFWHM);
6269 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6270 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6271 fNumberFitSuccess++;
6272 CalculChargeCoefMean(kTRUE);
6273 fCurrentCoef[0] = fp[1];
6274 fCurrentCoefE = fpe[1];
6275 //chargeCoefE2 = chisqr;
6278 CalculChargeCoefMean(kFALSE);
6279 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6281 if (fDebugLevel == 1) {
6282 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
6283 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6286 if(fitsnr) fitsnr->Draw("same");
6287 TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
6294 //_____________________________________________________________________________
6295 void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6298 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
6300 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
6301 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
6302 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
6307 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
6308 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
6312 //_____________________________________________________________________________
6313 void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6316 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
6318 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
6319 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
6320 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
6321 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
6325 c2 = -(x0*(x[1]+x[2]+x[3])
6326 +x1*(x[0]+x[2]+x[3])
6327 +x2*(x[0]+x[1]+x[3])
6328 +x3*(x[0]+x[1]+x[2]));
6329 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
6330 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
6331 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
6332 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
6334 c0 = -(x0*x[1]*x[2]*x[3]
6337 +x3*x[0]*x[1]*x[2]);
6342 //_____________________________________________________________________________
6343 void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
6346 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
6348 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
6349 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
6350 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
6351 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
6352 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
6355 c4 = x0+x1+x2+x3+x4;
6356 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
6357 +x1*(x[0]+x[2]+x[3]+x[4])
6358 +x2*(x[0]+x[1]+x[3]+x[4])
6359 +x3*(x[0]+x[1]+x[2]+x[4])
6360 +x4*(x[0]+x[1]+x[2]+x[3]));
6361 c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
6362 +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])
6363 +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])
6364 +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])
6365 +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]));
6367 c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
6368 +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])
6369 +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])
6370 +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])
6371 +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]));
6373 c0 = (x0*x[1]*x[2]*x[3]*x[4]
6374 +x1*x[0]*x[2]*x[3]*x[4]
6375 +x2*x[0]*x[1]*x[3]*x[4]
6376 +x3*x[0]*x[1]*x[2]*x[4]
6377 +x4*x[0]*x[1]*x[2]*x[3]);
6380 //_____________________________________________________________________________
6381 void AliTRDCalibraFit::NormierungCharge()
6384 // Normalisation of the gain factor resulting for the fits
6387 // Calcul of the mean of choosen method by fFitChargeNDB
6389 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
6390 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
6392 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
6393 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
6394 //printf("detector %d coef[0] %f\n",detector,coef[0]);
6395 if (GetStack(detector) == 2) {
6398 if (GetStack(detector) != 2) {
6401 for (Int_t j = 0; j < total; j++) {
6409 fScaleFitFactor = fScaleFitFactor / sum;
6412 fScaleFitFactor = 1.0;
6415 //methode de boeuf mais bon...
6416 Double_t scalefactor = fScaleFitFactor;
6418 if(fDebugLevel > 1){
6420 if ( !fDebugStreamer ) {
6422 TDirectory *backup = gDirectory;
6423 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
6424 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
6426 (* fDebugStreamer) << "NormierungCharge"<<
6427 "scalefactor="<<scalefactor<<
6431 //_____________________________________________________________________________
6432 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
6435 // Rebin of the 1D histo for the gain calibration if needed.
6436 // you have to choose fRebin, divider of fNumberBinCharge
6439 const TAxis *xhist = hist->GetXaxis();
6440 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6441 ,xhist->GetBinLowEdge(1)
6442 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6444 AliInfo(Form("fRebin: %d",fRebin));
6446 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6448 for (Int_t ji = i; ji < i+fRebin; ji++) {
6449 sum += hist->GetBinContent(ji);
6452 rehist->SetBinContent(k,sum);
6460 //_____________________________________________________________________________
6461 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
6464 // Rebin of the 1D histo for the gain calibration if needed
6465 // you have to choose fRebin divider of fNumberBinCharge
6468 const TAxis *xhist = hist->GetXaxis();
6469 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6470 ,xhist->GetBinLowEdge(1)
6471 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6473 AliInfo(Form("fRebin: %d",fRebin));
6475 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6477 for (Int_t ji = i; ji < i+fRebin; ji++) {
6478 sum += hist->GetBinContent(ji);
6481 rehist->SetBinContent(k,sum);
6489 //____________Some basic geometry function_____________________________________
6492 //_____________________________________________________________________________
6493 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
6496 // Reconstruct the plane number from the detector number
6499 return ((Int_t) (d % 6));
6503 //_____________________________________________________________________________
6504 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
6507 // Reconstruct the stack number from the detector number
6509 const Int_t kNlayer = 6;
6511 return ((Int_t) (d % 30) / kNlayer);
6515 //_____________________________________________________________________________
6516 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6519 // Reconstruct the sector number from the detector number
6523 return ((Int_t) (d / fg));
6528 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6530 //_______________________________________________________________________________
6531 void AliTRDCalibraFit::ResetVectorFit()
6534 // Reset the VectorFits
6537 fVectorFit.SetOwner();
6539 fVectorFit2.SetOwner();
6540 fVectorFit2.Clear();
6544 //____________Private Functions________________________________________________
6547 //_____________________________________________________________________________
6548 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
6551 // Function for the fit
6554 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6556 //PARAMETERS FOR FIT PH
6558 //fAsymmGauss->SetParameter(0,0.113755);
6559 //fAsymmGauss->SetParameter(1,0.350706);
6560 //fAsymmGauss->SetParameter(2,0.0604244);
6561 //fAsymmGauss->SetParameter(3,7.65596);
6562 //fAsymmGauss->SetParameter(4,1.00124);
6563 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6571 Double_t dx = 0.005;
6572 Double_t xs = par[1];
6574 Double_t paras[2] = { 0.0, 0.0 };
6577 if ((xs >= par[1]) &&
6578 (xs < (par[1]+par[2]))) {
6579 //fAsymmGauss->SetParameter(0,par[0]);
6580 //fAsymmGauss->SetParameter(1,xs);
6581 //ss += fAsymmGauss->Eval(xx);
6584 ss += AsymmGauss(&xx,paras);
6586 if ((xs >= (par[1]+par[2])) &&
6587 (xs < (par[1]+par[2]+par[3]))) {
6588 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6589 //fAsymmGauss->SetParameter(1,xs);
6590 //ss += fAsymmGauss->Eval(xx);
6591 paras[0] = par[0]*par[4];
6593 ss += AsymmGauss(&xx,paras);
6602 //_____________________________________________________________________________
6603 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
6606 // Function for the fit
6609 //par[0] = normalization
6617 Double_t par1save = par[1];
6618 //Double_t par2save = par[2];
6619 Double_t par2save = 0.0604244;
6620 //Double_t par3save = par[3];
6621 Double_t par3save = 7.65596;
6622 //Double_t par5save = par[5];
6623 Double_t par5save = 0.870597;
6624 Double_t dx = x[0] - par1save;
6626 Double_t sigma2 = par2save*par2save;
6627 Double_t sqrt2 = TMath::Sqrt(2.0);
6628 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6629 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6630 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6631 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6633 //return par[0]*(exp1+par[4]*exp2);
6634 return par[0] * (exp1 + 1.00124 * exp2);
6638 //_____________________________________________________________________________
6639 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
6642 // Sum Landau + Gaus with identical mean
6645 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6646 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6647 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6648 Double_t val = valLandau + valGaus;
6654 //_____________________________________________________________________________
6655 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6658 // Function for the fit
6661 // par[0]=Width (scale) parameter of Landau density
6662 // par[1]=Most Probable (MP, location) parameter of Landau density
6663 // par[2]=Total area (integral -inf to inf, normalization constant)
6664 // par[3]=Width (sigma) of convoluted Gaussian function
6666 // In the Landau distribution (represented by the CERNLIB approximation),
6667 // the maximum is located at x=-0.22278298 with the location parameter=0.
6668 // This shift is corrected within this function, so that the actual
6669 // maximum is identical to the MP parameter.
6672 // Numeric constants
6673 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6674 Double_t mpshift = -0.22278298; // Landau maximum location
6676 // Control constants
6677 Double_t np = 100.0; // Number of convolution steps
6678 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6683 Double_t fland = 0.0;
6685 Double_t xlow = 0.0;
6686 Double_t xupp = 0.0;
6687 Double_t step = 0.0;
6690 // MP shift correction
6691 mpc = par[1] - mpshift * par[0];
6693 // Range of convolution integral
6694 xlow = x[0] - sc * par[3];
6695 xupp = x[0] + sc * par[3];
6697 step = (xupp - xlow) / np;
6699 // Convolution integral of Landau and Gaussian by sum
6700 for (i = 1.0; i <= np/2; i++) {
6702 xx = xlow + (i-.5) * step;
6703 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6704 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6706 xx = xupp - (i-.5) * step;
6707 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6708 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6712 if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
6716 //_____________________________________________________________________________
6717 Double_t AliTRDCalibraFit::LanGauFunEx(const Double_t *x, const Double_t *par)
6720 // Function for the fit
6723 // par[0]=Width (scale) parameter of Landau density
6724 // par[1]=Most Probable (MP, location) parameter of Landau density
6725 // par[2]=Total area (integral -inf to inf, normalization constant)
6726 // par[3]=Width (sigma) of convoluted Gaussian function
6727 // par[4]=Exponential Slope Parameter
6729 // In the Landau distribution (represented by the CERNLIB approximation),
6730 // the maximum is located at x=-0.22278298 with the location parameter=0.
6731 // This shift is corrected within this function, so that the actual
6732 // maximum is identical to the MP parameter.
6735 // Numeric constants
6736 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6737 Double_t mpshift = -0.22278298; // Landau maximum location
6739 // Control constants
6740 Double_t np = 100.0; // Number of convolution steps
6741 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6746 Double_t fland = 0.0;
6753 // MP shift correction
6754 mpc = par[1] - mpshift * par[0];
6756 // Range of convolution integral
6757 xlow = x[0] - sc * par[3];
6758 xupp = x[0] + sc * par[3];
6760 step = (xupp - xlow) / np;
6762 // Convolution integral of Landau and Gaussian by sum
6763 for (i = 1.0; i <= np/2; i++) {
6765 xx = xlow + (i-.5) * step;
6766 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
6767 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6769 xx = xupp - (i-.5) * step;
6770 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
6771 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6775 if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
6779 //_____________________________________________________________________________
6780 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6781 , const Double_t *parlimitslo, const Double_t *parlimitshi
6782 , Double_t *fitparams, Double_t *fiterrors
6783 , Double_t *chiSqr, Int_t *ndf) const
6786 // Function for the fit
6790 Char_t funname[100];
6792 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6797 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6798 ffit->SetParameters(startvalues);
6799 ffit->SetParNames("Width","MP","Area","GSigma");
6801 for (i = 0; i < 4; i++) {
6802 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6805 his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
6807 ffit->GetParameters(fitparams); // Obtain fit parameters
6808 for (i = 0; i < 4; i++) {
6809 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6811 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6812 ndf[0] = ffit->GetNDF(); // Obtain ndf
6814 return (ffit); // Return fit function
6817 //_____________________________________________________________________________
6818 TF1 *AliTRDCalibraFit::LanGauFitEx(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6819 , const Double_t *parlimitslo, const Double_t *parlimitshi
6820 , Double_t *fitparams, Double_t *fiterrors
6821 , Double_t *chiSqr, Int_t *ndf) const
6824 // Function for the fit
6828 Char_t funname[100];
6830 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6835 TF1 *ffit = new TF1(funname,LanGauFunEx,fitrange[0],fitrange[1],5);
6836 ffit->SetParameters(startvalues);
6837 ffit->SetParNames("Width","MP","Area","GSigma","Ex");
6839 for (i = 0; i < 5; i++) {
6840 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6843 his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
6845 ffit->GetParameters(fitparams); // Obtain fit parameters
6846 for (i = 0; i < 5; i++) {
6847 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6849 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6850 ndf[0] = ffit->GetNDF(); // Obtain ndf
6852 return (ffit); // Return fit function