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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
721 fStatisticMean = fStatisticMean / fNumberFit;
724 AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
726 delete fDebugStreamer;
727 fDebugStreamer = 0x0;
730 //________________functions fit Online 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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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: %d over these fitted histograms and %d successfulled fits"
910 ,(Int_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 TAxis *xprf = prf->GetXaxis();
936 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: %d over these fitted histograms and %d successfulled fits"
1006 ,(Int_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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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 fEntriesCurrent = 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 fEntriesCurrent = (Int_t) error[2];
1227 //printf("Number of entries %d\n",fEntriesCurrent);
1228 // Nothing found or not enough statistic
1229 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1230 NotEnoughStatisticLinearFitter();
1237 fStatisticMean += fEntriesCurrent;
1240 if((-(param[1])) <= 0.0) {
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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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 fEntriesCurrent = 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 fEntriesCurrent = (Int_t) error[2];
1310 //printf("Number of entries %d\n",fEntriesCurrent);
1311 // Nothing found or not enough statistic
1312 if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1313 NotEnoughStatisticExbAlt();
1320 fStatisticMean += fEntriesCurrent;
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: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_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);
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);
1424 linearfitter.AddPoint(&x,y);
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.0) 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 Float_t value = 0.0;
1731 /////////////////////////////////
1732 // Calculate the mean values
1733 ////////////////////////////////
1735 ////////////////////////
1736 Double_t meanAll = 0.0;
1737 Double_t rmsAll = 0.0;
1742 for (Int_t k = 0; k < loop; k++) {
1743 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1744 sector = GetSector(detector);
1746 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1748 rmsAll += value*value;
1754 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1755 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1756 for (Int_t row = 0; row < rowMax; row++) {
1757 for (Int_t col = 0; col < colMax; col++) {
1758 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1760 rmsAll += value*value;
1770 meanAll = meanAll/countAll;
1771 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1773 //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
1774 /////////////////////////////////////////////////
1776 ////////////////////////////////////////////////
1777 Double_t defaultvalue = -1.0;
1778 if(type==1) defaultvalue = -1.5;
1779 for (Int_t k = 0; k < loop; k++) {
1780 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1781 sector = GetSector(detector);
1782 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1783 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1784 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1786 // remove the results too far away
1787 for (Int_t row = 0; row < rowMax; row++) {
1788 for (Int_t col = 0; col < colMax; col++) {
1789 value = coef[(Int_t)(col*rowMax+row)];
1790 if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
1791 coef[(Int_t)(col*rowMax+row)] = defaultvalue;
1797 //______________________________________________________________________
1798 void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
1800 // Remove the results too far from the mean and rms
1804 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1806 AliInfo("The Vector Fit is not complete!");
1809 Int_t detector = -1;
1811 Float_t value = 0.0;
1813 /////////////////////////////////
1814 // Calculate the mean values
1815 ////////////////////////////////
1817 ////////////////////////
1818 Double_t meanAll = 0.0;
1819 Double_t rmsAll = 0.0;
1824 for (Int_t k = 0; k < loop; k++) {
1825 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1826 sector = GetSector(detector);
1828 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1831 rmsAll += value*value;
1836 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1837 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1838 for (Int_t row = 0; row < rowMax; row++) {
1839 for (Int_t col = 0; col < colMax; col++) {
1840 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1842 rmsAll += value*value;
1851 meanAll = meanAll/countAll;
1852 rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
1854 //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1855 /////////////////////////////////////////////////
1857 ////////////////////////////////////////////////
1858 for (Int_t k = 0; k < loop; k++) {
1859 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1860 sector = GetSector(detector);
1861 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1862 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1863 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1865 // remove the results too far away
1866 for (Int_t row = 0; row < rowMax; row++) {
1867 for (Int_t col = 0; col < colMax; col++) {
1868 value = coef[(Int_t)(col*rowMax+row)];
1869 if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) {
1870 //printf("value outlier %f\n",value);
1871 coef[(Int_t)(col*rowMax+row)] = 100.0;
1877 //______________________________________________________________________
1878 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1880 // ofwhat is equaled to 0: mean value of all passing detectors
1881 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1884 Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1886 AliInfo("The Vector Fit is not complete!");
1889 Int_t detector = -1;
1891 Float_t value = 0.0;
1893 /////////////////////////////////
1894 // Calculate the mean values
1895 ////////////////////////////////
1897 ////////////////////////
1898 Double_t meanAll = 0.0;
1899 Double_t meanSupermodule[18];
1900 Double_t meanDetector[540];
1901 Double_t rmsAll = 0.0;
1902 Double_t rmsSupermodule[18];
1903 Double_t rmsDetector[540];
1905 Int_t countSupermodule[18];
1906 Int_t countDetector[540];
1907 for(Int_t sm = 0; sm < 18; sm++){
1908 rmsSupermodule[sm] = 0.0;
1909 meanSupermodule[sm] = 0.0;
1910 countSupermodule[sm] = 0;
1912 for(Int_t det = 0; det < 540; det++){
1913 rmsDetector[det] = 0.0;
1914 meanDetector[det] = 0.0;
1915 countDetector[det] = 0;
1920 for (Int_t k = 0; k < loop; k++) {
1921 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1922 sector = GetSector(detector);
1924 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1926 rmsDetector[detector] += value*value;
1927 meanDetector[detector] += value;
1928 countDetector[detector]++;
1929 rmsSupermodule[sector] += value*value;
1930 meanSupermodule[sector] += value;
1931 countSupermodule[sector]++;
1932 rmsAll += value*value;
1938 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1939 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1940 for (Int_t row = 0; row < rowMax; row++) {
1941 for (Int_t col = 0; col < colMax; col++) {
1942 value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1944 rmsDetector[detector] += value*value;
1945 meanDetector[detector] += value;
1946 countDetector[detector]++;
1947 rmsSupermodule[sector] += value*value;
1948 meanSupermodule[sector] += value;
1949 countSupermodule[sector]++;
1950 rmsAll += value*value;
1960 meanAll = meanAll/countAll;
1961 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
1963 for(Int_t sm = 0; sm < 18; sm++){
1964 if(countSupermodule[sm] > 0) {
1965 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1966 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
1969 for(Int_t det = 0; det < 540; det++){
1970 if(countDetector[det] > 0) {
1971 meanDetector[det] = meanDetector[det]/countDetector[det];
1972 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
1975 //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
1976 ///////////////////////////////////////////////
1977 // Put the mean value for the no-fitted
1978 /////////////////////////////////////////////
1979 for (Int_t k = 0; k < loop; k++) {
1980 detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1981 sector = GetSector(detector);
1982 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1983 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
1984 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1986 for (Int_t row = 0; row < rowMax; row++) {
1987 for (Int_t col = 0; col < colMax; col++) {
1988 value = coef[(Int_t)(col*rowMax+row)];
1990 if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1992 if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1993 else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1994 else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1998 if(fDebugLevel > 1){
2000 if ( !fDebugStreamer ) {
2002 TDirectory *backup = gDirectory;
2003 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2004 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2007 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2009 (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
2010 "detector="<<detector<<
2022 //______________________________________________________________________
2023 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
2025 // ofwhat is equaled to 0: mean value of all passing detectors
2026 // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
2029 Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
2031 AliInfo("The Vector Fit is not complete!");
2034 Int_t detector = -1;
2036 Float_t value = 0.0;
2038 /////////////////////////////////
2039 // Calculate the mean values
2040 ////////////////////////////////
2042 ////////////////////////
2043 Double_t meanAll = 0.0;
2044 Double_t rmsAll = 0.0;
2045 Double_t meanSupermodule[18];
2046 Double_t rmsSupermodule[18];
2047 Double_t meanDetector[540];
2048 Double_t rmsDetector[540];
2050 Int_t countSupermodule[18];
2051 Int_t countDetector[540];
2052 for(Int_t sm = 0; sm < 18; sm++){
2053 rmsSupermodule[sm] = 0.0;
2054 meanSupermodule[sm] = 0.0;
2055 countSupermodule[sm] = 0;
2057 for(Int_t det = 0; det < 540; det++){
2058 rmsDetector[det] = 0.0;
2059 meanDetector[det] = 0.0;
2060 countDetector[det] = 0;
2064 for (Int_t k = 0; k < loop; k++) {
2065 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2066 sector = GetSector(detector);
2068 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
2070 rmsDetector[detector] += value*value;
2071 meanDetector[detector] += value;
2072 countDetector[detector]++;
2073 rmsSupermodule[sector] += value*value;
2074 meanSupermodule[sector] += value;
2075 countSupermodule[sector]++;
2077 rmsAll += value*value;
2082 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2083 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2084 for (Int_t row = 0; row < rowMax; row++) {
2085 for (Int_t col = 0; col < colMax; col++) {
2086 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2088 rmsDetector[detector] += value*value;
2089 meanDetector[detector] += value;
2090 countDetector[detector]++;
2091 rmsSupermodule[sector] += value*value;
2092 meanSupermodule[sector] += value;
2093 countSupermodule[sector]++;
2094 rmsAll += value*value;
2104 meanAll = meanAll/countAll;
2105 rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
2107 for(Int_t sm = 0; sm < 18; sm++){
2108 if(countSupermodule[sm] > 0) {
2109 meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
2110 rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
2113 for(Int_t det = 0; det < 540; det++){
2114 if(countDetector[det] > 0) {
2115 meanDetector[det] = meanDetector[det]/countDetector[det];
2116 rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
2119 //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
2120 ////////////////////////////////////////////
2121 // Put the mean value for the no-fitted
2122 /////////////////////////////////////////////
2123 for (Int_t k = 0; k < loop; k++) {
2124 detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
2125 sector = GetSector(detector);
2126 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2127 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2128 Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
2130 for (Int_t row = 0; row < rowMax; row++) {
2131 for (Int_t col = 0; col < colMax; col++) {
2132 value = coef[(Int_t)(col*rowMax+row)];
2134 if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2136 if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
2137 else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
2138 else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
2142 if(fDebugLevel > 1){
2144 if ( !fDebugStreamer ) {
2146 TDirectory *backup = gDirectory;
2147 fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2148 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2151 Float_t coefnow = coef[(Int_t)(col*rowMax+row)];
2153 (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
2154 "detector="<<detector<<
2167 //_____________________________________________________________________________
2168 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
2171 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2172 // It takes the mean value of the coefficients per detector
2173 // This object has to be written in the database
2176 // Create the DetObject
2177 AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2179 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2180 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2181 Int_t detector = -1;
2182 Float_t value = 0.0;
2185 for (Int_t k = 0; k < loop; k++) {
2186 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2189 mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
2193 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2194 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2195 for (Int_t row = 0; row < rowMax; row++) {
2196 for (Int_t col = 0; col < colMax; col++) {
2197 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2198 mean += TMath::Abs(value);
2202 if(count > 0) mean = mean/count;
2204 object->SetValue(detector,mean);
2209 //_____________________________________________________________________________
2210 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
2213 // It creates the AliTRDCalDet object from the AliTRDFitInfo
2214 // It takes the mean value of the coefficients per detector
2215 // This object has to be written in the database
2218 // Create the DetObject
2219 AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2221 fScaleGain = scaleFitFactor;
2223 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2224 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2225 Int_t detector = -1;
2226 Float_t value = 0.0;
2228 for (Int_t k = 0; k < loop; k++) {
2229 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2232 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2233 if(!meanOtherBefore){
2234 if(value > 0) value = value*scaleFitFactor;
2236 else value = value*scaleFitFactor;
2237 mean = TMath::Abs(value);
2241 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2242 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2243 for (Int_t row = 0; row < rowMax; row++) {
2244 for (Int_t col = 0; col < colMax; col++) {
2245 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2246 if(!meanOtherBefore) {
2247 if(value > 0) value = value*scaleFitFactor;
2249 else value = value*scaleFitFactor;
2250 mean += TMath::Abs(value);
2254 if(count > 0) mean = mean/count;
2256 if(mean < 0.1) mean = 0.1;
2257 object->SetValue(detector,mean);
2262 //_____________________________________________________________________________
2263 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
2266 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2267 // It takes the min value of the coefficients per detector
2268 // This object has to be written in the database
2271 // Create the DetObject
2272 AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2274 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2275 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2276 Int_t detector = -1;
2277 Float_t value = 0.0;
2279 for (Int_t k = 0; k < loop; k++) {
2280 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2281 Float_t min = 100.0;
2283 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2284 //printf("Create det object %f for %d\n",value,k);
2286 if(value > 70.0) value = value-100.0;
2291 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2292 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2293 for (Int_t row = 0; row < rowMax; row++) {
2294 for (Int_t col = 0; col < colMax; col++) {
2295 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2297 if(value > 70.0) value = value-100.0;
2299 if(min > value) min = value;
2303 object->SetValue(detector,min);
2309 //_____________________________________________________________________________
2310 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2313 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2314 // It takes the min value of the coefficients per detector
2315 // This object has to be written in the database
2318 // Create the DetObject
2319 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2322 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2323 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2324 Int_t detector = -1;
2325 Float_t value = 0.0;
2327 for (Int_t k = 0; k < loop; k++) {
2328 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2330 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2331 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2332 Float_t min = 100.0;
2333 for (Int_t row = 0; row < rowMax; row++) {
2334 for (Int_t col = 0; col < colMax; col++) {
2335 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2336 mean += -TMath::Abs(value);
2340 if(count > 0) mean = mean/count;
2342 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2343 if(value > 70.0) value = value-100.0;
2344 object->SetValue(detector,value);
2350 //_____________________________________________________________________________
2351 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit)
2354 // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2355 // It takes the min value of the coefficients per detector
2356 // This object has to be written in the database
2359 // Create the DetObject
2360 AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2363 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2364 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2365 Int_t detector = -1;
2366 Float_t value = 0.0;
2368 for (Int_t k = 0; k < loop; k++) {
2369 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2371 Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2372 Int_t colMax = fGeo->GetColMax(GetLayer(detector));
2373 Float_t min = 100.0;
2374 for (Int_t row = 0; row < rowMax; row++) {
2375 for (Int_t col = 0; col < colMax; col++) {
2376 value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2377 mean += -TMath::Abs(value);
2381 if(count > 0) mean = mean/count;
2383 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2384 //if(value > 70.0) value = value-100.0;
2385 object->SetValue(detector,value);
2391 //_____________________________________________________________________________
2392 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2395 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2396 // You need first to create the object for the detectors,
2397 // where the mean value is put.
2398 // This object has to be written in the database
2401 // Create the DetObject
2402 AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2405 for(Int_t k = 0; k < 540; k++){
2406 AliTRDCalROC *calROC = object->GetCalROC(k);
2407 Int_t nchannels = calROC->GetNchannels();
2408 for(Int_t ch = 0; ch < nchannels; ch++){
2409 calROC->SetValue(ch,1.0);
2415 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2416 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2417 Int_t detector = -1;
2418 Float_t value = 0.0;
2420 for (Int_t k = 0; k < loop; k++) {
2421 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2422 AliTRDCalROC *calROC = object->GetCalROC(detector);
2423 Float_t mean = detobject->GetValue(detector);
2424 if(TMath::Abs(mean) <= 0.0000000001) continue;
2425 Int_t rowMax = calROC->GetNrows();
2426 Int_t colMax = calROC->GetNcols();
2427 for (Int_t row = 0; row < rowMax; row++) {
2428 for (Int_t col = 0; col < colMax; col++) {
2429 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2430 if(value > 0) value = value*scaleFitFactor;
2431 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2439 //_____________________________________________________________________________
2440 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2443 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2444 // You need first to create the object for the detectors,
2445 // where the mean value is put.
2446 // This object has to be written in the database
2449 // Create the DetObject
2450 AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2453 for(Int_t k = 0; k < 540; k++){
2454 AliTRDCalROC *calROC = object->GetCalROC(k);
2455 Int_t nchannels = calROC->GetNchannels();
2456 for(Int_t ch = 0; ch < nchannels; ch++){
2457 calROC->SetValue(ch,1.0);
2463 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2464 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2465 Int_t detector = -1;
2466 Float_t value = 0.0;
2468 for (Int_t k = 0; k < loop; k++) {
2469 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2470 AliTRDCalROC *calROC = object->GetCalROC(detector);
2471 Float_t mean = detobject->GetValue(detector);
2472 if(mean == 0) continue;
2473 Int_t rowMax = calROC->GetNrows();
2474 Int_t colMax = calROC->GetNcols();
2475 for (Int_t row = 0; row < rowMax; row++) {
2476 for (Int_t col = 0; col < colMax; col++) {
2477 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2478 calROC->SetValue(col,row,TMath::Abs(value)/mean);
2486 //_____________________________________________________________________________
2487 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2490 // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2491 // You need first to create the object for the detectors,
2492 // where the mean value is put.
2493 // This object has to be written in the database
2496 // Create the DetObject
2497 AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2500 for(Int_t k = 0; k < 540; k++){
2501 AliTRDCalROC *calROC = object->GetCalROC(k);
2502 Int_t nchannels = calROC->GetNchannels();
2503 for(Int_t ch = 0; ch < nchannels; ch++){
2504 calROC->SetValue(ch,0.0);
2510 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2511 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2512 Int_t detector = -1;
2513 Float_t value = 0.0;
2515 for (Int_t k = 0; k < loop; k++) {
2516 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2517 AliTRDCalROC *calROC = object->GetCalROC(detector);
2518 Float_t min = detobject->GetValue(detector);
2519 Int_t rowMax = calROC->GetNrows();
2520 Int_t colMax = calROC->GetNcols();
2521 for (Int_t row = 0; row < rowMax; row++) {
2522 for (Int_t col = 0; col < colMax; col++) {
2523 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2525 if(value > 70.0) value = value - 100.0;
2527 calROC->SetValue(col,row,value-min);
2535 //_____________________________________________________________________________
2536 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2539 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2540 // This object has to be written in the database
2543 // Create the DetObject
2544 AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2546 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2547 if(loop != 540) AliInfo("The Vector Fit is not complete!");
2548 Int_t detector = -1;
2549 Float_t value = 0.0;
2551 for (Int_t k = 0; k < loop; k++) {
2552 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2553 AliTRDCalROC *calROC = object->GetCalROC(detector);
2554 Int_t rowMax = calROC->GetNrows();
2555 Int_t colMax = calROC->GetNcols();
2556 for (Int_t row = 0; row < rowMax; row++) {
2557 for (Int_t col = 0; col < colMax; col++) {
2558 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2559 calROC->SetValue(col,row,TMath::Abs(value));
2567 //_____________________________________________________________________________
2568 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2571 // It Creates the AliTRDCalDet object from AliTRDFitInfo
2572 // 0 successful fit 1 not successful fit
2573 // mean is the mean value over the successful fit
2574 // do not use it for t0: no meaning
2577 // Create the CalObject
2578 AliTRDCalDet *object = new AliTRDCalDet(name,name);
2582 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2584 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2585 for(Int_t k = 0; k < 540; k++){
2586 object->SetValue(k,1.0);
2589 Int_t detector = -1;
2590 Float_t value = 0.0;
2592 for (Int_t k = 0; k < loop; k++) {
2593 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2594 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2595 if(value <= 0) object->SetValue(detector,1.0);
2597 object->SetValue(detector,0.0);
2602 if(count > 0) mean /= count;
2605 //_____________________________________________________________________________
2606 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2609 // It Creates the AliTRDCalPad object from AliTRDFitInfo
2610 // 0 not successful fit 1 successful fit
2611 // mean mean value over the successful fit
2614 // Create the CalObject
2615 AliTRDCalPad *object = new AliTRDCalPad(name,name);
2619 Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2621 AliInfo("The Vector Fit is not complete! We initialise all outliers");
2622 for(Int_t k = 0; k < 540; k++){
2623 AliTRDCalROC *calROC = object->GetCalROC(k);
2624 Int_t nchannels = calROC->GetNchannels();
2625 for(Int_t ch = 0; ch < nchannels; ch++){
2626 calROC->SetValue(ch,1.0);
2630 Int_t detector = -1;
2631 Float_t value = 0.0;
2633 for (Int_t k = 0; k < loop; k++) {
2634 detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2635 AliTRDCalROC *calROC = object->GetCalROC(detector);
2636 Int_t nchannels = calROC->GetNchannels();
2637 for (Int_t ch = 0; ch < nchannels; ch++) {
2638 value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2639 if(value <= 0) calROC->SetValue(ch,1.0);
2641 calROC->SetValue(ch,0.0);
2647 if(count > 0) mean /= count;
2650 //_____________________________________________________________________________
2651 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2654 // Set FitPH if 1 then each detector will be fitted
2657 if (periodeFitPH > 0) {
2658 fFitPHPeriode = periodeFitPH;
2661 AliInfo("periodeFitPH must be higher than 0!");
2665 //_____________________________________________________________________________
2666 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2669 // The fit of the deposited charge distribution begins at
2670 // histo->Mean()/beginFitCharge
2671 // You can here set beginFitCharge
2674 if (beginFitCharge > 0) {
2675 fBeginFitCharge = beginFitCharge;
2678 AliInfo("beginFitCharge must be strict positif!");
2683 //_____________________________________________________________________________
2684 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2687 // The t0 calculated with the maximum positif slope is shift from t0Shift0
2688 // You can here set t0Shift0
2692 fT0Shift0 = t0Shift;
2695 AliInfo("t0Shift0 must be strict positif!");
2700 //_____________________________________________________________________________
2701 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2704 // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2705 // You can here set t0Shift1
2709 fT0Shift1 = t0Shift;
2712 AliInfo("t0Shift must be strict positif!");
2717 //_____________________________________________________________________________
2718 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2721 // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2722 // You can here set rangeFitPRF
2725 if ((rangeFitPRF > 0) &&
2726 (rangeFitPRF <= 1.5)) {
2727 fRangeFitPRF = rangeFitPRF;
2730 AliInfo("rangeFitPRF must be between 0 and 1.0");
2735 //_____________________________________________________________________________
2736 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2739 // Minimum entries for fitting
2742 if (minEntries > 0) {
2743 fMinEntries = minEntries;
2746 AliInfo("fMinEntries must be >= 0.");
2751 //_____________________________________________________________________________
2752 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2755 // Rebin with rebin time less bins the Ch histo
2756 // You can set here rebin that should divide the number of bins of CH histo
2761 AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2764 AliInfo("You have to choose a positiv value!");
2768 //_____________________________________________________________________________
2769 Bool_t AliTRDCalibraFit::FillVectorFit()
2772 // For the Fit functions fill the vector Fit
2775 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2778 if (GetStack(fCountDet) == 2) {
2785 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2786 Float_t *coef = new Float_t[ntotal];
2787 for (Int_t i = 0; i < ntotal; i++) {
2788 coef[i] = fCurrentCoefDetector[i];
2791 Int_t detector = fCountDet;
2793 fitInfo->SetCoef(coef);
2794 fitInfo->SetDetector(detector);
2795 fVectorFit.Add((TObject *) fitInfo);
2800 //_____________________________________________________________________________
2801 Bool_t AliTRDCalibraFit::FillVectorFit2()
2804 // For the Fit functions fill the vector Fit
2807 AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2810 if (GetStack(fCountDet) == 2) {
2817 //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2818 Float_t *coef = new Float_t[ntotal];
2819 for (Int_t i = 0; i < ntotal; i++) {
2820 coef[i] = fCurrentCoefDetector2[i];
2823 Int_t detector = fCountDet;
2825 fitInfo->SetCoef(coef);
2826 fitInfo->SetDetector(detector);
2827 fVectorFit2.Add((TObject *) fitInfo);
2832 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2833 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2836 // Init the number of expected bins and fDect1[i] fDect2[i]
2839 gStyle->SetPalette(1);
2840 gStyle->SetOptStat(1111);
2841 gStyle->SetPadBorderMode(0);
2842 gStyle->SetCanvasColor(10);
2843 gStyle->SetPadLeftMargin(0.13);
2844 gStyle->SetPadRightMargin(0.01);
2846 // Mode groups of pads: the total number of bins!
2847 CalculNumberOfBinsExpected(i);
2849 // Quick verification that we have the good pad calibration mode!
2850 if (fNumberOfBinsExpected != nbins) {
2851 AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2855 // Security for fDebug 3 and 4
2856 if ((fDebugLevel >= 3) &&
2860 AliInfo("This detector doesn't exit!");
2864 // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2865 CalculDect1Dect2(i);
2870 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2871 Bool_t AliTRDCalibraFit::InitFitCH()
2874 // Init the fVectorFitCH for normalisation
2875 // Init the histo for debugging
2880 fScaleFitFactor = 0.0;
2881 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2882 fCurrentCoefDetector = new Float_t[2304];
2883 for (Int_t k = 0; k < 2304; k++) {
2884 fCurrentCoefDetector[k] = 0.0;
2886 fVectorFit.SetName("gainfactorscoefficients");
2888 // fDebug == 0 nothing
2889 // fDebug == 1 and fFitVoir no histo
2890 if (fDebugLevel == 1) {
2891 if(!CheckFitVoir()) return kFALSE;
2893 //Get the CalDet object
2895 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2897 AliInfo("Could not get calibDB");
2900 if(fCalDet) delete fCalDet;
2901 fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2904 Float_t devalue = 1.0;
2905 if(fCalDet) delete fCalDet;
2906 fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2907 for(Int_t k = 0; k < 540; k++){
2908 fCalDet->SetValue(k,devalue);
2914 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2915 Bool_t AliTRDCalibraFit::InitFitPH()
2918 // Init the arrays of results
2919 // Init the histos for debugging
2923 fVectorFit.SetName("driftvelocitycoefficients");
2924 fVectorFit2.SetName("t0coefficients");
2926 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2927 fCurrentCoefDetector = new Float_t[2304];
2928 for (Int_t k = 0; k < 2304; k++) {
2929 fCurrentCoefDetector[k] = 0.0;
2931 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
2932 fCurrentCoefDetector2 = new Float_t[2304];
2933 for (Int_t k = 0; k < 2304; k++) {
2934 fCurrentCoefDetector2[k] = 0.0;
2937 //fDebug == 0 nothing
2938 // fDebug == 1 and fFitVoir no histo
2939 if (fDebugLevel == 1) {
2940 if(!CheckFitVoir()) return kFALSE;
2942 //Get the CalDet object
2944 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2946 AliInfo("Could not get calibDB");
2949 if(fCalDet) delete fCalDet;
2950 if(fCalDet2) delete fCalDet2;
2951 fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
2952 fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2955 Float_t devalue = 1.5;
2956 Float_t devalue2 = 0.0;
2957 if(fCalDet) delete fCalDet;
2958 if(fCalDet2) delete fCalDet2;
2959 fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2960 fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2961 for(Int_t k = 0; k < 540; k++){
2962 fCalDet->SetValue(k,devalue);
2963 fCalDet2->SetValue(k,devalue2);
2968 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2969 Bool_t AliTRDCalibraFit::InitFitPRF()
2972 // Init the calibration mode (Nz, Nrphi), the histograms for
2973 // debugging the fit methods if fDebug > 0,
2977 fVectorFit.SetName("prfwidthcoefficients");
2979 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
2980 fCurrentCoefDetector = new Float_t[2304];
2981 for (Int_t k = 0; k < 2304; k++) {
2982 fCurrentCoefDetector[k] = 0.0;
2985 // fDebug == 0 nothing
2986 // fDebug == 1 and fFitVoir no histo
2987 if (fDebugLevel == 1) {
2988 if(!CheckFitVoir()) return kFALSE;
2992 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2993 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2996 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3001 if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
3002 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3003 fCurrentCoefDetector = new Float_t[2304];
3004 fCurrentCoefDetector2 = new Float_t[2304];
3005 for (Int_t k = 0; k < 2304; k++) {
3006 fCurrentCoefDetector[k] = 0.0;
3007 fCurrentCoefDetector2[k] = 0.0;
3010 if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) return kFALSE;
3014 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3015 Bool_t AliTRDCalibraFit::InitFitExbAlt()
3018 // Init the fCalDet, fVectorFit fCurrentCoefDetector
3023 if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
3024 fCurrentCoefDetector2 = new Float_t[2304];
3025 for (Int_t k = 0; k < 2304; k++) {
3026 fCurrentCoefDetector2[k] = 0.0;
3031 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3032 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
3035 // Init the current detector where we are fCountDet and the
3036 // next fCount for the functions Fit...
3039 // Loop on the Xbins of ch!!
3040 fCountDet = -1; // Current detector
3041 fCount = 0; // To find the next detector
3044 if (fDebugLevel >= 3) {
3045 // Set countdet to the detector
3046 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3047 // Set counter to write at the end of the detector
3049 // Get the right calib objects
3052 if(fDebugLevel == 1) {
3054 fCalibraMode->CalculXBins(fCountDet,i);
3055 if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
3056 while(fCalibraMode->GetXbins(i) <=fFitVoir){
3058 fCalibraMode->CalculXBins(fCountDet,i);
3059 //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
3065 fCount = fCalibraMode->GetXbins(i);
3067 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3068 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3069 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3070 ,(Int_t) GetStack(fCountDet)
3071 ,(Int_t) GetSector(fCountDet),i);
3074 //_______________________________________________________________________________
3075 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
3078 // Calculate the number of bins expected (calibration groups)
3081 fNumberOfBinsExpected = 0;
3083 if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
3084 fNumberOfBinsExpected = 1;
3088 if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
3089 fNumberOfBinsExpected = 18;
3093 fCalibraMode->ModePadCalibration(2,i);
3094 fCalibraMode->ModePadFragmentation(0,2,0,i);
3095 fCalibraMode->SetDetChamb2(i);
3096 if (fDebugLevel > 1) {
3097 AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
3099 fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
3100 fCalibraMode->ModePadCalibration(0,i);
3101 fCalibraMode->ModePadFragmentation(0,0,0,i);
3102 fCalibraMode->SetDetChamb0(i);
3103 if (fDebugLevel > 1) {
3104 AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
3106 fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
3109 //_______________________________________________________________________________
3110 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
3113 // Calculate the range of fits
3118 if (fDebugLevel == 1) {
3122 if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
3124 fDect2 = fNumberOfBinsExpected;
3126 if (fDebugLevel >= 3) {
3127 fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
3128 fCalibraMode->CalculXBins(fCountDet,i);
3129 fDect1 = fCalibraMode->GetXbins(i);
3130 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3131 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3132 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3133 ,(Int_t) GetStack(fCountDet)
3134 ,(Int_t) GetSector(fCountDet),i);
3135 // Set for the next detector
3136 fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3139 //_______________________________________________________________________________
3140 Bool_t AliTRDCalibraFit::CheckFitVoir()
3143 // Check if fFitVoir is in the range
3146 if (fFitVoir < fNumberOfBinsExpected) {
3147 AliInfo(Form("We will see the fit of the object %d",fFitVoir));
3150 AliInfo("fFitVoir is out of range of the histo!");
3155 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3156 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
3159 // See if we are in a new detector and update the
3160 // variables fNfragZ and fNfragRphi if yes
3161 // Will never happen for only one detector (3 and 4)
3162 // Doesn't matter for 2
3164 if (fCount == idect) {
3165 // On en est au detector (or first detector in the group)
3167 AliDebug(2,Form("We are at the detector %d\n",fCountDet));
3168 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3169 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
3170 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3171 ,(Int_t) GetStack(fCountDet)
3172 ,(Int_t) GetSector(fCountDet),i);
3173 // Set for the next detector
3174 fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
3179 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3180 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
3183 // Reconstruct the min pad row, max pad row, min pad col and
3184 // max pad col of the calibration group for the Fit functions
3185 // idect is the calibration group inside the detector
3187 if (fDebugLevel != 1) {
3188 fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
3190 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
3191 AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
3193 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3194 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
3197 // For the case where there are not enough entries in the histograms
3198 // of the calibration group, the value present in the choosen database
3199 // will be put. A negativ sign enables to know that a fit was not possible.
3202 if (fDebugLevel == 1) {
3203 AliInfo("The element has not enough statistic to be fitted");
3205 else if (fNbDet > 0){
3206 Int_t firstdetector = fCountDet;
3207 Int_t lastdetector = fCountDet+fNbDet;
3208 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3209 // loop over detectors
3210 for(Int_t det = firstdetector; det < lastdetector; det++){
3212 //Set the calibration object again
3216 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3218 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3219 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3220 ,(Int_t) GetStack(fCountDet)
3221 ,(Int_t) GetSector(fCountDet),0);
3222 // Reconstruct row min row max
3223 ReconstructFitRowMinRowMax(idect,0);
3225 // Calcul the coef from the database choosen for the detector
3226 CalculChargeCoefMean(kFALSE);
3228 //stack 2, not stack 2
3230 if(GetStack(fCountDet) == 2) factor = 12;
3233 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3234 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3235 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3236 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3240 //Put default value negative
3241 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3242 fCurrentCoefE = 0.0;
3247 if(fDebugLevel > 1){
3249 if ( !fDebugStreamer ) {
3251 TDirectory *backup = gDirectory;
3252 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3253 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3256 Int_t detector = fCountDet;
3257 Int_t caligroup = idect;
3258 Short_t rowmin = fCalibraMode->GetRowMin(0);
3259 Short_t rowmax = fCalibraMode->GetRowMax(0);
3260 Short_t colmin = fCalibraMode->GetColMin(0);
3261 Short_t colmax = fCalibraMode->GetColMax(0);
3262 Float_t gf = fCurrentCoef[0];
3263 Float_t gfs = fCurrentCoef[1];
3264 Float_t gfE = fCurrentCoefE;
3266 (*fDebugStreamer) << "FillFillCH" <<
3267 "detector=" << detector <<
3268 "caligroup=" << caligroup <<
3269 "rowmin=" << rowmin <<
3270 "rowmax=" << rowmax <<
3271 "colmin=" << colmin <<
3272 "colmax=" << colmax <<
3280 for (Int_t k = 0; k < 2304; k++) {
3281 fCurrentCoefDetector[k] = 0.0;
3285 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3289 //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));
3291 // Calcul the coef from the database choosen
3292 CalculChargeCoefMean(kFALSE);
3294 //stack 2, not stack 2
3296 if(GetStack(fCountDet) == 2) factor = 12;
3299 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3300 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3301 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3302 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3306 //Put default value negative
3307 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3308 fCurrentCoefE = 0.0;
3317 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3318 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
3321 // For the case where there are not enough entries in the histograms
3322 // of the calibration group, the value present in the choosen database
3323 // will be put. A negativ sign enables to know that a fit was not possible.
3325 if (fDebugLevel == 1) {
3326 AliInfo("The element has not enough statistic to be fitted");
3328 else if (fNbDet > 0) {
3330 Int_t firstdetector = fCountDet;
3331 Int_t lastdetector = fCountDet+fNbDet;
3332 //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3333 // loop over detectors
3334 for(Int_t det = firstdetector; det < lastdetector; det++){
3336 //Set the calibration object again
3340 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3342 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3343 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3344 ,(Int_t) GetStack(fCountDet)
3345 ,(Int_t) GetSector(fCountDet),1);
3346 // Reconstruct row min row max
3347 ReconstructFitRowMinRowMax(idect,1);
3349 // Calcul the coef from the database choosen for the detector
3350 CalculVdriftCoefMean();
3353 //stack 2, not stack 2
3355 if(GetStack(fCountDet) == 2) factor = 12;
3358 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3359 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3360 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3361 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3362 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3366 //Put default value negative
3367 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3368 fCurrentCoefE = 0.0;
3369 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3370 fCurrentCoefE2 = 0.0;
3376 if(fDebugLevel > 1){
3378 if ( !fDebugStreamer ) {
3380 TDirectory *backup = gDirectory;
3381 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3382 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3386 Int_t detector = fCountDet;
3387 Int_t caligroup = idect;
3388 Short_t rowmin = fCalibraMode->GetRowMin(1);
3389 Short_t rowmax = fCalibraMode->GetRowMax(1);
3390 Short_t colmin = fCalibraMode->GetColMin(1);
3391 Short_t colmax = fCalibraMode->GetColMax(1);
3392 Float_t vf = fCurrentCoef[0];
3393 Float_t vs = fCurrentCoef[1];
3394 Float_t vfE = fCurrentCoefE;
3395 Float_t t0f = fCurrentCoef2[0];
3396 Float_t t0s = fCurrentCoef2[1];
3397 Float_t t0E = fCurrentCoefE2;
3401 (* fDebugStreamer) << "FillFillPH"<<
3402 "detector="<<detector<<
3403 "nentries="<<nentries<<
3404 "caligroup="<<caligroup<<
3418 for (Int_t k = 0; k < 2304; k++) {
3419 fCurrentCoefDetector[k] = 0.0;
3420 fCurrentCoefDetector2[k] = 0.0;
3424 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3428 //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));
3430 CalculVdriftCoefMean();
3433 //stack 2 and not stack 2
3435 if(GetStack(fCountDet) == 2) factor = 12;
3439 // Fill the fCurrentCoefDetector 2
3440 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3441 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3442 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3443 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3447 // Put the default value
3448 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3449 fCurrentCoefE = 0.0;
3450 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3451 fCurrentCoefE2 = 0.0;
3453 FillFillPH(idect,nentries);
3462 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3463 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3466 // For the case where there are not enough entries in the histograms
3467 // of the calibration group, the value present in the choosen database
3468 // will be put. A negativ sign enables to know that a fit was not possible.
3471 if (fDebugLevel == 1) {
3472 AliInfo("The element has not enough statistic to be fitted");
3474 else if (fNbDet > 0){
3476 Int_t firstdetector = fCountDet;
3477 Int_t lastdetector = fCountDet+fNbDet;
3478 // AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
3480 // loop over detectors
3481 for(Int_t det = firstdetector; det < lastdetector; det++){
3483 //Set the calibration object again
3487 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3489 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3490 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3491 ,(Int_t) GetStack(fCountDet)
3492 ,(Int_t) GetSector(fCountDet),2);
3493 // Reconstruct row min row max
3494 ReconstructFitRowMinRowMax(idect,2);
3496 // Calcul the coef from the database choosen for the detector
3497 CalculPRFCoefMean();
3499 //stack 2, not stack 2
3501 if(GetStack(fCountDet) == 2) factor = 12;
3504 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3505 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3506 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3507 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3511 //Put default value negative
3512 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3513 fCurrentCoefE = 0.0;
3518 if(fDebugLevel > 1){
3520 if ( !fDebugStreamer ) {
3522 TDirectory *backup = gDirectory;
3523 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3524 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3527 Int_t detector = fCountDet;
3528 Int_t layer = GetLayer(fCountDet);
3529 Int_t caligroup = idect;
3530 Short_t rowmin = fCalibraMode->GetRowMin(2);
3531 Short_t rowmax = fCalibraMode->GetRowMax(2);
3532 Short_t colmin = fCalibraMode->GetColMin(2);
3533 Short_t colmax = fCalibraMode->GetColMax(2);
3534 Float_t widf = fCurrentCoef[0];
3535 Float_t wids = fCurrentCoef[1];
3536 Float_t widfE = fCurrentCoefE;
3538 (* fDebugStreamer) << "FillFillPRF"<<
3539 "detector="<<detector<<
3541 "caligroup="<<caligroup<<
3552 for (Int_t k = 0; k < 2304; k++) {
3553 fCurrentCoefDetector[k] = 0.0;
3557 AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3561 // 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));
3563 CalculPRFCoefMean();
3565 // stack 2 and not stack 2
3567 if(GetStack(fCountDet) == 2) factor = 12;
3571 // Fill the fCurrentCoefDetector
3572 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3573 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3574 fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3578 // Put the default value
3579 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3580 fCurrentCoefE = 0.0;
3588 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3589 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3592 // For the case where there are not enough entries in the histograms
3593 // of the calibration group, the value present in the choosen database
3594 // will be put. A negativ sign enables to know that a fit was not possible.
3597 // Calcul the coef from the database choosen
3598 CalculVdriftLorentzCoef();
3601 if(GetStack(fCountDet) == 2) factor = 1728;
3605 // Fill the fCurrentCoefDetector
3606 for (Int_t k = 0; k < factor; k++) {
3607 fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3608 // should be negative
3609 fCurrentCoefDetector2[k] = fCurrentCoef2[1]+100.0;
3613 //Put default opposite sign only for vdrift
3614 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3615 fCurrentCoefE = 0.0;
3616 fCurrentCoef2[0] = fCurrentCoef2[1]+100.0;
3617 fCurrentCoefE2 = 0.0;
3619 FillFillLinearFitter();
3624 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3625 Bool_t AliTRDCalibraFit::NotEnoughStatisticExbAlt()
3628 // For the case where there are not enough entries in the histograms
3629 // of the calibration group, the value present in the choosen database
3630 // will be put. A negativ sign enables to know that a fit was not possible.
3634 if(GetStack(fCountDet) == 2) factor = 1728;
3638 // Fill the fCurrentCoefDetector
3639 for (Int_t k = 0; k < factor; k++) {
3640 fCurrentCoefDetector2[k] = 100.0;
3643 fCurrentCoef2[0] = 100.0;
3644 fCurrentCoefE2 = 0.0;
3651 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3652 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3655 // Fill the coefficients found with the fits or other
3656 // methods from the Fit functions
3659 if (fDebugLevel != 1) {
3661 Int_t firstdetector = fCountDet;
3662 Int_t lastdetector = fCountDet+fNbDet;
3663 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3664 // loop over detectors
3665 for(Int_t det = firstdetector; det < lastdetector; det++){
3667 //Set the calibration object again
3671 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3673 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3674 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3675 ,(Int_t) GetStack(fCountDet)
3676 ,(Int_t) GetSector(fCountDet),0);
3677 // Reconstruct row min row max
3678 ReconstructFitRowMinRowMax(idect,0);
3680 // Calcul the coef from the database choosen for the detector
3681 if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3682 else CalculChargeCoefMean(kTRUE);
3684 //stack 2, not stack 2
3686 if(GetStack(fCountDet) == 2) factor = 12;
3689 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3690 Double_t coeftoput = 1.0;
3691 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3692 else coeftoput = fCurrentCoef[0];
3693 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3694 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3695 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3702 if(fDebugLevel > 1){
3704 if ( !fDebugStreamer ) {
3706 TDirectory *backup = gDirectory;
3707 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3708 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3711 Int_t detector = fCountDet;
3712 Int_t caligroup = idect;
3713 Short_t rowmin = fCalibraMode->GetRowMin(0);
3714 Short_t rowmax = fCalibraMode->GetRowMax(0);
3715 Short_t colmin = fCalibraMode->GetColMin(0);
3716 Short_t colmax = fCalibraMode->GetColMax(0);
3717 Float_t gf = fCurrentCoef[0];
3718 Float_t gfs = fCurrentCoef[1];
3719 Float_t gfE = fCurrentCoefE;
3721 (*fDebugStreamer) << "FillFillCH" <<
3722 "detector=" << detector <<
3723 "caligroup=" << caligroup <<
3724 "rowmin=" << rowmin <<
3725 "rowmax=" << rowmax <<
3726 "colmin=" << colmin <<
3727 "colmax=" << colmax <<
3735 for (Int_t k = 0; k < 2304; k++) {
3736 fCurrentCoefDetector[k] = 0.0;
3740 //printf("Check the count now: fCountDet %d\n",fCountDet);
3745 if(GetStack(fCountDet) == 2) factor = 12;
3748 for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3749 for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3750 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3761 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3762 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3765 // Fill the coefficients found with the fits or other
3766 // methods from the Fit functions
3769 if (fDebugLevel != 1) {
3772 Int_t firstdetector = fCountDet;
3773 Int_t lastdetector = fCountDet+fNbDet;
3774 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3776 // loop over detectors
3777 for(Int_t det = firstdetector; det < lastdetector; det++){
3779 //Set the calibration object again
3783 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3785 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3786 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3787 ,(Int_t) GetStack(fCountDet)
3788 ,(Int_t) GetSector(fCountDet),1);
3789 // Reconstruct row min row max
3790 ReconstructFitRowMinRowMax(idect,1);
3792 // Calcul the coef from the database choosen for the detector
3793 CalculVdriftCoefMean();
3796 //stack 2, not stack 2
3798 if(GetStack(fCountDet) == 2) factor = 12;
3801 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3802 Double_t coeftoput = 1.5;
3803 Double_t coeftoput2 = 0.0;
3805 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3806 else coeftoput = fCurrentCoef[0];
3808 if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3809 else coeftoput2 = fCurrentCoef2[0];
3811 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3812 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3813 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3814 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3822 if(fDebugLevel > 1){
3824 if ( !fDebugStreamer ) {
3826 TDirectory *backup = gDirectory;
3827 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3828 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3832 Int_t detector = fCountDet;
3833 Int_t caligroup = idect;
3834 Short_t rowmin = fCalibraMode->GetRowMin(1);
3835 Short_t rowmax = fCalibraMode->GetRowMax(1);
3836 Short_t colmin = fCalibraMode->GetColMin(1);
3837 Short_t colmax = fCalibraMode->GetColMax(1);
3838 Float_t vf = fCurrentCoef[0];
3839 Float_t vs = fCurrentCoef[1];
3840 Float_t vfE = fCurrentCoefE;
3841 Float_t t0f = fCurrentCoef2[0];
3842 Float_t t0s = fCurrentCoef2[1];
3843 Float_t t0E = fCurrentCoefE2;
3847 (* fDebugStreamer) << "FillFillPH"<<
3848 "detector="<<detector<<
3849 "nentries="<<nentries<<
3850 "caligroup="<<caligroup<<
3864 for (Int_t k = 0; k < 2304; k++) {
3865 fCurrentCoefDetector[k] = 0.0;
3866 fCurrentCoefDetector2[k] = 0.0;
3870 //printf("Check the count now: fCountDet %d\n",fCountDet);
3875 if(GetStack(fCountDet) == 2) factor = 12;
3878 for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3879 for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3880 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3881 fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3885 FillFillPH(idect,nentries);
3890 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3891 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3894 // Fill the coefficients found with the fits or other
3895 // methods from the Fit functions
3898 if (fDebugLevel != 1) {
3901 Int_t firstdetector = fCountDet;
3902 Int_t lastdetector = fCountDet+fNbDet;
3903 // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
3905 // loop over detectors
3906 for(Int_t det = firstdetector; det < lastdetector; det++){
3908 //Set the calibration object again
3912 // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3914 fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3915 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3916 ,(Int_t) GetStack(fCountDet)
3917 ,(Int_t) GetSector(fCountDet),2);
3918 // Reconstruct row min row max
3919 ReconstructFitRowMinRowMax(idect,2);
3921 // Calcul the coef from the database choosen for the detector
3922 CalculPRFCoefMean();
3924 //stack 2, not stack 2
3926 if(GetStack(fCountDet) == 2) factor = 12;
3929 // Fill the fCurrentCoefDetector with negative value to say: not fitted
3930 Double_t coeftoput = 1.0;
3931 if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3932 else coeftoput = fCurrentCoef[0];
3933 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3934 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3935 fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3942 if(fDebugLevel > 1){
3944 if ( !fDebugStreamer ) {
3946 TDirectory *backup = gDirectory;
3947 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3948 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3951 Int_t detector = fCountDet;
3952 Int_t layer = GetLayer(fCountDet);
3953 Int_t caligroup = idect;
3954 Short_t rowmin = fCalibraMode->GetRowMin(2);
3955 Short_t rowmax = fCalibraMode->GetRowMax(2);
3956 Short_t colmin = fCalibraMode->GetColMin(2);
3957 Short_t colmax = fCalibraMode->GetColMax(2);
3958 Float_t widf = fCurrentCoef[0];
3959 Float_t wids = fCurrentCoef[1];
3960 Float_t widfE = fCurrentCoefE;
3962 (* fDebugStreamer) << "FillFillPRF"<<
3963 "detector="<<detector<<
3965 "caligroup="<<caligroup<<
3976 for (Int_t k = 0; k < 2304; k++) {
3977 fCurrentCoefDetector[k] = 0.0;
3981 //printf("Check the count now: fCountDet %d\n",fCountDet);
3986 if(GetStack(fCountDet) == 2) factor = 12;
3989 // Pointer to the branch
3990 for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3991 for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3992 fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
4002 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4003 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
4006 // Fill the coefficients found with the fits or other
4007 // methods from the Fit functions
4011 if(GetStack(fCountDet) == 2) factor = 1728;
4014 // Pointer to the branch
4015 for (Int_t k = 0; k < factor; k++) {
4016 fCurrentCoefDetector[k] = fCurrentCoef[0];
4017 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4020 FillFillLinearFitter();
4025 //____________Functions for initialising the AliTRDCalibraFit in the code_________
4026 Bool_t AliTRDCalibraFit::FillInfosFitExbAlt()
4029 // Fill the coefficients found with the fits or other
4030 // methods from the Fit functions
4034 if(GetStack(fCountDet) == 2) factor = 1728;
4037 // Pointer to the branch
4038 for (Int_t k = 0; k < factor; k++) {
4039 fCurrentCoefDetector2[k] = fCurrentCoef2[0];
4047 //________________________________________________________________________________
4048 void AliTRDCalibraFit::FillFillCH(Int_t idect)
4051 // DebugStream and fVectorFit
4054 // End of one detector
4055 if (idect == (fCount-1)) {
4058 for (Int_t k = 0; k < 2304; k++) {
4059 fCurrentCoefDetector[k] = 0.0;
4063 if(fDebugLevel > 1){
4065 if ( !fDebugStreamer ) {
4067 TDirectory *backup = gDirectory;
4068 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
4069 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4072 Int_t detector = fCountDet;
4073 Int_t caligroup = idect;
4074 Short_t rowmin = fCalibraMode->GetRowMin(0);
4075 Short_t rowmax = fCalibraMode->GetRowMax(0);
4076 Short_t colmin = fCalibraMode->GetColMin(0);
4077 Short_t colmax = fCalibraMode->GetColMax(0);
4078 Float_t gf = fCurrentCoef[0];
4079 Float_t gfs = fCurrentCoef[1];
4080 Float_t gfE = fCurrentCoefE;
4082 (*fDebugStreamer) << "FillFillCH" <<
4083 "detector=" << detector <<
4084 "caligroup=" << caligroup <<
4085 "rowmin=" << rowmin <<
4086 "rowmax=" << rowmax <<
4087 "colmin=" << colmin <<
4088 "colmax=" << colmax <<
4096 //________________________________________________________________________________
4097 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
4100 // DebugStream and fVectorFit and fVectorFit2
4103 // End of one detector
4104 if (idect == (fCount-1)) {
4108 for (Int_t k = 0; k < 2304; k++) {
4109 fCurrentCoefDetector[k] = 0.0;
4110 fCurrentCoefDetector2[k] = 0.0;
4114 if(fDebugLevel > 1){
4116 if ( !fDebugStreamer ) {
4118 TDirectory *backup = gDirectory;
4119 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
4120 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4124 Int_t detector = fCountDet;
4125 Int_t caligroup = idect;
4126 Short_t rowmin = fCalibraMode->GetRowMin(1);
4127 Short_t rowmax = fCalibraMode->GetRowMax(1);
4128 Short_t colmin = fCalibraMode->GetColMin(1);
4129 Short_t colmax = fCalibraMode->GetColMax(1);
4130 Float_t vf = fCurrentCoef[0];
4131 Float_t vs = fCurrentCoef[1];
4132 Float_t vfE = fCurrentCoefE;
4133 Float_t t0f = fCurrentCoef2[0];
4134 Float_t t0s = fCurrentCoef2[1];
4135 Float_t t0E = fCurrentCoefE2;
4139 (* fDebugStreamer) << "FillFillPH"<<
4140 "detector="<<detector<<
4141 "nentries="<<nentries<<
4142 "caligroup="<<caligroup<<
4157 //________________________________________________________________________________
4158 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
4161 // DebugStream and fVectorFit
4164 // End of one detector
4165 if (idect == (fCount-1)) {
4168 for (Int_t k = 0; k < 2304; k++) {
4169 fCurrentCoefDetector[k] = 0.0;
4174 if(fDebugLevel > 1){
4176 if ( !fDebugStreamer ) {
4178 TDirectory *backup = gDirectory;
4179 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4180 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4183 Int_t detector = fCountDet;
4184 Int_t layer = GetLayer(fCountDet);
4185 Int_t caligroup = idect;
4186 Short_t rowmin = fCalibraMode->GetRowMin(2);
4187 Short_t rowmax = fCalibraMode->GetRowMax(2);
4188 Short_t colmin = fCalibraMode->GetColMin(2);
4189 Short_t colmax = fCalibraMode->GetColMax(2);
4190 Float_t widf = fCurrentCoef[0];
4191 Float_t wids = fCurrentCoef[1];
4192 Float_t widfE = fCurrentCoefE;
4194 (* fDebugStreamer) << "FillFillPRF"<<
4195 "detector="<<detector<<
4197 "caligroup="<<caligroup<<
4209 //________________________________________________________________________________
4210 void AliTRDCalibraFit::FillFillLinearFitter()
4213 // DebugStream and fVectorFit
4216 // End of one detector
4222 for (Int_t k = 0; k < 2304; k++) {
4223 fCurrentCoefDetector[k] = 0.0;
4224 fCurrentCoefDetector2[k] = 0.0;
4228 if(fDebugLevel > 1){
4230 if ( !fDebugStreamer ) {
4232 TDirectory *backup = gDirectory;
4233 fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
4234 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4237 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4238 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4239 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4240 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4241 Float_t tiltangle = padplane->GetTiltingAngle();
4242 Int_t detector = fCountDet;
4243 Int_t stack = GetStack(fCountDet);
4244 Int_t layer = GetLayer(fCountDet);
4245 Float_t vf = fCurrentCoef[0];
4246 Float_t vs = fCurrentCoef[1];
4247 Float_t vfE = fCurrentCoefE;
4248 Float_t lorentzangler = fCurrentCoef2[0];
4249 Float_t elorentzangler = fCurrentCoefE2;
4250 Float_t lorentzangles = fCurrentCoef2[1];
4252 (* fDebugStreamer) << "FillFillLinearFitter"<<
4253 "detector="<<detector<<
4258 "tiltangle="<<tiltangle<<
4262 "lorentzangler="<<lorentzangler<<
4263 "Elorentzangler="<<elorentzangler<<
4264 "lorentzangles="<<lorentzangles<<
4269 //________________________________________________________________________________
4270 void AliTRDCalibraFit::FillFillExbAlt()
4273 // DebugStream and fVectorFit
4276 // End of one detector
4281 for (Int_t k = 0; k < 2304; k++) {
4282 fCurrentCoefDetector2[k] = 0.0;
4286 if(fDebugLevel > 1){
4288 if ( !fDebugStreamer ) {
4290 TDirectory *backup = gDirectory;
4291 fDebugStreamer = new TTreeSRedirector("TRDDebugFitExbAlt.root");
4292 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
4295 //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
4296 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
4297 Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
4298 Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
4299 Float_t tiltangle = padplane->GetTiltingAngle();
4300 Int_t detector = fCountDet;
4301 Int_t stack = GetStack(fCountDet);
4302 Int_t layer = GetLayer(fCountDet);
4303 Float_t vf = fCurrentCoef2[0];
4304 Float_t vfE = fCurrentCoefE2;
4306 (* fDebugStreamer) << "FillFillLinearFitter"<<
4307 "detector="<<detector<<
4312 "tiltangle="<<tiltangle<<
4320 //____________Calcul Coef Mean_________________________________________________
4322 //_____________________________________________________________________________
4323 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
4326 // For the detector Dect calcul the mean time 0
4327 // for the calibration group idect from the choosen database
4330 fCurrentCoef2[1] = 0.0;
4331 if(fDebugLevel != 1){
4332 if(((fCalibraMode->GetNz(1) > 0) ||
4333 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4335 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4336 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4337 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4341 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4347 fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4351 for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
4352 for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
4353 fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
4356 fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
4364 //_____________________________________________________________________________
4365 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
4368 // For the detector Dect calcul the mean gain factor
4369 // for the calibration group idect from the choosen database
4372 fCurrentCoef[1] = 0.0;
4373 if(fDebugLevel != 1){
4374 if (((fCalibraMode->GetNz(0) > 0) ||
4375 (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) {
4376 for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
4377 for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
4378 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4379 if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4382 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
4386 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4387 if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
4392 //_____________________________________________________________________________
4393 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
4396 // For the detector Dect calcul the mean sigma of pad response
4397 // function for the calibration group idect from the choosen database
4400 fCurrentCoef[1] = 0.0;
4401 if(fDebugLevel != 1){
4402 for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
4403 for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
4404 fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
4407 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
4411 //_____________________________________________________________________________
4412 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
4415 // For the detector dect calcul the mean drift velocity for the
4416 // calibration group idect from the choosen database
4419 fCurrentCoef[1] = 0.0;
4420 if(fDebugLevel != 1){
4421 if (((fCalibraMode->GetNz(1) > 0) ||
4422 (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) {
4424 for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
4425 for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
4426 fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
4430 fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4435 fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4440 //_____________________________________________________________________________
4441 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4444 // For the detector fCountDet, mean drift velocity and tan lorentzangle
4447 fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet);
4448 fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet);
4452 //_____________________________________________________________________________
4453 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4456 // Default width of the PRF if there is no database as reference
4461 //case 0: return 0.515;
4462 //case 1: return 0.502;
4463 //case 2: return 0.491;
4464 //case 3: return 0.481;
4465 //case 4: return 0.471;
4466 //case 5: return 0.463;
4467 //default: return 0.0;
4470 case 0: return 0.538429;
4471 case 1: return 0.524302;
4472 case 2: return 0.511591;
4473 case 3: return 0.500140;
4474 case 4: return 0.489821;
4475 case 5: return 0.480524;
4476 default: return 0.0;
4479 //________________________________________________________________________________
4480 void AliTRDCalibraFit::SetCalROC(Int_t i)
4483 // Set the calib object for fCountDet
4486 Float_t value = 0.0;
4488 //Get the CalDet object
4490 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
4492 AliInfo("Could not get calibDB");
4499 fCalROC->~AliTRDCalROC();
4500 new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4501 }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4505 fCalROC->~AliTRDCalROC();
4506 new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4507 }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4509 fCalROC2->~AliTRDCalROC();
4510 new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4511 }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4515 fCalROC->~AliTRDCalROC();
4516 new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4517 }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4526 if(fCalROC) delete fCalROC;
4527 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4528 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4529 fCalROC->SetValue(k,1.0);
4533 if(fCalROC) delete fCalROC;
4534 if(fCalROC2) delete fCalROC2;
4535 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4536 fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4537 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4538 fCalROC->SetValue(k,1.0);
4539 fCalROC2->SetValue(k,0.0);
4543 if(fCalROC) delete fCalROC;
4544 value = GetPRFDefault(GetLayer(fCountDet));
4545 fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4546 for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4547 fCalROC->SetValue(k,value);
4555 //____________Fit Methods______________________________________________________
4557 //_____________________________________________________________________________
4558 void AliTRDCalibraFit::FitPente(TH1* projPH)
4561 // Slope methode for the drift velocity
4565 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4572 fCurrentCoefE = 0.0;
4573 fCurrentCoefE2 = 0.0;
4574 fCurrentCoef[0] = 0.0;
4575 fCurrentCoef2[0] = 0.0;
4576 TLine *line = new TLine();
4579 TAxis *xpph = projPH->GetXaxis();
4580 Int_t nbins = xpph->GetNbins();
4581 Double_t lowedge = xpph->GetBinLowEdge(1);
4582 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4583 Double_t widbins = (upedge - lowedge) / nbins;
4584 Double_t limit = upedge + 0.5 * widbins;
4587 // Beginning of the signal
4588 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4589 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4590 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4592 binmax = (Int_t) pentea->GetMaximumBin();
4595 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4597 if (binmax >= nbins) {
4600 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4602 pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4603 Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4604 Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4605 Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4606 Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4607 if (TMath::Abs(l3P2am) > 0.00000001) {
4608 fPhd[0] = -(l3P1am / (2 * l3P2am));
4611 if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4612 fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4615 // Amplification region
4618 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4619 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4626 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4628 if (binmax >= nbins) {
4631 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4633 projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4634 Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4635 Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4636 Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4637 Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4638 if (TMath::Abs(l3P2amf) > 0.00000000001) {
4639 fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4641 if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4642 fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4645 fCurrentCoefE2 = fCurrentCoefE;
4648 TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4649 for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
4650 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4653 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4656 AliInfo("Put the binmax from 1 to 2 to enable the fit");
4658 if (binmin >= nbins) {
4661 AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4666 ,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
4667 ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4668 Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4669 Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4670 Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4671 Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4672 if (TMath::Abs(l3P2dr) > 0.00000001) {
4673 fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4675 if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4676 fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4678 Float_t fPhdt0 = 0.0;
4679 Float_t t0Shift = 0.0;
4682 t0Shift = fT0Shift1;
4686 t0Shift = fT0Shift0;
4689 if ((fPhd[2] > fPhd[0]) &&
4690 (fPhd[2] > fPhd[1]) &&
4691 (fPhd[1] > fPhd[0]) &&
4693 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4694 fNumberFitSuccess++;
4696 if (fPhdt0 >= 0.0) {
4697 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4698 if (fCurrentCoef2[0] < -3.0) {
4699 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4703 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4708 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4709 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4712 if (fDebugLevel == 1) {
4713 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4716 line->SetLineColor(2);
4717 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4718 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4719 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4720 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
4721 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
4722 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
4723 AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4724 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4727 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4736 //_____________________________________________________________________________
4737 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4740 // Slope methode but with polynomes de Lagrange
4744 const Float_t kDrWidth = AliTRDgeometry::DrThick();
4747 //Double_t *x = new Double_t[5];
4748 //Double_t *y = new Double_t[5];
4765 fCurrentCoefE = 0.0;
4766 fCurrentCoefE2 = 1.0;
4767 fCurrentCoef[0] = 0.0;
4768 fCurrentCoef2[0] = 0.0;
4769 TLine *line = new TLine();
4770 TF1 * polynome = 0x0;
4771 TF1 * polynomea = 0x0;
4772 TF1 * polynomeb = 0x0;
4780 TAxis *xpph = projPH->GetXaxis();
4781 Int_t nbins = xpph->GetNbins();
4782 Double_t lowedge = xpph->GetBinLowEdge(1);
4783 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4784 Double_t widbins = (upedge - lowedge) / nbins;
4785 Double_t limit = upedge + 0.5 * widbins;
4790 // Beginning of the signal
4791 TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4792 for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
4793 pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4796 binmax = (Int_t) pentea->GetMaximumBin();
4798 Double_t minnn = 0.0;
4799 Double_t maxxx = 0.0;
4801 Int_t kase = nbins-binmax;
4809 minnn = pentea->GetBinCenter(binmax-2);
4810 maxxx = pentea->GetBinCenter(binmax);
4811 x[0] = pentea->GetBinCenter(binmax-2);
4812 x[1] = pentea->GetBinCenter(binmax-1);
4813 x[2] = pentea->GetBinCenter(binmax);
4814 y[0] = pentea->GetBinContent(binmax-2);
4815 y[1] = pentea->GetBinContent(binmax-1);
4816 y[2] = pentea->GetBinContent(binmax);
4817 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4818 //AliInfo("At the limit for beginning!");
4821 minnn = pentea->GetBinCenter(binmax-2);
4822 maxxx = pentea->GetBinCenter(binmax+1);
4823 x[0] = pentea->GetBinCenter(binmax-2);
4824 x[1] = pentea->GetBinCenter(binmax-1);
4825 x[2] = pentea->GetBinCenter(binmax);
4826 x[3] = pentea->GetBinCenter(binmax+1);
4827 y[0] = pentea->GetBinContent(binmax-2);
4828 y[1] = pentea->GetBinContent(binmax-1);
4829 y[2] = pentea->GetBinContent(binmax);
4830 y[3] = pentea->GetBinContent(binmax+1);
4831 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4839 minnn = pentea->GetBinCenter(binmax);
4840 maxxx = pentea->GetBinCenter(binmax+2);
4841 x[0] = pentea->GetBinCenter(binmax);
4842 x[1] = pentea->GetBinCenter(binmax+1);
4843 x[2] = pentea->GetBinCenter(binmax+2);
4844 y[0] = pentea->GetBinContent(binmax);
4845 y[1] = pentea->GetBinContent(binmax+1);
4846 y[2] = pentea->GetBinContent(binmax+2);
4847 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4850 minnn = pentea->GetBinCenter(binmax-1);
4851 maxxx = pentea->GetBinCenter(binmax+2);
4852 x[0] = pentea->GetBinCenter(binmax-1);
4853 x[1] = pentea->GetBinCenter(binmax);
4854 x[2] = pentea->GetBinCenter(binmax+1);
4855 x[3] = pentea->GetBinCenter(binmax+2);
4856 y[0] = pentea->GetBinContent(binmax-1);
4857 y[1] = pentea->GetBinContent(binmax);
4858 y[2] = pentea->GetBinContent(binmax+1);
4859 y[3] = pentea->GetBinContent(binmax+2);
4860 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4863 minnn = pentea->GetBinCenter(binmax-2);
4864 maxxx = pentea->GetBinCenter(binmax+2);
4865 x[0] = pentea->GetBinCenter(binmax-2);
4866 x[1] = pentea->GetBinCenter(binmax-1);
4867 x[2] = pentea->GetBinCenter(binmax);
4868 x[3] = pentea->GetBinCenter(binmax+1);
4869 x[4] = pentea->GetBinCenter(binmax+2);
4870 y[0] = pentea->GetBinContent(binmax-2);
4871 y[1] = pentea->GetBinContent(binmax-1);
4872 y[2] = pentea->GetBinContent(binmax);
4873 y[3] = pentea->GetBinContent(binmax+1);
4874 y[4] = pentea->GetBinContent(binmax+2);
4875 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4883 polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4884 polynomeb->SetParameters(c0,c1,c2,c3,c4);
4886 Double_t step = (maxxx-minnn)/10000;
4888 Double_t maxvalue = 0.0;
4889 Double_t placemaximum = minnn;
4890 for(Int_t o = 0; o < 10000; o++){
4891 if(o == 0) maxvalue = polynomeb->Eval(l);
4892 if(maxvalue < (polynomeb->Eval(l))){
4893 maxvalue = polynomeb->Eval(l);
4898 fPhd[0] = placemaximum;
4901 // Amplification region
4904 for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4905 if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4911 Double_t minn = 0.0;
4912 Double_t maxx = 0.0;
4924 Int_t kase1 = nbins - binmax;
4926 //Determination of minn and maxx
4927 //case binmax = nbins
4932 minn = projPH->GetBinCenter(binmax-2);
4933 maxx = projPH->GetBinCenter(binmax);
4934 x[0] = projPH->GetBinCenter(binmax-2);
4935 x[1] = projPH->GetBinCenter(binmax-1);
4936 x[2] = projPH->GetBinCenter(binmax);
4937 y[0] = projPH->GetBinContent(binmax-2);
4938 y[1] = projPH->GetBinContent(binmax-1);
4939 y[2] = projPH->GetBinContent(binmax);
4940 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4941 //AliInfo("At the limit for the drift!");
4944 minn = projPH->GetBinCenter(binmax-2);
4945 maxx = projPH->GetBinCenter(binmax+1);
4946 x[0] = projPH->GetBinCenter(binmax-2);
4947 x[1] = projPH->GetBinCenter(binmax-1);
4948 x[2] = projPH->GetBinCenter(binmax);
4949 x[3] = projPH->GetBinCenter(binmax+1);
4950 y[0] = projPH->GetBinContent(binmax-2);
4951 y[1] = projPH->GetBinContent(binmax-1);
4952 y[2] = projPH->GetBinContent(binmax);
4953 y[3] = projPH->GetBinContent(binmax+1);
4954 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4963 minn = projPH->GetBinCenter(binmax);
4964 maxx = projPH->GetBinCenter(binmax+2);
4965 x[0] = projPH->GetBinCenter(binmax);
4966 x[1] = projPH->GetBinCenter(binmax+1);
4967 x[2] = projPH->GetBinCenter(binmax+2);
4968 y[0] = projPH->GetBinContent(binmax);
4969 y[1] = projPH->GetBinContent(binmax+1);
4970 y[2] = projPH->GetBinContent(binmax+2);
4971 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4974 minn = projPH->GetBinCenter(binmax-1);
4975 maxx = projPH->GetBinCenter(binmax+2);
4976 x[0] = projPH->GetBinCenter(binmax-1);
4977 x[1] = projPH->GetBinCenter(binmax);
4978 x[2] = projPH->GetBinCenter(binmax+1);
4979 x[3] = projPH->GetBinCenter(binmax+2);
4980 y[0] = projPH->GetBinContent(binmax-1);
4981 y[1] = projPH->GetBinContent(binmax);
4982 y[2] = projPH->GetBinContent(binmax+1);
4983 y[3] = projPH->GetBinContent(binmax+2);
4984 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4987 minn = projPH->GetBinCenter(binmax-2);
4988 maxx = projPH->GetBinCenter(binmax+2);
4989 x[0] = projPH->GetBinCenter(binmax-2);
4990 x[1] = projPH->GetBinCenter(binmax-1);
4991 x[2] = projPH->GetBinCenter(binmax);
4992 x[3] = projPH->GetBinCenter(binmax+1);
4993 x[4] = projPH->GetBinCenter(binmax+2);
4994 y[0] = projPH->GetBinContent(binmax-2);
4995 y[1] = projPH->GetBinContent(binmax-1);
4996 y[2] = projPH->GetBinContent(binmax);
4997 y[3] = projPH->GetBinContent(binmax+1);
4998 y[4] = projPH->GetBinContent(binmax+2);
4999 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5006 polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
5007 polynomea->SetParameters(c0,c1,c2,c3,c4);
5009 Double_t step = (maxx-minn)/1000;
5011 Double_t maxvalue = 0.0;
5012 Double_t placemaximum = minn;
5013 for(Int_t o = 0; o < 1000; o++){
5014 if(o == 0) maxvalue = polynomea->Eval(l);
5015 if(maxvalue < (polynomea->Eval(l))){
5016 maxvalue = polynomea->Eval(l);
5021 fPhd[1] = placemaximum;
5025 TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
5026 for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
5027 pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
5030 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5036 //AliInfo("Put the binmax from 1 to 2 to enable the fit");
5040 if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
5041 //AliInfo("Too many fluctuations at the end!");
5044 if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
5045 //AliInfo("Too many fluctuations at the end!");
5048 if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
5049 //AliInfo("No entries for the next bin!");
5050 pente->SetBinContent(binmin,0);
5051 if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
5067 Bool_t case1 = kFALSE;
5068 Bool_t case2 = kFALSE;
5069 Bool_t case4 = kFALSE;
5071 //Determination of min and max
5072 //case binmin <= nbins-3
5074 if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5075 min = pente->GetBinCenter(binmin-2);
5076 max = pente->GetBinCenter(binmin+2);
5077 x[0] = pente->GetBinCenter(binmin-2);
5078 x[1] = pente->GetBinCenter(binmin-1);
5079 x[2] = pente->GetBinCenter(binmin);
5080 x[3] = pente->GetBinCenter(binmin+1);
5081 x[4] = pente->GetBinCenter(binmin+2);
5082 y[0] = pente->GetBinContent(binmin-2);
5083 y[1] = pente->GetBinContent(binmin-1);
5084 y[2] = pente->GetBinContent(binmin);
5085 y[3] = pente->GetBinContent(binmin+1);
5086 y[4] = pente->GetBinContent(binmin+2);
5087 //Calcul the polynome de Lagrange
5088 CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
5090 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5091 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5092 //AliInfo("polynome 4 false 1");
5095 if(((binmin+3) <= (nbins-1)) &&
5096 (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
5097 ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
5098 (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
5099 //AliInfo("polynome 4 false 2");
5103 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
5104 (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
5105 //AliInfo("polynome 4 case 1");
5108 if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
5109 (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5110 //AliInfo("polynome 4 case 4");
5115 //case binmin = nbins-2
5117 if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5119 min = pente->GetBinCenter(binmin-2);
5120 max = pente->GetBinCenter(binmin+1);
5121 x[0] = pente->GetBinCenter(binmin-2);
5122 x[1] = pente->GetBinCenter(binmin-1);
5123 x[2] = pente->GetBinCenter(binmin);
5124 x[3] = pente->GetBinCenter(binmin+1);
5125 y[0] = pente->GetBinContent(binmin-2);
5126 y[1] = pente->GetBinContent(binmin-1);
5127 y[2] = pente->GetBinContent(binmin);
5128 y[3] = pente->GetBinContent(binmin+1);
5129 //Calcul the polynome de Lagrange
5130 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5131 //richtung +: nothing
5133 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5134 //AliInfo("polynome 3- case 2");
5139 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5141 min = pente->GetBinCenter(binmin-1);
5142 max = pente->GetBinCenter(binmin+2);
5143 x[0] = pente->GetBinCenter(binmin-1);
5144 x[1] = pente->GetBinCenter(binmin);
5145 x[2] = pente->GetBinCenter(binmin+1);
5146 x[3] = pente->GetBinCenter(binmin+2);
5147 y[0] = pente->GetBinContent(binmin-1);
5148 y[1] = pente->GetBinContent(binmin);
5149 y[2] = pente->GetBinContent(binmin+1);
5150 y[3] = pente->GetBinContent(binmin+2);
5151 //Calcul the polynome de Lagrange
5152 CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
5154 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5155 //AliInfo("polynome 3+ case 2");
5160 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
5161 min = pente->GetBinCenter(binmin);
5162 max = pente->GetBinCenter(binmin+2);
5163 x[0] = pente->GetBinCenter(binmin);
5164 x[1] = pente->GetBinCenter(binmin+1);
5165 x[2] = pente->GetBinCenter(binmin+2);
5166 y[0] = pente->GetBinContent(binmin);
5167 y[1] = pente->GetBinContent(binmin+1);
5168 y[2] = pente->GetBinContent(binmin+2);
5169 //Calcul the polynome de Lagrange
5170 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5172 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
5173 //AliInfo("polynome 2+ false");
5178 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
5180 min = pente->GetBinCenter(binmin-1);
5181 max = pente->GetBinCenter(binmin+1);
5182 x[0] = pente->GetBinCenter(binmin-1);
5183 x[1] = pente->GetBinCenter(binmin);
5184 x[2] = pente->GetBinCenter(binmin+1);
5185 y[0] = pente->GetBinContent(binmin-1);
5186 y[1] = pente->GetBinContent(binmin);
5187 y[2] = pente->GetBinContent(binmin+1);
5188 //Calcul the polynome de Lagrange
5189 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5190 //richtung +: nothing
5191 //richtung -: nothing
5193 //case binmin = nbins-1
5195 if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
5196 min = pente->GetBinCenter(binmin-2);
5197 max = pente->GetBinCenter(binmin);
5198 x[0] = pente->GetBinCenter(binmin-2);
5199 x[1] = pente->GetBinCenter(binmin-1);
5200 x[2] = pente->GetBinCenter(binmin);
5201 y[0] = pente->GetBinContent(binmin-2);
5202 y[1] = pente->GetBinContent(binmin-1);
5203 y[2] = pente->GetBinContent(binmin);
5204 //Calcul the polynome de Lagrange
5205 CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
5206 //AliInfo("At the limit for the drift!");
5207 //fluctuation too big!
5208 //richtung +: nothing
5210 if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
5211 //AliInfo("polynome 2- false ");
5215 if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
5217 //AliInfo("At the limit for the drift and not usable!");
5221 if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
5223 //AliInfo("For the drift...problem!");
5225 //pass but should not happen
5226 if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
5228 //AliInfo("For the drift...problem!");
5232 polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
5233 polynome->SetParameters(c0,c1,c2,c3,c4);
5234 //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
5235 Double_t step = (max-min)/1000;
5237 Double_t minvalue = 0.0;
5238 Double_t placeminimum = min;
5239 for(Int_t o = 0; o < 1000; o++){
5240 if(o == 0) minvalue = polynome->Eval(l);
5241 if(minvalue > (polynome->Eval(l))){
5242 minvalue = polynome->Eval(l);
5247 fPhd[2] = placeminimum;
5249 //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
5250 if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
5251 if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
5253 Float_t fPhdt0 = 0.0;
5254 Float_t t0Shift = 0.0;
5257 t0Shift = fT0Shift1;
5261 t0Shift = fT0Shift0;
5264 if ((fPhd[2] > fPhd[0]) &&
5265 (fPhd[2] > fPhd[1]) &&
5266 (fPhd[1] > fPhd[0]) &&
5268 fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
5269 if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5270 else fNumberFitSuccess++;
5271 if (fPhdt0 >= 0.0) {
5272 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5273 //printf("Value of timeoffset %f\n",fCurrentCoef2[0]);
5274 if (fCurrentCoef2[0] < -3.0) {
5275 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5279 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5283 ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
5284 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5286 if((fPhd[1] > fPhd[0]) &&
5288 if (fPhdt0 >= 0.0) {
5289 fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
5290 if (fCurrentCoef2[0] < -3.0) {
5291 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5295 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5299 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5300 //printf("Fit failed!\n");
5304 if (fDebugLevel == 1) {
5305 TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
5308 line->SetLineColor(2);
5309 line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
5310 line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
5311 line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
5312 AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
5313 AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
5314 AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
5315 AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
5316 TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
5319 TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
5326 if(polynome) delete polynome;
5327 if(polynomea) delete polynomea;
5328 if(polynomeb) delete polynomeb;
5329 //if(x) delete [] x;
5330 //if(y) delete [] y;
5331 if(line) delete line;
5336 //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5337 //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5338 //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]);
5339 projPH->SetDirectory(0);
5343 //_____________________________________________________________________________
5344 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
5347 // Fit methode for the drift velocity
5351 const Float_t kDrWidth = AliTRDgeometry::DrThick();
5354 TAxis *xpph = projPH->GetXaxis();
5355 Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
5357 TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
5358 fPH->SetParameter(0,0.469); // Scaling
5359 fPH->SetParameter(1,0.18); // Start
5360 fPH->SetParameter(2,0.0857325); // AR
5361 fPH->SetParameter(3,1.89); // DR
5362 fPH->SetParameter(4,0.08); // QA/QD
5363 fPH->SetParameter(5,0.0); // Baseline
5365 TLine *line = new TLine();
5367 fCurrentCoef[0] = 0.0;
5368 fCurrentCoef2[0] = 0.0;
5369 fCurrentCoefE = 0.0;
5370 fCurrentCoefE2 = 0.0;
5372 if (idect%fFitPHPeriode == 0) {
5374 AliInfo(Form("The detector %d will be fitted",idect));
5375 fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
5376 fPH->SetParameter(1,fPhd[0] - 0.1); // Start
5377 fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
5378 fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
5379 fPH->SetParameter(4,0.225); // QA/QD
5380 fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
5382 if (fDebugLevel != 1) {
5383 projPH->Fit(fPH,"0M","",0.0,upedge);
5386 TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
5388 projPH->Fit(fPH,"M+","",0.0,upedge);
5390 line->SetLineColor(4);
5391 line->DrawLine(fPH->GetParameter(1)
5393 ,fPH->GetParameter(1)
5394 ,projPH->GetMaximum());
5395 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
5397 ,fPH->GetParameter(1)+fPH->GetParameter(2)
5398 ,projPH->GetMaximum());
5399 line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5401 ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
5402 ,projPH->GetMaximum());
5405 if (fPH->GetParameter(3) != 0) {
5406 fNumberFitSuccess++;
5407 fCurrentCoef[0] = kDrWidth / (fPH->GetParameter(3));
5408 fCurrentCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
5409 fCurrentCoef2[0] = fPH->GetParameter(1);
5410 fCurrentCoefE2 = fPH->GetParError(1);
5413 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5414 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5420 // Put the default value
5421 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5422 fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
5425 if (fDebugLevel != 1) {
5430 //_____________________________________________________________________________
5431 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5434 // Fit methode for the sigma of the pad response function
5439 fCurrentCoef[0] = 0.0;
5440 fCurrentCoefE = 0.0;
5442 Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
5444 if(TMath::Abs(ret+4) <= 0.000000001){
5445 fCurrentCoef[0] = -fCurrentCoef[1];
5449 fNumberFitSuccess++;
5450 fCurrentCoef[0] = param[2];
5451 fCurrentCoefE = ret;
5455 //_____________________________________________________________________________
5456 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)
5459 // Fit methode for the sigma of the pad response function
5462 //We should have at least 3 points
5463 if(nBins <=3) return -4.0;
5465 TLinearFitter fitter(3,"pol2");
5466 fitter.StoreData(kFALSE);
5467 fitter.ClearPoints();
5469 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5470 Float_t entries = 0;
5471 Int_t nbbinwithentries = 0;
5472 for (Int_t i=0; i<nBins; i++){
5474 if(arraye[i] > 15) nbbinwithentries++;
5475 //printf("entries for i %d: %f\n",i,arraye[i]);
5477 if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5478 //printf("entries %f\n",entries);
5479 //printf("nbbinwithentries %d\n",nbbinwithentries);
5482 Float_t errorm = 0.0;
5483 Float_t errorn = 0.0;
5484 Float_t error = 0.0;
5487 for (Int_t ibin=0;ibin<nBins; ibin++){
5488 Float_t entriesI = arraye[ibin];
5489 Float_t valueI = arraym[ibin];
5490 Double_t xcenter = 0.0;
5492 if ((entriesI>15) && (valueI>0.0)){
5493 xcenter = xMin+(ibin+0.5)*binWidth;
5498 if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5499 if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5500 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5503 if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5504 if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5505 error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5507 if(TMath::Abs(error) < 0.000000001) continue;
5508 val = TMath::Log(Float_t(valueI));
5509 fitter.AddPoint(&xcenter,val,error);
5513 if(fDebugLevel > 1){
5515 if ( !fDebugStreamer ) {
5517 TDirectory *backup = gDirectory;
5518 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5519 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5522 Int_t detector = fCountDet;
5523 Int_t layer = GetLayer(fCountDet);
5526 (* fDebugStreamer) << "FitGausMIFill"<<
5527 "detector="<<detector<<
5531 "entriesI="<<entriesI<<
5534 "xcenter="<<xcenter<<
5544 if(npoints <=3) return -4.0;
5549 fitter.GetParameters(par);
5550 chi2 = fitter.GetChisquare()/Float_t(npoints);
5553 if (!param) param = new TVectorD(3);
5554 if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5555 Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
5556 Double_t deltax = (fitter.GetParError(2))/x;
5557 Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5560 (*param)[1] = par[1]/(-2.*par[2]);
5561 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5562 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
5563 if ( lnparam0>307 ) return -4;
5564 (*param)[0] = TMath::Exp(lnparam0);
5566 if(fDebugLevel > 1){
5568 if ( !fDebugStreamer ) {
5570 TDirectory *backup = gDirectory;
5571 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5572 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5575 Int_t detector = fCountDet;
5576 Int_t layer = GetLayer(fCountDet);
5579 (* fDebugStreamer) << "FitGausMIFit"<<
5580 "detector="<<detector<<
5583 "errorsigma="<<chi2<<
5584 "mean="<<(*param)[1]<<
5585 "sigma="<<(*param)[2]<<
5586 "constant="<<(*param)[0]<<
5591 if((chi2/(*param)[2]) > 0.1){
5593 chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5598 if(fDebugLevel == 1){
5599 TString name("PRF");
5600 name += (Int_t)xMin;
5601 name += (Int_t)xMax;
5602 TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5605 TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5606 for(Int_t k = 0; k < nBins; k++){
5607 histo->SetBinContent(k+1,arraym[k]);
5608 histo->SetBinError(k+1,arrayme[k]);
5611 name += "functionf";
5612 TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5613 f1->SetParameter(0, (*param)[0]);
5614 f1->SetParameter(1, (*param)[1]);
5615 f1->SetParameter(2, (*param)[2]);
5623 //_____________________________________________________________________________
5624 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5627 // Fit methode for the sigma of the pad response function
5630 fCurrentCoef[0] = 0.0;
5631 fCurrentCoefE = 0.0;
5633 if (fDebugLevel != 1) {
5634 projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5637 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5639 projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5642 fCurrentCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
5643 fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5644 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5646 fNumberFitSuccess++;
5649 //_____________________________________________________________________________
5650 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5653 // Fit methode for the sigma of the pad response function
5655 fCurrentCoef[0] = 0.0;
5656 fCurrentCoefE = 0.0;
5657 if (fDebugLevel == 1) {
5658 TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5662 fCurrentCoef[0] = projPRF->GetRMS();
5663 if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5665 fNumberFitSuccess++;
5668 //_____________________________________________________________________________
5669 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5672 // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5675 TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5678 Int_t nbins = (Int_t)(nybins/(2*nbg));
5679 Float_t lowedge = -3.0*nbg;
5680 Float_t upedge = lowedge + 3.0;
5683 Double_t xvalues = -0.2*nbg+0.1;
5685 Int_t total = 2*nbg;
5688 for(Int_t k = 0; k < total; k++){
5689 if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5691 y = fCurrentCoef[0]*fCurrentCoef[0];
5692 linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5695 if(fDebugLevel > 1){
5697 if ( !fDebugStreamer ) {
5699 TDirectory *backup = gDirectory;
5700 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5701 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5704 Int_t detector = fCountDet;
5705 Int_t layer = GetLayer(fCountDet);
5706 Int_t nbtotal = total;
5708 Float_t low = lowedge;
5709 Float_t up = upedge;
5710 Float_t tnp = xvalues;
5711 Float_t wid = fCurrentCoef[0];
5712 Float_t widfE = fCurrentCoefE;
5714 (* fDebugStreamer) << "FitTnpRange0"<<
5715 "detector="<<detector<<
5717 "nbtotal="<<nbtotal<<
5735 fCurrentCoefE = 0.0;
5736 fCurrentCoef[0] = 0.0;
5738 //printf("npoints\n",npoints);
5741 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5746 linearfitter.Eval();
5747 linearfitter.GetParameters(pars0);
5748 Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5749 Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
5750 Double_t min0 = 0.0;
5751 Double_t ermin0 = 0.0;
5752 //Double_t prfe0 = 0.0;
5753 Double_t prf0 = 0.0;
5754 if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5755 min0 = -pars0[1]/(2*pars0[2]);
5756 ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5757 prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5760 prfe0 = linearfitter->GetParError(0)*pointError0
5761 +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5762 +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5763 prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5764 fCurrentCoefE = (Float_t) prfe0;
5766 fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5769 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5773 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5776 if(fDebugLevel > 1){
5778 if ( !fDebugStreamer ) {
5780 TDirectory *backup = gDirectory;
5781 fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5782 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
5785 Int_t detector = fCountDet;
5786 Int_t layer = GetLayer(fCountDet);
5787 Int_t nbtotal = total;
5788 Double_t colsize[6] = {0.635,0.665,0.695,0.725,0.755,0.785};
5789 Double_t sigmax = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5791 (* fDebugStreamer) << "FitTnpRange1"<<
5792 "detector="<<detector<<
5794 "nbtotal="<<nbtotal<<
5798 "npoints="<<npoints<<
5801 "sigmaprf="<<fCurrentCoef[0]<<
5802 "sigprf="<<fCurrentCoef[1]<<
5809 //_____________________________________________________________________________
5810 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5813 // Only mean methode for the gain factor
5816 fCurrentCoef[0] = mean;
5817 fCurrentCoefE = 0.0;
5818 if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5819 if (fDebugLevel == 1) {
5820 TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5824 CalculChargeCoefMean(kTRUE);
5825 fNumberFitSuccess++;
5827 //_____________________________________________________________________________
5828 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5831 // mean w methode for the gain factor
5835 Int_t nybins = projch->GetNbinsX();
5837 //The weight function
5838 Double_t a = 0.00228515;
5839 Double_t b = -0.00231487;
5840 Double_t c = 0.00044298;
5841 Double_t d = -0.00379239;
5842 Double_t e = 0.00338349;
5852 //A arbitrary error for the moment
5853 fCurrentCoefE = 0.0;
5854 fCurrentCoef[0] = 0.0;
5857 Double_t sumw = 0.0;
5859 Float_t sumAll = (Float_t) nentries;
5860 Int_t sumCurrent = 0;
5861 for(Int_t k = 0; k <nybins; k++){
5862 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5863 if (fraction<fOutliersFitChargeLow) {
5864 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5865 //printf("Take only after bin %d\n",k);
5868 if (fraction>fOutliersFitChargeHigh) {
5869 //printf("Break by the bin %d\n",k);
5872 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5873 e*fraction*fraction*fraction*fraction;
5874 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5875 sum += weight*projch->GetBinContent(k+1);
5876 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5877 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5879 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5881 if (fDebugLevel == 1) {
5882 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5885 TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
5888 fNumberFitSuccess++;
5889 CalculChargeCoefMean(kTRUE);
5891 //_____________________________________________________________________________
5892 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5895 // mean w methode for the gain factor
5899 Int_t nybins = projch->GetNbinsX();
5901 //The weight function
5902 Double_t a = 0.00228515;
5903 Double_t b = -0.00231487;
5904 Double_t c = 0.00044298;
5905 Double_t d = -0.00379239;
5906 Double_t e = 0.00338349;
5916 //A arbitrary error for the moment
5917 fCurrentCoefE = 0.0;
5918 fCurrentCoef[0] = 0.0;
5921 Double_t sumw = 0.0;
5923 Int_t sumCurrent = 0;
5924 for(Int_t k = 0; k <nybins; k++){
5925 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5926 if (fraction>0.95) break;
5927 Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5928 e*fraction*fraction*fraction*fraction;
5929 sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5930 sum += weight*projch->GetBinContent(k+1);
5931 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5932 //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5934 if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5936 if (fDebugLevel == 1) {
5937 TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5940 TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
5943 fNumberFitSuccess++;
5945 //_____________________________________________________________________________
5946 void AliTRDCalibraFit::FitLandau(TH1 *projch, Double_t mean, Double_t nentries)
5949 // Fit methode for the gain factor
5953 //Calcul Range of the fit
5954 Double_t lastvalue = 0.0;
5955 Float_t sumAll = (Float_t) nentries;
5956 Int_t sumCurrent = 0;
5957 //printf("There are %d bins\n",nybins);
5958 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
5959 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5960 if (fraction>fOutliersFitChargeHigh) {
5961 lastvalue = projch->GetBinCenter(k+1);
5962 //printf("Break by %f\n",lastvalue);
5965 sumCurrent += (Int_t) projch->GetBinContent(k+1);
5969 fCurrentCoef[0] = 0.0;
5970 fCurrentCoefE = 0.0;
5971 Double_t chisqrl = 0.0;
5973 projch->Fit("landau","WWQ+",""
5974 ,(Double_t) mean/fBeginFitCharge
5976 chisqrl = projch->GetFunction("landau")->GetChisquare();
5978 if (fDebugLevel == 1) {
5979 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5982 TLine *line = new TLine( projch->GetFunction("landau")->GetParameter(1),0.0,projch->GetFunction("landau")->GetParameter(1),20000.0);
5986 if ((projch->GetFunction("landau")->GetParameter(1) > 0) && (projch->GetFunction("landau")->GetParError(1) < (0.05*projch->GetFunction("landau")->GetParameter(1)))) {
5987 fNumberFitSuccess++;
5988 CalculChargeCoefMean(kTRUE);
5989 fCurrentCoef[0] = projch->GetFunction("landau")->GetParameter(1);
5990 fCurrentCoefE = projch->GetFunction("landau")->GetParError(1);
5993 CalculChargeCoefMean(kFALSE);
5994 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6000 //_____________________________________________________________________________
6001 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean, Double_t nentries)
6004 // Fit methode for the gain factor
6007 //Calcul Range of the fit
6008 Double_t lastvalue = 0.0;
6009 Float_t sumAll = (Float_t) nentries;
6010 Int_t sumCurrent = 0;
6011 //printf("There are %d bins\n",nybins);
6012 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6013 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6014 if (fraction>fOutliersFitChargeHigh) {
6015 lastvalue = projch->GetBinCenter(k+1);
6016 //printf("Break by %f\n",lastvalue);
6019 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6023 fCurrentCoef[0] = 0.0;
6024 fCurrentCoefE = 0.0;
6025 Double_t chisqrl = 0.0;
6026 Double_t chisqrg = 0.0;
6027 Double_t chisqr = 0.0;
6028 TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,(Double_t) mean/fBeginFitCharge,lastvalue,5);
6030 projch->Fit("landau","WWQ0",""
6031 ,(Double_t) mean/fBeginFitCharge
6033 Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
6034 Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
6035 Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
6036 chisqrl = projch->GetFunction("landau")->GetChisquare();
6038 projch->Fit("gaus","WWQ0",""
6039 ,(Double_t) mean/fBeginFitCharge
6041 Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
6042 Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
6043 chisqrg = projch->GetFunction("gaus")->GetChisquare();
6045 fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
6046 if (fDebugLevel != 1) {
6047 projch->Fit("fLandauGaus","WWQ0",""
6048 ,(Double_t) mean/fBeginFitCharge
6050 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6053 TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
6055 projch->Fit("fLandauGaus","WWQ+",""
6056 ,(Double_t) mean/fBeginFitCharge
6058 chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
6060 fLandauGaus->Draw("same");
6061 TLine *line = new TLine(projch->GetFunction("fLandauGaus")->GetParameter(1),0.0,projch->GetFunction("fLandauGaus")->GetParameter(1),20000.0);
6065 if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
6066 fNumberFitSuccess++;
6067 CalculChargeCoefMean(kTRUE);
6068 fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
6069 fCurrentCoefE = projch->GetFunction("fLandauGaus")->GetParError(1);
6072 CalculChargeCoefMean(kFALSE);
6073 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6076 if (fDebugLevel != 1) {
6081 //_____________________________________________________________________________
6082 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean, Double_t nentries)
6085 // Fit methode for the gain factor more time consuming
6088 //Calcul Range of the fit
6089 Double_t lastvalue = 0.0;
6090 Float_t sumAll = (Float_t) nentries;
6091 Int_t sumCurrent = 0;
6092 //printf("There are %d bins\n",nybins);
6093 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6094 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6095 if (fraction>fOutliersFitChargeHigh) {
6096 lastvalue = projch->GetBinCenter(k+1);
6097 //printf("Break by %f\n",lastvalue);
6100 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6104 //Some parameters to initialise
6105 Double_t widthLandau, widthGaus, mPV, integral;
6106 Double_t chisquarel = 0.0;
6107 Double_t chisquareg = 0.0;
6108 projch->Fit("landau","WWQ0M+",""
6109 ,(Double_t) mean/fBeginFitCharge
6111 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6112 chisquarel = projch->GetFunction("landau")->GetChisquare();
6113 projch->Fit("gaus","WWQ0M+",""
6114 ,(Double_t) mean/fBeginFitCharge
6116 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6117 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6119 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6120 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
6122 // Setting fit range and start values
6124 Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
6125 Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
6126 Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
6127 Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
6128 Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
6129 fr[0] = mean/fBeginFitCharge;
6131 fCurrentCoef[0] = 0.0;
6132 fCurrentCoefE = 0.0;
6136 TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
6141 //Double_t projchPeak;
6142 //Double_t projchFWHM;
6143 //LanGauPro(fp,projchPeak,projchFWHM);
6145 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6146 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6147 fNumberFitSuccess++;
6148 CalculChargeCoefMean(kTRUE);
6149 fCurrentCoef[0] = fp[1];
6150 fCurrentCoefE = fpe[1];
6151 //chargeCoefE2 = chisqr;
6154 CalculChargeCoefMean(kFALSE);
6155 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6157 if (fDebugLevel == 1) {
6158 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
6159 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6162 fitsnr->Draw("same");
6163 TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
6170 //_____________________________________________________________________________
6171 void AliTRDCalibraFit::FitBisCHEx(TH1* projch, Double_t mean, Double_t nentries)
6174 // Fit methode for the gain factor more time consuming
6177 //Calcul Range of the fit
6178 Double_t lastvalue = 0.0;
6179 Float_t sumAll = (Float_t) nentries;
6180 Int_t sumCurrent = 0;
6181 //printf("There are %d bins\n",nybins);
6182 for(Int_t k = 0; k <projch->GetNbinsX(); k++){
6183 Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
6184 if (fraction>fOutliersFitChargeHigh) {
6185 lastvalue = projch->GetBinCenter(k+1);
6186 //printf("Break by %f\n",lastvalue);
6189 sumCurrent += (Int_t) projch->GetBinContent(k+1);
6194 //Some parameters to initialise
6195 Double_t widthLandau, widthGaus, mPV, integral;
6196 Double_t chisquarel = 0.0;
6197 Double_t chisquareg = 0.0;
6198 projch->Fit("landau","WWQM+",""
6199 ,(Double_t) mean/fBeginFitCharge
6201 widthLandau = projch->GetFunction("landau")->GetParameter(2);
6202 chisquarel = projch->GetFunction("landau")->GetChisquare();
6203 projch->Fit("gaus","WWQM+",""
6204 ,(Double_t) mean/fBeginFitCharge
6206 widthGaus = projch->GetFunction("gaus")->GetParameter(2);
6207 chisquareg = projch->GetFunction("gaus")->GetChisquare();
6209 mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
6210 integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
6212 // Setting fit range and start values
6214 //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
6215 //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
6216 Double_t sv[5] = { widthLandau, mPV, integral, widthGaus, 0.0};
6217 Double_t pllo[5] = { 0.001, 0.001, projch->Integral()/3, 0.001, 0.0};
6218 Double_t plhi[5] = { 300.0, 300.0, 30*projch->Integral(), 300.0, 2.0};
6219 Double_t fp[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
6220 Double_t fpe[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
6222 //fr[0] = 0.3 * mean;
6223 //fr[1] = 3.0 * mean;
6225 fr[0] = mean/fBeginFitCharge;
6228 fCurrentCoef[0] = 0.0;
6229 fCurrentCoefE = 0.0;
6231 Double_t chisqr = 100.0;
6236 if((mPV > 0.0) && (projch->GetFunction("gaus")->GetParameter(1) > 0.0)) {
6237 fitsnr = LanGauFitEx(projch,&fr[0],&sv[0]
6243 //Double_t projchPeak;
6244 //Double_t projchFWHM;
6245 //LanGauProEx(fp,projchPeak,projchFWHM);
6247 if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
6248 //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
6249 fNumberFitSuccess++;
6250 CalculChargeCoefMean(kTRUE);
6251 fCurrentCoef[0] = fp[1];
6252 fCurrentCoefE = fpe[1];
6253 //chargeCoefE2 = chisqr;
6256 CalculChargeCoefMean(kFALSE);
6257 fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
6259 if (fDebugLevel == 1) {
6260 AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
6261 TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
6264 if(fitsnr) fitsnr->Draw("same");
6265 TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
6272 //_____________________________________________________________________________
6273 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
6276 // Calcul the coefficients of the polynome passant par ces trois points de degre 2
6278 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
6279 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
6280 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
6285 c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
6286 c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
6290 //_____________________________________________________________________________
6291 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
6294 // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
6296 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
6297 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
6298 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
6299 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
6303 c2 = -(x0*(x[1]+x[2]+x[3])
6304 +x1*(x[0]+x[2]+x[3])
6305 +x2*(x[0]+x[1]+x[3])
6306 +x3*(x[0]+x[1]+x[2]));
6307 c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
6308 +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
6309 +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
6310 +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
6312 c0 = -(x0*x[1]*x[2]*x[3]
6315 +x3*x[0]*x[1]*x[2]);
6320 //_____________________________________________________________________________
6321 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
6324 // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
6326 Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
6327 Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
6328 Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
6329 Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
6330 Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
6333 c4 = x0+x1+x2+x3+x4;
6334 c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
6335 +x1*(x[0]+x[2]+x[3]+x[4])
6336 +x2*(x[0]+x[1]+x[3]+x[4])
6337 +x3*(x[0]+x[1]+x[2]+x[4])
6338 +x4*(x[0]+x[1]+x[2]+x[3]));
6339 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])
6340 +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])
6341 +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])
6342 +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])
6343 +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]));
6345 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])
6346 +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])
6347 +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])
6348 +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])
6349 +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]));
6351 c0 = (x0*x[1]*x[2]*x[3]*x[4]
6352 +x1*x[0]*x[2]*x[3]*x[4]
6353 +x2*x[0]*x[1]*x[3]*x[4]
6354 +x3*x[0]*x[1]*x[2]*x[4]
6355 +x4*x[0]*x[1]*x[2]*x[3]);
6358 //_____________________________________________________________________________
6359 void AliTRDCalibraFit::NormierungCharge()
6362 // Normalisation of the gain factor resulting for the fits
6365 // Calcul of the mean of choosen method by fFitChargeNDB
6367 //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
6368 for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
6370 Int_t detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
6371 Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
6372 //printf("detector %d coef[0] %f\n",detector,coef[0]);
6373 if (GetStack(detector) == 2) {
6376 if (GetStack(detector) != 2) {
6379 for (Int_t j = 0; j < total; j++) {
6387 fScaleFitFactor = fScaleFitFactor / sum;
6390 fScaleFitFactor = 1.0;
6393 //methode de boeuf mais bon...
6394 Double_t scalefactor = fScaleFitFactor;
6396 if(fDebugLevel > 1){
6398 if ( !fDebugStreamer ) {
6400 TDirectory *backup = gDirectory;
6401 fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
6402 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
6404 (* fDebugStreamer) << "NormierungCharge"<<
6405 "scalefactor="<<scalefactor<<
6409 //_____________________________________________________________________________
6410 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
6413 // Rebin of the 1D histo for the gain calibration if needed.
6414 // you have to choose fRebin, divider of fNumberBinCharge
6417 TAxis *xhist = hist->GetXaxis();
6418 TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6419 ,xhist->GetBinLowEdge(1)
6420 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6422 AliInfo(Form("fRebin: %d",fRebin));
6424 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6426 for (Int_t ji = i; ji < i+fRebin; ji++) {
6427 sum += hist->GetBinContent(ji);
6430 rehist->SetBinContent(k,sum);
6438 //_____________________________________________________________________________
6439 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
6442 // Rebin of the 1D histo for the gain calibration if needed
6443 // you have to choose fRebin divider of fNumberBinCharge
6446 TAxis *xhist = hist->GetXaxis();
6447 TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
6448 ,xhist->GetBinLowEdge(1)
6449 ,xhist->GetBinUpEdge(xhist->GetNbins()));
6451 AliInfo(Form("fRebin: %d",fRebin));
6453 for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
6455 for (Int_t ji = i; ji < i+fRebin; ji++) {
6456 sum += hist->GetBinContent(ji);
6459 rehist->SetBinContent(k,sum);
6467 //____________Some basic geometry function_____________________________________
6470 //_____________________________________________________________________________
6471 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
6474 // Reconstruct the plane number from the detector number
6477 return ((Int_t) (d % 6));
6481 //_____________________________________________________________________________
6482 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
6485 // Reconstruct the stack number from the detector number
6487 const Int_t kNlayer = 6;
6489 return ((Int_t) (d % 30) / kNlayer);
6493 //_____________________________________________________________________________
6494 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
6497 // Reconstruct the sector number from the detector number
6501 return ((Int_t) (d / fg));
6506 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
6508 //_______________________________________________________________________________
6509 void AliTRDCalibraFit::ResetVectorFit()
6512 // Reset the VectorFits
6515 fVectorFit.SetOwner();
6517 fVectorFit2.SetOwner();
6518 fVectorFit2.Clear();
6522 //____________Private Functions________________________________________________
6525 //_____________________________________________________________________________
6526 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
6529 // Function for the fit
6532 //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
6534 //PARAMETERS FOR FIT PH
6536 //fAsymmGauss->SetParameter(0,0.113755);
6537 //fAsymmGauss->SetParameter(1,0.350706);
6538 //fAsymmGauss->SetParameter(2,0.0604244);
6539 //fAsymmGauss->SetParameter(3,7.65596);
6540 //fAsymmGauss->SetParameter(4,1.00124);
6541 //fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
6549 Double_t dx = 0.005;
6550 Double_t xs = par[1];
6552 Double_t paras[2] = { 0.0, 0.0 };
6555 if ((xs >= par[1]) &&
6556 (xs < (par[1]+par[2]))) {
6557 //fAsymmGauss->SetParameter(0,par[0]);
6558 //fAsymmGauss->SetParameter(1,xs);
6559 //ss += fAsymmGauss->Eval(xx);
6562 ss += AsymmGauss(&xx,paras);
6564 if ((xs >= (par[1]+par[2])) &&
6565 (xs < (par[1]+par[2]+par[3]))) {
6566 //fAsymmGauss->SetParameter(0,par[0]*par[4]);
6567 //fAsymmGauss->SetParameter(1,xs);
6568 //ss += fAsymmGauss->Eval(xx);
6569 paras[0] = par[0]*par[4];
6571 ss += AsymmGauss(&xx,paras);
6580 //_____________________________________________________________________________
6581 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
6584 // Function for the fit
6587 //par[0] = normalization
6595 Double_t par1save = par[1];
6596 //Double_t par2save = par[2];
6597 Double_t par2save = 0.0604244;
6598 //Double_t par3save = par[3];
6599 Double_t par3save = 7.65596;
6600 //Double_t par5save = par[5];
6601 Double_t par5save = 0.870597;
6602 Double_t dx = x[0] - par1save;
6604 Double_t sigma2 = par2save*par2save;
6605 Double_t sqrt2 = TMath::Sqrt(2.0);
6606 Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
6607 * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
6608 Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
6609 * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
6611 //return par[0]*(exp1+par[4]*exp2);
6612 return par[0] * (exp1 + 1.00124 * exp2);
6616 //_____________________________________________________________________________
6617 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
6620 // Sum Landau + Gaus with identical mean
6623 Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
6624 //Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6625 Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6626 Double_t val = valLandau + valGaus;
6632 //_____________________________________________________________________________
6633 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6636 // Function for the fit
6639 // par[0]=Width (scale) parameter of Landau density
6640 // par[1]=Most Probable (MP, location) parameter of Landau density
6641 // par[2]=Total area (integral -inf to inf, normalization constant)
6642 // par[3]=Width (sigma) of convoluted Gaussian function
6644 // In the Landau distribution (represented by the CERNLIB approximation),
6645 // the maximum is located at x=-0.22278298 with the location parameter=0.
6646 // This shift is corrected within this function, so that the actual
6647 // maximum is identical to the MP parameter.
6650 // Numeric constants
6651 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6652 Double_t mpshift = -0.22278298; // Landau maximum location
6654 // Control constants
6655 Double_t np = 100.0; // Number of convolution steps
6656 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6661 Double_t fland = 0.0;
6663 Double_t xlow = 0.0;
6664 Double_t xupp = 0.0;
6665 Double_t step = 0.0;
6668 // MP shift correction
6669 mpc = par[1] - mpshift * par[0];
6671 // Range of convolution integral
6672 xlow = x[0] - sc * par[3];
6673 xupp = x[0] + sc * par[3];
6675 step = (xupp - xlow) / np;
6677 // Convolution integral of Landau and Gaussian by sum
6678 for (i = 1.0; i <= np/2; i++) {
6680 xx = xlow + (i-.5) * step;
6681 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6682 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6684 xx = xupp - (i-.5) * step;
6685 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6686 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6690 if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
6694 //_____________________________________________________________________________
6695 Double_t AliTRDCalibraFit::LanGauFunEx(const Double_t *x, const Double_t *par)
6698 // Function for the fit
6701 // par[0]=Width (scale) parameter of Landau density
6702 // par[1]=Most Probable (MP, location) parameter of Landau density
6703 // par[2]=Total area (integral -inf to inf, normalization constant)
6704 // par[3]=Width (sigma) of convoluted Gaussian function
6705 // par[4]=Exponential Slope Parameter
6707 // In the Landau distribution (represented by the CERNLIB approximation),
6708 // the maximum is located at x=-0.22278298 with the location parameter=0.
6709 // This shift is corrected within this function, so that the actual
6710 // maximum is identical to the MP parameter.
6713 // Numeric constants
6714 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
6715 Double_t mpshift = -0.22278298; // Landau maximum location
6717 // Control constants
6718 Double_t np = 100.0; // Number of convolution steps
6719 Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
6724 Double_t fland = 0.0;
6731 // MP shift correction
6732 mpc = par[1] - mpshift * par[0];
6734 // Range of convolution integral
6735 xlow = x[0] - sc * par[3];
6736 xupp = x[0] + sc * par[3];
6738 step = (xupp - xlow) / np;
6740 // Convolution integral of Landau and Gaussian by sum
6741 for (i = 1.0; i <= np/2; i++) {
6743 xx = xlow + (i-.5) * step;
6744 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
6745 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6747 xx = xupp - (i-.5) * step;
6748 if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
6749 sum += fland * TMath::Gaus(x[0],xx,par[3]);
6753 if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
6757 //_____________________________________________________________________________
6758 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6759 , const Double_t *parlimitslo, const Double_t *parlimitshi
6760 , Double_t *fitparams, Double_t *fiterrors
6761 , Double_t *chiSqr, Int_t *ndf) const
6764 // Function for the fit
6768 Char_t funname[100];
6770 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6775 TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6776 ffit->SetParameters(startvalues);
6777 ffit->SetParNames("Width","MP","Area","GSigma");
6779 for (i = 0; i < 4; i++) {
6780 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6783 his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
6785 ffit->GetParameters(fitparams); // Obtain fit parameters
6786 for (i = 0; i < 4; i++) {
6787 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6789 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6790 ndf[0] = ffit->GetNDF(); // Obtain ndf
6792 return (ffit); // Return fit function
6795 //_____________________________________________________________________________
6796 TF1 *AliTRDCalibraFit::LanGauFitEx(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6797 , const Double_t *parlimitslo, const Double_t *parlimitshi
6798 , Double_t *fitparams, Double_t *fiterrors
6799 , Double_t *chiSqr, Int_t *ndf) const
6802 // Function for the fit
6806 Char_t funname[100];
6808 TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6813 TF1 *ffit = new TF1(funname,LanGauFunEx,fitrange[0],fitrange[1],5);
6814 ffit->SetParameters(startvalues);
6815 ffit->SetParNames("Width","MP","Area","GSigma","Ex");
6817 for (i = 0; i < 5; i++) {
6818 ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6821 his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
6823 ffit->GetParameters(fitparams); // Obtain fit parameters
6824 for (i = 0; i < 5; i++) {
6825 fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
6827 chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
6828 ndf[0] = ffit->GetNDF(); // Obtain ndf
6830 return (ffit); // Return fit function