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 /////////////////////////////////////////////////////////////////////////////////
20 // AliTRDCalibraFillHisto
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 fills histos or vectors.
24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
25 // The user has to choose with the functions SetNz and SetNrphi the precision of the
26 // calibration (see AliTRDCalibraMode).
27 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
28 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
29 // in the function "FollowBackProlongation" (AliTRDtracker)
30 // Per default the functions to fill are off.
33 // R. Bailhache (R.Bailhache@gsi.de)
35 //////////////////////////////////////////////////////////////////////////////////////
38 #include <TProfile2D.h>
44 #include <TGraphErrors.h>
45 #include <TObjArray.h>
49 #include <TStopwatch.h>
51 #include <TDirectory.h>
55 #include "AliCDBManager.h"
57 #include "AliTRDCalibraFillHisto.h"
58 #include "AliTRDCalibraMode.h"
59 #include "AliTRDCalibraVector.h"
60 #include "AliTRDcalibDB.h"
61 #include "AliTRDCommonParam.h"
62 #include "AliTRDmcmTracklet.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcluster.h"
65 #include "AliTRDtrack.h"
68 ClassImp(AliTRDCalibraFillHisto)
70 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
71 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
73 //_____________singleton implementation_________________________________________________
74 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
77 // Singleton implementation
80 if (fgTerminated != kFALSE) {
84 if (fgInstance == 0) {
85 fgInstance = new AliTRDCalibraFillHisto();
92 //______________________________________________________________________________________
93 void AliTRDCalibraFillHisto::Terminate()
96 // Singleton implementation
97 // Deletes the instance of this class
100 fgTerminated = kTRUE;
102 if (fgInstance != 0) {
109 //______________________________________________________________________________________
110 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
113 ,fMcmTracking(kFALSE)
114 ,fMcmCorrectAngle(kFALSE)
121 ,fCountRelativeScale(0)
122 ,fRelativeScaleAuto(kFALSE)
123 ,fThresholdClusterPRF1(0.0)
124 ,fThresholdClusterPRF2(0.0)
125 ,fCenterOfflineCluster(kFALSE)
128 ,fDetectorAliTRDtrack(kFALSE)
129 ,fChamberAliTRDtrack(-1)
130 ,fDetectorPreviousTrack(-1)
149 // Default constructor
152 fCalibraMode = new AliTRDCalibraMode();
155 for (Int_t i = 0; i < 3; i++) {
164 //______________________________________________________________________________________
165 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
168 ,fMcmTracking(kFALSE)
169 ,fMcmCorrectAngle(kFALSE)
176 ,fCountRelativeScale(0)
177 ,fRelativeScaleAuto(kFALSE)
178 ,fThresholdClusterPRF1(0.0)
179 ,fThresholdClusterPRF2(0.0)
180 ,fCenterOfflineCluster(kFALSE)
183 ,fDetectorAliTRDtrack(kFALSE)
184 ,fChamberAliTRDtrack(-1)
185 ,fDetectorPreviousTrack(-1)
209 //____________________________________________________________________________________
210 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
213 // AliTRDCalibraFillHisto destructor
220 //_____________________________________________________________________________
221 void AliTRDCalibraFillHisto::Destroy()
234 //_____________________________________________________________________________
235 void AliTRDCalibraFillHisto::ClearHistos()
256 //_____________________________________________________________________________
257 void AliTRDCalibraFillHisto::Init()
260 // Init some default values
263 // How to fill the 2D
264 fThresholdClusterPRF1 = 2.0;
265 fThresholdClusterPRF2 = 15.0;
268 fNumberBinCharge = 100;
272 fWriteName = "TRD.calibration.root";
274 // Internal variables
276 // Fill the 2D histos in the offline tracking
277 fDetectorPreviousTrack = -1;
278 fChamberAliTRDtrack = -1;
283 fNumberClusters = 18;
285 fNumberUsedCh[0] = 0;
286 fNumberUsedCh[1] = 0;
287 fNumberUsedPh[0] = 0;
288 fNumberUsedPh[1] = 0;
292 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
293 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
296 // For the offline tracking
297 // This function will be called in the function AliReconstruction::Run()
298 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
303 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
305 AliInfo("Could not get calibDB");
308 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
310 AliInfo("Could not get CommonParam");
315 fTimeMax = cal->GetNumberOfTimeBins();
316 fSf = parCom->GetSamplingFrequency();
317 if (fRelativeScaleAuto) {
324 //If vector method On initialised all the stuff
326 fCalibraVector = new AliTRDCalibraVector();
330 // Create the 2D histos corresponding to the pad groupCalibration mode
333 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
334 ,fCalibraMode->GetNz(0)
335 ,fCalibraMode->GetNrphi(0)));
337 // Calcul the number of Xbins
339 fCalibraMode->ModePadCalibration(2,0);
340 fCalibraMode->ModePadFragmentation(0,2,0,0);
341 fCalibraMode->SetDetChamb2(0);
342 Ntotal0 += 6 * 18 * fCalibraMode->GetDetChamb2(0);
343 fCalibraMode->ModePadCalibration(0,0);
344 fCalibraMode->ModePadFragmentation(0,0,0,0);
345 fCalibraMode->SetDetChamb0(0);
346 Ntotal0 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(0);
347 AliInfo(Form("Total number of Xbins: %d",Ntotal0));
349 // Create the 2D histo
354 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
358 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
359 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
367 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
368 ,fCalibraMode->GetNz(1)
369 ,fCalibraMode->GetNrphi(1)));
371 // Calcul the number of Xbins
373 fCalibraMode->ModePadCalibration(2,1);
374 fCalibraMode->ModePadFragmentation(0,2,0,1);
375 fCalibraMode->SetDetChamb2(1);
376 Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
377 fCalibraMode->ModePadCalibration(0,1);
378 fCalibraMode->ModePadFragmentation(0,0,0,1);
379 fCalibraMode->SetDetChamb0(1);
380 Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
381 AliInfo(Form("Total number of Xbins: %d",Ntotal1));
383 // Create the 2D histo
388 fCalibraVector->SetTimeMax(fTimeMax);
392 fPHPlace = new Short_t[fTimeMax];
393 for (Int_t k = 0; k < fTimeMax; k++) {
396 fPHValue = new Float_t[fTimeMax];
397 for (Int_t k = 0; k < fTimeMax; k++) {
405 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
406 ,fCalibraMode->GetNz(2)
407 ,fCalibraMode->GetNrphi(2)));
409 // Calcul the number of Xbins
411 fCalibraMode->ModePadCalibration(2,2);
412 fCalibraMode->ModePadFragmentation(0,2,0,2);
413 fCalibraMode->SetDetChamb2(2);
414 Ntotal2 += 6 * 18 * fCalibraMode->GetDetChamb2(2);
415 fCalibraMode->ModePadCalibration(0,2);
416 fCalibraMode->ModePadFragmentation(0,0,0,2);
417 fCalibraMode->SetDetChamb0(2);
418 Ntotal2 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(2);
419 AliInfo(Form("Total number of Xbins: %d",Ntotal2));
421 // Create the 2D histo
423 CreatePRF2d(Ntotal2);
426 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
435 //____________Functions for filling the histos in the code_____________________
437 //____________Offine tracking in the AliTRDtracker_____________________________
438 Bool_t AliTRDCalibraFillHisto::ResetTrack()
441 // For the offline tracking
442 // This function will be called in the function
443 // AliTRDtracker::FollowBackPropagation() at the beginning
444 // Reset the parameter to know we have a new TRD track
447 fDetectorAliTRDtrack = kFALSE;
452 //____________Offline tracking in the AliTRDtracker____________________________
453 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
456 // For the offline tracking
457 // This function will be called in the function
458 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
460 // Fill the 2D histos or the vectors with the info of the clusters at
461 // the end of a detectors if the track is "good"
464 // Get the parameter object
465 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
467 AliInfo("Could not get CommonParam");
471 // Get the parameter object
472 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
474 AliInfo("Could not get calibDB");
478 // Localisation of the detector
479 Int_t detector = cl->GetDetector();
480 Int_t chamber = GetChamber(detector);
481 Int_t plane = GetPlane(detector);
483 // Fill the infos for the previous clusters if not the same
484 // detector anymore or if not the same track
485 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
486 (fDetectorPreviousTrack != -1)) {
490 // If the same track, then look if the previous detector is in
491 // the same plane, if yes: not a good track
492 if (fDetectorAliTRDtrack &&
493 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
497 // Fill only if the track doesn't touch a masked pad or doesn't
498 // appear in the middle (fGoodTrack)
503 FillTheInfoOfTheTrackCH();
508 FillTheInfoOfTheTrackPH();
515 } // Fill at the end the charge
517 // Calcul the position of the detector
518 if (detector != fDetectorPreviousTrack) {
519 LocalisationDetectorXbins(detector);
522 // Reset the good track for the PRF
525 // Localisation of the cluster
526 Double_t pos[3] = { 0.0, 0.0, 0.0 };
530 Int_t time = cl->GetLocalTimeBin();
532 // Reset the detector
533 fDetectorPreviousTrack = detector;
534 fDetectorAliTRDtrack = kTRUE;
536 // Position of the cluster
537 AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
538 Int_t row = padplane->GetPadRowNumber(pos[2]);
539 Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
540 Double_t offsettilt = padplane->GetTiltOffset(offsetz);
541 Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt);
543 // See if we are not near a masked pad
544 if (!IsPadOn(detector,col,row)) {
550 if (!IsPadOn(detector,col-1,row)) {
557 if (!IsPadOn(detector,col+1,row)) {
563 // Row of the cluster and position in the pad groups
564 Int_t posr[3] = { 0, 0, 0 };
565 if ((fCH2dOn) && (fCalibraMode->GetNnZ(0) != 0)) {
566 posr[0] = (Int_t) row / fCalibraMode->GetNnZ(0);
568 if ((fPH2dOn) && (fCalibraMode->GetNnZ(1) != 0)) {
569 posr[1] = (Int_t) row / fCalibraMode->GetNnZ(1);
571 if ((fPRF2dOn) && (fCalibraMode->GetNnZ(2) != 0)) {
572 posr[2] = (Int_t) row / fCalibraMode->GetNnZ(2);
575 // Col of the cluster and position in the pad groups
576 Int_t posc[3] = { 0, 0, 0 };
577 if ((fCH2dOn) && (fCalibraMode->GetNnRphi(0) != 0)) {
578 posc[0] = (Int_t) col / fCalibraMode->GetNnRphi(0);
580 if ((fPH2dOn) && (fCalibraMode->GetNnRphi(1) != 0)) {
581 posc[1] = (Int_t) col / fCalibraMode->GetNnRphi(1);
583 if ((fPRF2dOn) && (fCalibraMode->GetNnRphi(2) != 0)) {
584 posc[2] = (Int_t) col / fCalibraMode->GetNnRphi(2);
587 // Charge in the cluster
588 // For the moment take the abs
589 Float_t q = TMath::Abs(cl->GetQ());
590 Short_t *signals = cl->GetSignals();
592 // Correction due to the track angle
593 Float_t correction = 1.0;
594 Float_t normalisation = 6.67;
595 if ((q >0) && (t->GetNdedx() > 0)) {
596 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
599 // Fill the fAmpTotal with the charge
601 fAmpTotal[(Int_t) (posc[0]*fCalibraMode->GetNfragZ(0)+posr[0])] += correction;
604 // Fill the fPHPlace and value
606 fPHPlace[time] = posc[1]*fCalibraMode->GetNfragZ(1)+posr[1];
607 fPHValue[time] = correction;
610 // Fill direct the PRF
611 if ((fPRF2dOn) && (good)) {
613 Float_t yminus = 0.0;
614 Float_t xcenter = 0.0;
615 Float_t ycenter = 0.0;
617 Bool_t echec = kFALSE;
619 if ((cl->From3pad()) && (!cl->IsUsed())) {
621 // Center 3 balanced and cut on the cluster shape
622 if ((((Float_t) signals[3]) > fThresholdClusterPRF2) &&
623 (((Float_t) signals[2]) > fThresholdClusterPRF2) &&
624 (((Float_t) signals[4]) > fThresholdClusterPRF2) &&
625 (((Float_t) signals[1]) < fThresholdClusterPRF1) &&
626 (((Float_t) signals[5]) < fThresholdClusterPRF1) &&
627 ((((Float_t) signals[2])*((Float_t) signals[4])/(((Float_t) signals[3])*((Float_t) signals[3]))) < 0.06)) {
630 //First calculate the position of the cluster and the y
631 //echec enables to repair cases where it fails
633 // Col correspond to signals[3]
634 if (fCenterOfflineCluster) {
635 xcenter = cl->GetCenter();
638 // Security if the denomiateur is 0
639 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
640 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
641 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
642 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
643 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
649 //after having calculating the position calculate the y
650 if (TMath::Abs(xcenter) < 0.5) {
651 ycenter = (Float_t) (((Float_t) signals[3])
652 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
653 yminus = (Float_t) (((Float_t) signals[2])
654 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
655 ymax = (Float_t) (((Float_t) signals[4])
656 / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
657 //If the charge of the cluster is too far away from the corrected one cut
658 if ((TMath::Abs(((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4])) - q) > 10.0)) {
666 //Then Fill the histo if no echec
668 // Fill only if it is in the drift region!
669 if ((((Float_t) (((Float_t) time) / fSf)) > 0.3) && (!echec)) {
671 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),xcenter,ycenter);
672 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),-(xcenter+1.0),yminus);
673 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),1.0-xcenter,ymax);
676 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],xcenter,ycenter);
677 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],-(xcenter+1.0),yminus);
678 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],1.0-xcenter,ymax);
680 } // If in the drift region
681 } // center 3 balanced and cut on the cluster shape
689 //____________Online trackling in AliTRDtrigger________________________________
690 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
694 // This function will be called in the function AliTRDtrigger::TestTracklet
695 // before applying the pt cut on the tracklets
696 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
699 // Localisation of the Xbins involved
700 Int_t idect = trk->GetDetector();
701 LocalisationDetectorXbins(idect);
703 // Get the parameter object
704 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
706 AliInfo("Could not get calibDB");
713 // Row of the tracklet and position in the pad groups
714 Int_t row = trk->GetRow();
715 Int_t posr[3] = { 0, 0, 0 };
716 if ((fCH2dOn) && (fCalibraMode->GetNnZ(0) != 0)) {
717 posr[0] = (Int_t) row / fCalibraMode->GetNnZ(0);
719 if ((fPH2dOn) && (fCalibraMode->GetNnZ(1) != 0)) {
720 posr[1] = (Int_t) row / fCalibraMode->GetNnZ(1);
722 if ((fPRF2dOn) && (fCalibraMode->GetNnZ(2) != 0)) {
723 posr[2] = (Int_t) row / fCalibraMode->GetNnZ(2);
726 // Eventuelle correction due to track angle in z direction
727 Float_t correction = 1.0;
728 if (fMcmCorrectAngle) {
729 Float_t z = trk->GetRowz();
730 Float_t r = trk->GetTime0();
731 correction = r / TMath::Sqrt((r*r+z*z));
734 // Boucle sur les clusters
735 // Condition on number of cluster: don't come from the middle of the detector
736 if (trk->GetNclusters() >= fNumberClusters) {
738 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
740 Float_t amp[3] = { 0.0, 0.0, 0.0 };
741 Int_t time = trk->GetClusterTime(icl);
742 Int_t col = trk->GetClusterCol(icl);
744 amp[0] = trk->GetClusterADC(icl)[0] * correction;
745 amp[1] = trk->GetClusterADC(icl)[1] * correction;
746 amp[2] = trk->GetClusterADC(icl)[2] * correction;
749 if ((amp[0] < 0.0) ||
755 // Col of cluster and position in the pad groups
756 Int_t posc[3] = { 0, 0, 0 };
757 if ((fCH2dOn) && (fCalibraMode->GetNnRphi(0) != 0)) {
758 posc[0] = (Int_t) col / fCalibraMode->GetNnRphi(0);
760 if ((fPH2dOn) && (fCalibraMode->GetNnRphi(1) != 0)) {
761 posc[1] = (Int_t) col / fCalibraMode->GetNnRphi(1);
763 if ((fPRF2dOn) && (fCalibraMode->GetNnRphi(2) != 0)) {
764 posc[2] = (Int_t) col / fCalibraMode->GetNnRphi(2);
767 // See if we are not near a masked pad
769 if (!IsPadOn(idect,col,row)) {
775 if (!IsPadOn(idect,col-1,row)) {
782 if (!IsPadOn(idect,col+1,row)) {
790 fPHPlace[time] = posc[1] * fCalibraMode->GetNfragZ(1) + posr[1];
794 fAmpTotal[(Int_t) (posc[0]*fCalibraMode->GetNfragZ(0)+posr[0])] += (Float_t) (amp[0]+amp[1]+amp[2]);
797 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
802 if (fPRF2dOn && good) {
804 if ((amp[0] > fThresholdClusterPRF2) &&
805 (amp[1] > fThresholdClusterPRF2) &&
806 (amp[2] > fThresholdClusterPRF2) &&
807 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
809 // Security of the denomiateur is 0
810 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
811 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
812 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
813 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
814 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
816 if (TMath::Abs(xcenter) < 0.5) {
817 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
818 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
819 // Fill only if it is in the drift region!
820 if (((Float_t) time / fSf) > 0.3) {
822 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),xcenter,ycenter);
823 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),-(xcenter+1.0),yminus);
824 fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),(1.0-xcenter),ymax);
827 fCalibraVector->UpdateVectorPRF((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]),xcenter,ycenter);
828 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],-(xcenter+1.0),yminus);
829 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],(1.0-xcenter),ymax);
831 }//in the drift region
833 }//denominateur security
834 }//cluster shape and thresholds
840 if (fCH2dOn && fGoodTrack) {
841 FillTheInfoOfTheTrackCH();
845 if (fPH2dOn && fGoodTrack) {
846 FillTheInfoOfTheTrackPH();
851 } // Condition on number of clusters
857 //_____________________________________________________________________________
858 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
861 // Look in the choosen database if the pad is On.
862 // If no the track will be "not good"
865 // Get the parameter object
866 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
868 AliInfo("Could not get calibDB");
872 if (!cal->IsChamberInstalled(detector) ||
873 cal->IsChamberMasked(detector) ||
874 cal->IsPadMasked(detector,col,row)) {
883 //____________Functions for plotting the 2D____________________________________
885 //_____________________________________________________________________________
886 void AliTRDCalibraFillHisto::Plot2d()
889 // Plot the 2D histos
893 TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
898 TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
903 TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
905 fPRF2d->Draw("LEGO");
910 //____________Writing the 2D___________________________________________________
912 //_____________________________________________________________________________
913 Bool_t AliTRDCalibraFillHisto::Write2d()
916 // Write the 2D histograms or the vectors converted in trees in the file
917 // "TRD.calibration.root"
920 TFile *fout = TFile::Open(fWriteName,"RECREATE");
921 // Check if the file could be opened
922 if (!fout || !fout->IsOpen()) {
923 AliInfo("No File found!");
926 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
933 TStopwatch stopwatch;
937 if ((fCH2dOn ) && (fWrite[0])) {
939 fout->WriteTObject(fCH2d);
943 name += fCalibraMode->GetNz(0);
945 name += fCalibraMode->GetNrphi(0);
946 TTree *treeCH2d = fCalibraVector->ConvertVectorCTTreeHisto(fCalibraVector->GetVectorCH(),fCalibraVector->GetPlaCH(),"treeCH2d",(const char *) name);
947 fout->WriteTObject(treeCH2d);
950 if ((fPH2dOn ) && (fWrite[1])) {
952 fout->WriteTObject(fPH2d);
956 name += fCalibraMode->GetNz(1);
958 name += fCalibraMode->GetNrphi(1);
959 TTree *treePH2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPH(),fCalibraVector->GetPlaPH(),"treePH2d",(const char *) name);
960 fout->WriteTObject(treePH2d);
963 if ((fPRF2dOn ) && (fWrite[2])) {
965 fout->WriteTObject(fPRF2d);
969 name += fCalibraMode->GetNz(2);
971 name += fCalibraMode->GetNrphi(2);
972 TTree *treePRF2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPRF(),fCalibraVector->GetPlaPRF(),"treePRF2d",(const char *) name);
973 fout->WriteTObject(treePRF2d);
979 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
980 ,stopwatch.RealTime(),stopwatch.CpuTime()));
986 //____________Probe the histos_________________________________________________
987 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
990 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
991 // debug mode with 2 for TH2I and 3 for TProfile2D
992 // It gives a pointer to a Double_t[7] with the info following...
993 // [0] : number of calibration groups with entries
994 // [1] : minimal number of entries found
995 // [2] : calibration group number of the min
996 // [3] : maximal number of entries found
997 // [4] : calibration group number of the max
998 // [5] : mean number of entries found
999 // [6] : mean relativ error
1002 Double_t *info = new Double_t[7];
1004 // Number of Xbins (detectors or groups of pads)
1005 Int_t nbins = h->GetNbinsX(); //number of calibration groups
1006 Int_t nybins = h->GetNbinsY(); //number of bins per histo
1009 Double_t nbwe = 0; //number of calibration groups with entries
1010 Double_t minentries = 0; //minimal number of entries found
1011 Double_t maxentries = 0; //maximal number of entries found
1012 Double_t placemin = 0; //calibration group number of the min
1013 Double_t placemax = -1; //calibration group number of the max
1014 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1015 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1017 Double_t counter = 0;
1020 TH1F *NbEntries = 0x0;//distribution of the number of entries
1021 TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1022 TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1024 // Beginning of the loop over the calibration groups
1025 for (Int_t idect = 0; idect < nbins; idect++) {
1027 TH1I *projch = (TH1I *) h->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
1028 projch->SetDirectory(0);
1030 // Number of entries for this calibration group
1031 Double_t nentries = 0.0;
1033 for (Int_t k = 0; k < nybins; k++) {
1034 nentries += h->GetBinContent(h->GetBin(idect+1,k+1));
1038 for (Int_t k = 0; k < nybins; k++) {
1039 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(idect+1,k+1));
1040 if(h->GetBinContent(h->GetBin(idect+1,k+1)) != 0) {
1041 meanrelativerror += (h->GetBinError(h->GetBin(idect+1,k+1))
1042 / (TMath::Abs(h->GetBinContent(h->GetBin(idect+1,k+1)))));
1050 if((!((Bool_t)NbEntries)) && (nentries > 0)){
1051 NbEntries = new TH1F("Number of entries","Number of entries"
1052 ,100,(Int_t)nentries/2,nentries*2);
1053 NbEntries->SetDirectory(0);
1054 NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1056 NbEntriesPerGroup->SetDirectory(0);
1057 NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1058 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1059 NbEntriesPerSp->SetDirectory(0);
1062 if(nentries > 0) NbEntries->Fill(nentries);
1063 NbEntriesPerGroup->Fill(idect+0.5,nentries);
1064 NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1069 if(nentries > maxentries){
1070 maxentries = nentries;
1074 minentries = nentries;
1076 if(nentries < minentries){
1077 minentries = nentries;
1083 meanstats += nentries;
1086 }//calibration groups loop
1088 if(nbwe > 0) meanstats /= nbwe;
1089 if(counter > 0) meanrelativerror /= counter;
1091 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1092 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1093 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1094 AliInfo(Form("The mean number of entries is %f",meanstats));
1095 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1098 info[1] = minentries;
1100 info[3] = maxentries;
1102 info[5] = meanstats;
1103 info[6] = meanrelativerror;
1106 gStyle->SetPalette(1);
1107 gStyle->SetOptStat(1111);
1108 gStyle->SetPadBorderMode(0);
1109 gStyle->SetCanvasColor(10);
1110 gStyle->SetPadLeftMargin(0.13);
1111 gStyle->SetPadRightMargin(0.01);
1112 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1115 NbEntries->Draw("");
1117 NbEntriesPerSp->SetStats(0);
1118 NbEntriesPerSp->Draw("");
1119 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1121 NbEntriesPerGroup->SetStats(0);
1122 NbEntriesPerGroup->Draw("");
1129 //_____________________________________________________________________________
1130 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1133 // Set the factor that will divide the deposited charge
1134 // to fit in the histo range [0,300]
1137 if (RelativeScale > 0.0) {
1138 fRelativeScale = RelativeScale;
1141 AliInfo("RelativeScale must be strict positif!");
1146 //_____________________________________________________________________________
1147 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1150 // Set the mode of calibration group in the z direction for the parameter i
1155 fCalibraMode->SetNz(i, Nz);
1158 AliInfo("You have to choose between 0 and 4");
1163 //_____________________________________________________________________________
1164 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1167 // Set the mode of calibration group in the rphi direction for the parameter i
1172 fCalibraMode->SetNrphi(i ,Nrphi);
1175 AliInfo("You have to choose between 0 and 6");
1180 //____________Protected Functions______________________________________________
1181 //____________Create the 2D histo to be filled online__________________________
1184 //_____________________________________________________________________________
1185 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1188 // Create the 2D histos
1192 name += fCalibraMode->GetNz(2);
1194 name += fCalibraMode->GetNrphi(2);
1196 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1197 ,nn,0,nn,fNumberBinPRF,-1.5,1.5);
1198 fPRF2d->SetXTitle("Det/pad groups");
1199 fPRF2d->SetYTitle("Position x/W [pad width units]");
1200 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1201 fPRF2d->SetStats(0);
1205 //_____________________________________________________________________________
1206 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1209 // Create the 2D histos
1213 name += fCalibraMode->GetNz(1);
1215 name += fCalibraMode->GetNrphi(1);
1217 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1219 ,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf);
1220 fPH2d->SetXTitle("Det/pad groups");
1221 fPH2d->SetYTitle("time [#mus]");
1222 fPH2d->SetZTitle("<PH> [a.u.]");
1227 //_____________________________________________________________________________
1228 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1231 // Create the 2D histos
1235 name += fCalibraMode->GetNz(0);
1237 name += fCalibraMode->GetNrphi(0);
1239 fCH2d = new TH2I("CH2d",(const Char_t *) name
1240 ,nn,0,nn,fNumberBinCharge,0,300);
1241 fCH2d->SetXTitle("Det/pad groups");
1242 fCH2d->SetYTitle("charge deposit [a.u]");
1243 fCH2d->SetZTitle("counts");
1249 //____________Offine tracking in the AliTRDtracker_____________________________
1250 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1253 // For the offline tracking or mcm tracklets
1254 // This function will be called in the functions UpdateHistogram...
1255 // to fill the info of a track for the relativ gain calibration
1258 Int_t nb = 0; // Nombre de zones traversees
1259 Int_t fd = -1; // Premiere zone non nulle
1262 // See if the track goes through different zones
1263 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1264 if (fAmpTotal[k] > 0.0) {
1272 // If automatic scale
1273 if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
1274 // Take only the one zone track
1276 fRelativeScale += fAmpTotal[fd] * 0.014 * 0.01;
1277 fCountRelativeScale++;
1281 // We fill the CH2d after having scale with the first 100
1282 if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
1283 // Case of track with only one zone
1286 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1289 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1292 // Case of track with two zones
1294 // Two zones voisines sinon rien!
1295 if ((fAmpTotal[fd] > 0.0) &&
1296 (fAmpTotal[fd+1] > 0.0)) {
1297 // One of the two very big
1298 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1300 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1303 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1306 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1308 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
1311 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd+1]/fRelativeScale);
1318 // Fill with no automatic scale
1319 if (!fRelativeScaleAuto) {
1320 // Case of track with only one zone
1324 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1327 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1330 // Case of track with two zones
1332 // Two zones voisines sinon rien!
1334 if ((fAmpTotal[fd] > 0.0) &&
1335 (fAmpTotal[fd+1] > 0.0)) {
1336 // One of the two very big
1337 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1339 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1342 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1346 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1348 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
1351 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1357 if (fCalibraMode->GetNfragZ(0) > 1) {
1358 if (fAmpTotal[fd] > 0.0) {
1359 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1360 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1361 // One of the two very big
1362 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1364 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1367 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1371 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1373 fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)
1374 + 0.5,fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1378 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)
1379 ,fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1392 //____________Offine tracking in the AliTRDtracker_____________________________
1393 void AliTRDCalibraFillHisto::ResetfVariables()
1396 // Reset values of fAmpTotal, fPHValue and fPHPlace for
1397 // the updateHistogram... functions
1400 // Reset the good track
1403 // Reset the fAmpTotal where we put value
1405 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1410 // Reset the fPHValue
1412 for (Int_t k = 0; k < fTimeMax; k++) {
1420 //____________Offine tracking in the AliTRDtracker_____________________________
1421 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1424 // For the offline tracking or mcm tracklets
1425 // This function will be called in the functions UpdateHistogram...
1426 // to fill the info of a track for the drift velocity calibration
1429 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1430 Int_t fd1 = -1; // Premiere zone non nulle
1431 Int_t fd2 = -1; // Deuxieme zone non nulle
1432 Int_t k1 = -1; // Debut de la premiere zone
1433 Int_t k2 = -1; // Debut de la seconde zone
1435 // See if the track goes through different zones
1436 for (Int_t k = 0; k < fTimeMax; k++) {
1437 if (fPHValue[k] > 0.0) {
1442 if (fPHPlace[k] != fd1) {
1448 if (fPHPlace[k] != fd2) {
1456 // Case of track with only one zone
1459 //fd1 is the only zone
1460 for (Int_t i = 0; i < fTimeMax; i++) {
1462 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1465 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1469 // Case of track with two zones
1471 // Two zones voisines sinon rien!
1473 if ((fd1 == fd2+1) ||
1475 // One of the two fast all the think
1476 if (k2 > (k1+fDifference)) {
1477 //we choose to fill the fd1 with all the values
1479 for (Int_t i = 0; i < fTimeMax; i++) {
1481 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1484 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1488 if ((k2+fDifference) < fTimeMax) {
1489 //we choose to fill the fd2 with all the values
1491 for (Int_t i = 0; i < fTimeMax; i++) {
1493 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1496 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1501 // Two zones voisines sinon rien!
1502 if (fCalibraMode->GetNfragZ(1) > 1) {
1504 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1505 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1506 // One of the two fast all the think
1507 if (k2 > (k1+fDifference)) {
1508 //we choose to fill the fd1 with all the values
1510 for (Int_t i = 0; i < fTimeMax; i++) {
1512 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1515 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1519 if ((k2+fDifference) < fTimeMax) {
1520 //we choose to fill the fd2 with all the values
1522 for (Int_t i = 0; i < fTimeMax; i++) {
1524 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1527 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1533 // Two zones voisines sinon rien!
1535 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1536 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1537 // One of the two fast all the think
1538 if (k2 > (k1 + fDifference)) {
1539 //we choose to fill the fd1 with all the values
1541 for (Int_t i = 0; i < fTimeMax; i++) {
1543 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1546 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1550 if ((k2+fDifference) < fTimeMax) {
1551 //we choose to fill the fd2 with all the values
1553 for (Int_t i = 0; i < fTimeMax; i++) {
1555 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1558 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1570 //____________Set the pad calibration variables for the detector_______________
1571 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1574 // For the detector calcul the first Xbins and set the number of row
1575 // and col pads per calibration groups, the number of calibration
1576 // groups in the detector.
1579 // first Xbins of the detector
1581 fCalibraMode->CalculXBins(detector,0);
1584 fCalibraMode->CalculXBins(detector,1);
1587 fCalibraMode->CalculXBins(detector,2);
1590 // fragmentation of idect
1591 for (Int_t i = 0; i < 3; i++) {
1592 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1593 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1594 , (Int_t) GetChamber(detector)
1595 , (Int_t) GetSector(detector),i);
1603 //____________Some basic geometry function_____________________________________
1606 //_____________________________________________________________________________
1607 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
1610 // Reconstruct the plane number from the detector number
1613 return ((Int_t) (d % 6));
1617 //_____________________________________________________________________________
1618 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
1621 // Reconstruct the chamber number from the detector number
1625 return ((Int_t) (d % 30) / fgkNplan);
1629 //_____________________________________________________________________________
1630 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
1633 // Reconstruct the sector number from the detector number
1637 return ((Int_t) (d / fg));