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 calibration (see AliTRDCalibraMode).
26 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
27 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
28 // in the function "FollowBackProlongation" (AliTRDtracker)
29 // Per default the functions to fill are off.
32 // R. Bailhache (R.Bailhache@gsi.de)
34 //////////////////////////////////////////////////////////////////////////////////////
37 #include <TProfile2D.h>
43 #include <TGraphErrors.h>
44 #include <TObjArray.h>
48 #include <TStopwatch.h>
50 #include <TDirectory.h>
52 #include <TTreeStream.h>
56 #include "AliCDBManager.h"
58 #include "AliTRDCalibraFillHisto.h"
59 #include "AliTRDCalibraMode.h"
60 #include "AliTRDCalibraVector.h"
61 #include "AliTRDcalibDB.h"
62 #include "AliTRDCommonParam.h"
63 #include "AliTRDmcmTracklet.h"
64 #include "AliTRDpadPlane.h"
65 #include "AliTRDcluster.h"
66 #include "AliTRDtrack.h"
67 #include "AliTRDRawStream.h"
68 #include "AliRawReader.h"
69 #include "AliRawReaderDate.h"
76 ClassImp(AliTRDCalibraFillHisto)
78 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
79 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
81 //_____________singleton implementation_________________________________________________
82 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
85 // Singleton implementation
88 if (fgTerminated != kFALSE) {
92 if (fgInstance == 0) {
93 fgInstance = new AliTRDCalibraFillHisto();
100 //______________________________________________________________________________________
101 void AliTRDCalibraFillHisto::Terminate()
104 // Singleton implementation
105 // Deletes the instance of this class
108 fgTerminated = kTRUE;
110 if (fgInstance != 0) {
116 //______________________________________________________________________________________
117 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
120 ,fMcmTracking(kFALSE)
121 ,fMcmCorrectAngle(kFALSE)
127 ,fLinearFitterOn(kFALSE)
128 ,fLinearFitterDebugOn(kFALSE)
130 ,fThresholdClusterPRF2(15.0)
131 ,fCalibraMode(new AliTRDCalibraMode())
134 ,fDetectorAliTRDtrack(kFALSE)
135 ,fDetectorPreviousTrack(-1)
142 ,fNumberBinCharge(100)
145 ,fListClusters(new TObjArray())
154 ,fGoodTracklet(kTRUE)
157 ,fEntriesLinearFitter(0x0)
162 ,fLinearFitterArray(0)
163 ,fLinearFitterHistoArray(0)
166 // Default constructor
170 // Init some default values
173 fNumberUsedCh[0] = 0;
174 fNumberUsedCh[1] = 0;
175 fNumberUsedPh[0] = 0;
176 fNumberUsedPh[1] = 0;
179 //______________________________________________________________________________________
180 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
182 ,fMITracking(c.fMITracking)
183 ,fMcmTracking(c.fMcmTracking)
184 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
187 ,fPRF2dOn(c.fPRF2dOn)
188 ,fHisto2d(c.fHisto2d)
189 ,fVector2d(c.fVector2d)
190 ,fLinearFitterOn(c.fLinearFitterOn)
191 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
192 ,fRelativeScale(c.fRelativeScale)
193 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
196 ,fDebugLevel(c.fDebugLevel)
197 ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
198 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
199 ,fNumberClusters(c.fNumberClusters)
200 ,fProcent(c.fProcent)
201 ,fDifference(c.fDifference)
202 ,fNumberTrack(c.fNumberTrack)
203 ,fTimeMax(c.fTimeMax)
205 ,fNumberBinCharge(c.fNumberBinCharge)
206 ,fNumberBinPRF(c.fNumberBinPRF)
207 ,fNgroupprf(c.fNgroupprf)
208 ,fListClusters(new TObjArray())
214 ,fAmpTotal(c.fAmpTotal)
215 ,fPHPlace(c.fPHPlace)
216 ,fPHValue(c.fPHValue)
217 ,fGoodTracklet(c.fGoodTracklet)
218 ,fGoodTrack(c.fGoodTrack)
219 ,fEntriesCH(c.fEntriesCH)
220 ,fEntriesLinearFitter(fEntriesLinearFitter)
225 ,fLinearFitterArray(0)
226 ,fLinearFitterHistoArray(0)
231 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
232 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
234 fPH2d = new TProfile2D(*c.fPH2d);
235 fPH2d->SetDirectory(0);
238 fPRF2d = new TProfile2D(*c.fPRF2d);
239 fPRF2d->SetDirectory(0);
242 fCH2d = new TH2I(*c.fCH2d);
243 fCH2d->SetDirectory(0);
245 if(c.fLinearFitterOn){
246 fLinearFitterArray.Expand(540);
247 for (Int_t cb = 0; cb < 540; cb++){
248 const TLinearFitter *linearFitter = (TLinearFitter*)c.fLinearFitterArray.UncheckedAt(cb);
249 if ( linearFitter != 0x0 ) fLinearFitterArray.AddAt(new TLinearFitter(*linearFitter), cb);
252 if(c.fLinearFitterDebugOn){
253 fLinearFitterHistoArray.Expand(540);
254 for (Int_t cb = 0; cb < 540; cb++){
255 const TH2F *linearfitterhisto = (TH2F*)c.fLinearFitterHistoArray.UncheckedAt(cb);
256 if ( linearfitterhisto != 0x0 ){
257 TH2F *hNew = new TH2F(*linearfitterhisto);
258 hNew->SetDirectory(0);
259 fLinearFitterHistoArray.AddAt(hNew,cb);
264 //____________________________________________________________________________________
265 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
268 // AliTRDCalibraFillHisto destructor
272 if ( fDebugStreamer ) delete fDebugStreamer;
276 //_____________________________________________________________________________
277 void AliTRDCalibraFillHisto::Destroy()
290 //_____________________________________________________________________________
291 void AliTRDCalibraFillHisto::ClearHistos()
311 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
312 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
315 // For the offline tracking
316 // This function will be called in the function AliReconstruction::Run()
317 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
322 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
324 AliInfo("Could not get calibDB");
327 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
329 AliInfo("Could not get CommonParam");
334 fTimeMax = cal->GetNumberOfTimeBins();
335 fSf = parCom->GetSamplingFrequency();
338 //Init the tracklet parameters
339 fPar0 = new Double_t[fTimeMax];
340 fPar1 = new Double_t[fTimeMax];
341 fPar2 = new Double_t[fTimeMax];
342 fPar3 = new Double_t[fTimeMax];
343 fPar4 = new Double_t[fTimeMax];
345 for(Int_t k = 0; k < fTimeMax; k++){
353 //If vector method On initialised all the stuff
355 fCalibraVector = new AliTRDCalibraVector();
356 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
357 fCalibraVector->SetTimeMax(fTimeMax);
358 if(fNgroupprf != 0) {
359 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
360 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
363 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
364 fCalibraVector->SetPRFRange(1.5);
368 // Create the 2D histos corresponding to the pad groupCalibration mode
371 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
372 ,fCalibraMode->GetNz(0)
373 ,fCalibraMode->GetNrphi(0)));
375 // Calcul the number of Xbins
376 Int_t Ntotal0 = CalculateTotalNumberOfBins(0);
377 // Create the 2D histo
382 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
383 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
387 fEntriesCH = new Int_t[Ntotal0];
388 for(Int_t k = 0; k < Ntotal0; k++){
394 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
395 ,fCalibraMode->GetNz(1)
396 ,fCalibraMode->GetNrphi(1)));
398 // Calcul the number of Xbins
399 Int_t Ntotal1 = CalculateTotalNumberOfBins(1);
400 // Create the 2D histo
405 fPHPlace = new Short_t[fTimeMax];
406 for (Int_t k = 0; k < fTimeMax; k++) {
409 fPHValue = new Float_t[fTimeMax];
410 for (Int_t k = 0; k < fTimeMax; k++) {
414 if (fLinearFitterOn) {
415 fLinearFitterArray.Expand(540);
416 fLinearFitterArray.SetName("ArrayLinearFitters");
417 fEntriesLinearFitter = new Int_t[540];
418 for(Int_t k = 0; k < 540; k++){
419 fEntriesLinearFitter[k] = 0;
421 if(fLinearFitterDebugOn) {
422 fLinearFitterHistoArray.Expand(540);
423 fLinearFitterHistoArray.SetName("ArrayHistos");
429 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
430 ,fCalibraMode->GetNz(2)
431 ,fCalibraMode->GetNrphi(2)));
433 // Calcul the number of Xbins
434 Int_t Ntotal2 = CalculateTotalNumberOfBins(2);
435 // Create the 2D histo
437 CreatePRF2d(Ntotal2);
445 //____________Functions for filling the histos in the code_____________________
447 //____________Offine tracking in the AliTRDtracker_____________________________
448 Bool_t AliTRDCalibraFillHisto::ResetTrack()
451 // For the offline tracking
452 // This function will be called in the function
453 // AliTRDtracker::FollowBackPropagation() at the beginning
454 // Reset the parameter to know we have a new TRD track
457 fDetectorAliTRDtrack = kFALSE;
462 //____________Offine tracking in the AliTRDtracker_____________________________
463 void AliTRDCalibraFillHisto::ResetfVariables()
466 // Reset values per tracklet
469 // Reset the list of clusters
470 fListClusters->Clear();
472 //Reset the tracklet parameters
473 for(Int_t k = 0; k < fTimeMax; k++){
482 //Reset good tracklet
483 fGoodTracklet = kTRUE;
485 // Reset the fPHValue
487 //Reset the fPHValue and fPHPlace
488 for (Int_t k = 0; k < fTimeMax; k++) {
494 // Reset the fAmpTotal where we put value
496 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
502 //____________Offline tracking in the AliTRDtracker____________________________
503 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
506 // For the offline tracking
507 // This function will be called in the function
508 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
510 // Fill the 2D histos or the vectors with the info of the clusters at
511 // the end of a detectors if the track is "good"
514 // Localisation of the detector
515 Int_t detector = cl->GetDetector();
517 // Fill the infos for the previous clusters if not the same
518 // detector anymore or if not the same track
519 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
520 (fDetectorPreviousTrack != -1)) {
524 // If the same track, then look if the previous detector is in
525 // the same plane, if yes: not a good track
526 if (fDetectorAliTRDtrack &&
527 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
531 // Fill only if the track doesn't touch a masked pad or doesn't
532 // appear in the middle (fGoodTrack)
533 if (fGoodTrack && fGoodTracklet) {
535 // drift velocity unables to cut bad tracklets
536 Bool_t pass = FindP1TrackPH();
540 FillTheInfoOfTheTrackCH();
545 FillTheInfoOfTheTrackPH();
548 if(pass && fPRF2dOn) HandlePRF();
555 } // Fill at the end the charge
557 // Calcul the position of the detector
558 if (detector != fDetectorPreviousTrack) {
559 LocalisationDetectorXbins(detector);
562 // Reset the detector
563 fDetectorPreviousTrack = detector;
564 fDetectorAliTRDtrack = kTRUE;
566 // Store the infos of the tracklets
567 AliTRDcluster *kcl = new AliTRDcluster(*cl);
568 fListClusters->Add((TObject *)kcl);
569 Int_t time = cl->GetLocalTimeBin();
570 fPar0[time] = t->GetY();
571 fPar1[time] = t->GetZ();
572 fPar2[time] = t->GetSnp();
573 fPar3[time] = t->GetTgl();
574 fPar4[time] = t->Get1Pt();
576 // Store the info bis of the tracklet
577 Int_t *rowcol = CalculateRowCol(cl);
578 CheckGoodTracklet(detector,rowcol);
579 Int_t group[2] = {0,0};
580 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
581 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
582 StoreInfoCHPH(cl,t,group);
587 //____________Online trackling in AliTRDtrigger________________________________
588 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
592 // This function will be called in the function AliTRDtrigger::TestTracklet
593 // before applying the pt cut on the tracklets
594 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
597 // Localisation of the Xbins involved
598 Int_t idect = trk->GetDetector();
599 LocalisationDetectorXbins(idect);
601 // Get the parameter object
602 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
604 AliInfo("Could not get calibDB");
611 // Row of the tracklet and position in the pad groups
612 Int_t *rowcol = new Int_t[2];
613 rowcol[0] = trk->GetRow();
614 Int_t group[3] = {-1,-1,-1};
616 // Eventuelle correction due to track angle in z direction
617 Float_t correction = 1.0;
618 if (fMcmCorrectAngle) {
619 Float_t z = trk->GetRowz();
620 Float_t r = trk->GetTime0();
621 correction = r / TMath::Sqrt((r*r+z*z));
624 // Boucle sur les clusters
625 // Condition on number of cluster: don't come from the middle of the detector
626 if (trk->GetNclusters() >= fNumberClusters) {
628 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
630 Float_t amp[3] = { 0.0, 0.0, 0.0 };
631 Int_t time = trk->GetClusterTime(icl);
632 rowcol[1] = trk->GetClusterCol(icl);
634 amp[0] = trk->GetClusterADC(icl)[0] * correction;
635 amp[1] = trk->GetClusterADC(icl)[1] * correction;
636 amp[2] = trk->GetClusterADC(icl)[2] * correction;
639 if ((amp[0] < 0.0) ||
645 // Col of cluster and position in the pad groups
647 group[0] = CalculateCalibrationGroup(0,rowcol);
648 fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
651 group[1] = CalculateCalibrationGroup(1,rowcol);
652 fPHPlace[time] = group[1];
653 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
655 if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
657 // See if we are not near a masked pad fGoodTracklet
658 CheckGoodTracklet(idect,rowcol);
660 // Fill PRF direct without tnp bins...only for monitoring...
661 if (fPRF2dOn && fGoodTracklet) {
663 if ((amp[0] > fThresholdClusterPRF2) &&
664 (amp[1] > fThresholdClusterPRF2) &&
665 (amp[2] > fThresholdClusterPRF2) &&
666 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
668 // Security of the denomiateur is 0
669 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
670 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
671 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
672 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
673 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
675 if (TMath::Abs(xcenter) < 0.5) {
676 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
677 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
678 // Fill only if it is in the drift region!
679 if (((Float_t) time / fSf) > 0.3) {
681 fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
682 fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
683 fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
686 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],xcenter,ycenter);
687 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],-(xcenter+1.0),yminus);
688 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],(1.0-xcenter),ymax);
690 }//in the drift region
692 }//denominateur security
693 }//cluster shape and thresholds
700 if (fCH2dOn) FillTheInfoOfTheTrackCH();
701 if (fPH2dOn) FillTheInfoOfTheTrackPH();
706 } // Condition on number of clusters
711 //_____________________________________________________________________________
712 Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
715 // Calculate the row and col number of the cluster
719 Int_t *rowcol = new Int_t[2];
723 // Get the parameter object
724 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
726 AliInfo("Could not get CommonParam");
730 // Localisation of the detector
731 Int_t detector = cl->GetDetector();
732 Int_t chamber = GetChamber(detector);
733 Int_t plane = GetPlane(detector);
735 // Localisation of the cluster
736 Double_t pos[3] = { 0.0, 0.0, 0.0 };
737 pos[0] = ((AliCluster *)cl)->GetX();
741 // Position of the cluster
742 AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
743 Int_t row = padplane->GetPadRowNumber(pos[2]);
744 //Do not take from here because it was corrected from ExB already....
745 //Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
746 //Double_t offsettilt = padplane->GetTiltOffset(offsetz);
747 //Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
748 //Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt);
749 Int_t col = cl->GetPad();
757 //_____________________________________________________________________________
758 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
761 // See if we are not near a masked pad
764 Int_t row = rowcol[0];
765 Int_t col = rowcol[1];
767 if (!IsPadOn(detector, col, row)) {
768 fGoodTracklet = kFALSE;
772 if (!IsPadOn(detector, col-1, row)) {
773 fGoodTracklet = kFALSE;
778 if (!IsPadOn(detector, col+1, row)) {
779 fGoodTracklet = kFALSE;
784 //_____________________________________________________________________________
785 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
788 // Look in the choosen database if the pad is On.
789 // If no the track will be "not good"
792 // Get the parameter object
793 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
795 AliInfo("Could not get calibDB");
799 if (!cal->IsChamberInstalled(detector) ||
800 cal->IsChamberMasked(detector) ||
801 cal->IsPadMasked(detector,col,row)) {
809 //_____________________________________________________________________________
810 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
813 // Calculate the calibration group number for i
816 // Row of the cluster and position in the pad groups
818 if (fCalibraMode->GetNnZ(i) != 0) {
819 posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
823 // Col of the cluster and position in the pad groups
825 if (fCalibraMode->GetNnRphi(i) != 0) {
826 posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
829 return posc*fCalibraMode->GetNfragZ(i)+posr;
832 //____________________________________________________________________________________
833 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
836 // Calculate the total number of calibration groups
840 fCalibraMode->ModePadCalibration(2,i);
841 fCalibraMode->ModePadFragmentation(0,2,0,i);
842 fCalibraMode->SetDetChamb2(i);
843 Ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
844 fCalibraMode->ModePadCalibration(0,i);
845 fCalibraMode->ModePadFragmentation(0,0,0,i);
846 fCalibraMode->SetDetChamb0(i);
847 Ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
848 AliInfo(Form("Total number of Xbins: %d",Ntotal));
852 //____________Set the pad calibration variables for the detector_______________
853 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
856 // For the detector calcul the first Xbins and set the number of row
857 // and col pads per calibration groups, the number of calibration
858 // groups in the detector.
861 // first Xbins of the detector
863 fCalibraMode->CalculXBins(detector,0);
866 fCalibraMode->CalculXBins(detector,1);
869 fCalibraMode->CalculXBins(detector,2);
872 // fragmentation of idect
873 for (Int_t i = 0; i < 3; i++) {
874 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
875 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
876 , (Int_t) GetChamber(detector)
877 , (Int_t) GetSector(detector),i);
883 //_____________________________________________________________________________
884 void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group)
887 // Store the infos in fAmpTotal, fPHPlace and fPHValue
890 // Charge in the cluster
891 Float_t q = TMath::Abs(cl->GetQ());
892 Int_t time = cl->GetLocalTimeBin();
894 // Correction due to the track angle
895 Float_t correction = 1.0;
896 Float_t normalisation = 6.67;
897 if ((q >0) && (t->GetNdedx() > 0)) {
898 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
901 // Fill the fAmpTotal with the charge
903 fAmpTotal[(Int_t) group[0]] += correction;
906 // Fill the fPHPlace and value
908 fPHPlace[time] = group[1];
909 fPHValue[time] = correction;
913 //_____________________________________________________________________
914 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDRawStream *rawStream)
917 // Event Processing loop - AliTRDRawStream
920 Bool_t withInput = kFALSE;
925 for(Int_t k = 0; k < 36; k++){
930 fDetectorPreviousTrack = -1;
933 while (rawStream->Next()) {
935 Int_t idetector = rawStream->GetDet(); // current detector
936 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
937 if(TMath::Mean(nbtimebin,phvalue)>0.0){
939 for(Int_t k = 0; k < nbtimebin; k++){
940 UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
947 fDetectorPreviousTrack = idetector;
948 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
949 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
950 //row[iTimeBin] = rawStream->GetRow(); // current row
951 //col[iTimeBin] = rawStream->GetCol(); // current col
952 Int_t *signal = rawStream->GetSignals(); // current ADC signal
954 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
956 for(Int_t itime = iTimeBin; itime < fin; itime++){
957 // should extract baseline here!
958 if(signal[n]>5.0) phvalue[itime] = signal[n];
964 if(fDetectorPreviousTrack != -1){
965 if(TMath::Mean(nbtimebin,phvalue)>0.0){
967 for(Int_t k = 0; k < nbtimebin; k++){
968 UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
979 //_____________________________________________________________________
980 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
983 // Event processing loop - AliRawReader
987 AliTRDRawStream rawStream(rawReader);
989 rawReader->Select("TRD");
991 return ProcessEventDAQ(&rawStream);
993 //_________________________________________________________________________
994 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(
996 eventHeaderStruct *event
998 eventHeaderStruct* /*event*/
1004 // process date event
1007 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1008 Bool_t result=ProcessEventDAQ(rawReader);
1012 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1017 //____________Online trackling in AliTRDtrigger________________________________
1018 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1022 // Fill a simple average pulse height
1025 // Localisation of the Xbins involved
1026 //LocalisationDetectorXbins(det);
1028 // Row and position in the pad groups
1030 //if (fCalibraMode->GetNnZ(1) != 0) {
1031 // posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1034 // Col of cluster and position in the pad groups
1036 //if (fCalibraMode->GetNnRphi(1) != 0) {
1037 // posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1040 //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1042 ((TProfile2D *)GetPH2d(nbtimebins,fSf,kTRUE))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1047 //____________Functions for plotting the 2D____________________________________
1049 //_____________________________________________________________________________
1050 void AliTRDCalibraFillHisto::Plot2d()
1053 // Plot the 2D histos
1057 TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
1059 fPH2d->Draw("LEGO");
1062 TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
1064 fCH2d->Draw("LEGO");
1067 TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
1069 fPRF2d->Draw("LEGO");
1073 //_____________________Reset____________________________________________________
1074 //_______________________________________________________________________________
1075 void AliTRDCalibraFillHisto::ResetLinearFitter()
1077 fLinearFitterArray.SetOwner();
1078 fLinearFitterArray.Clear();
1079 fLinearFitterHistoArray.SetOwner();
1080 fLinearFitterHistoArray.Clear();
1082 //____________Write_____________________________________________________
1083 //_____________________________________________________________________
1084 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1087 // Write infos to file
1091 if ( fDebugStreamer ) {
1092 delete fDebugStreamer;
1093 fDebugStreamer = 0x0;
1096 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1101 ,fNumberUsedPh[1]));
1103 TDirectory *backup = gDirectory;
1109 option = "recreate";
1111 TFile f(filename,option.Data());
1113 TStopwatch stopwatch;
1118 f.WriteTObject(fCH2d);
1122 name += fCalibraMode->GetNz(0);
1124 name += fCalibraMode->GetNrphi(0);
1125 TTree *treeCH2d = fCalibraVector->ConvertVectorCTTreeHisto(fCalibraVector->GetVectorCH(),fCalibraVector->GetPlaCH(),"treeCH2d",(const char *) name);
1126 f.WriteTObject(treeCH2d);
1131 f.WriteTObject(fPH2d);
1135 name += fCalibraMode->GetNz(1);
1137 name += fCalibraMode->GetNrphi(1);
1138 TTree *treePH2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPH(),fCalibraVector->GetPlaPH(),"treePH2d",(const char *) name);
1139 f.WriteTObject(treePH2d);
1144 f.WriteTObject(fPRF2d);
1148 name += fCalibraMode->GetNz(2);
1150 name += fCalibraMode->GetNrphi(2);
1153 TTree *treePRF2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPRF(),fCalibraVector->GetPlaPRF(),"treePRF2d",(const char *) name);
1154 f.WriteTObject(treePRF2d);
1157 if(fLinearFitterOn && fLinearFitterDebugOn){
1158 f.WriteTObject(&fLinearFitterHistoArray);
1163 if ( backup ) backup->cd();
1165 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1166 ,stopwatch.RealTime(),stopwatch.CpuTime()));
1168 //___________________________________________probe the histos__________________________________________________
1169 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1172 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1173 // debug mode with 2 for TH2I and 3 for TProfile2D
1174 // It gives a pointer to a Double_t[7] with the info following...
1175 // [0] : number of calibration groups with entries
1176 // [1] : minimal number of entries found
1177 // [2] : calibration group number of the min
1178 // [3] : maximal number of entries found
1179 // [4] : calibration group number of the max
1180 // [5] : mean number of entries found
1181 // [6] : mean relativ error
1184 Double_t *info = new Double_t[7];
1186 // Number of Xbins (detectors or groups of pads)
1187 Int_t nbins = h->GetNbinsY(); //number of calibration groups
1188 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
1191 Double_t nbwe = 0; //number of calibration groups with entries
1192 Double_t minentries = 0; //minimal number of entries found
1193 Double_t maxentries = 0; //maximal number of entries found
1194 Double_t placemin = 0; //calibration group number of the min
1195 Double_t placemax = -1; //calibration group number of the max
1196 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1197 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1199 Double_t counter = 0;
1202 TH1F *NbEntries = 0x0;//distribution of the number of entries
1203 TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1204 TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1206 // Beginning of the loop over the calibration groups
1207 for (Int_t idect = 0; idect < nbins; idect++) {
1209 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1210 projch->SetDirectory(0);
1212 // Number of entries for this calibration group
1213 Double_t nentries = 0.0;
1215 for (Int_t k = 0; k < nxbins; k++) {
1216 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1220 for (Int_t k = 0; k < nxbins; k++) {
1221 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1222 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1223 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1231 if((!((Bool_t)NbEntries)) && (nentries > 0)){
1232 NbEntries = new TH1F("Number of entries","Number of entries"
1233 ,100,(Int_t)nentries/2,nentries*2);
1234 NbEntries->SetDirectory(0);
1235 NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1237 NbEntriesPerGroup->SetDirectory(0);
1238 NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1239 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1240 NbEntriesPerSp->SetDirectory(0);
1243 if(nentries > 0) NbEntries->Fill(nentries);
1244 NbEntriesPerGroup->Fill(idect+0.5,nentries);
1245 NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1250 if(nentries > maxentries){
1251 maxentries = nentries;
1255 minentries = nentries;
1257 if(nentries < minentries){
1258 minentries = nentries;
1264 meanstats += nentries;
1266 }//calibration groups loop
1268 if(nbwe > 0) meanstats /= nbwe;
1269 if(counter > 0) meanrelativerror /= counter;
1271 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1272 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1273 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1274 AliInfo(Form("The mean number of entries is %f",meanstats));
1275 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1278 info[1] = minentries;
1280 info[3] = maxentries;
1282 info[5] = meanstats;
1283 info[6] = meanrelativerror;
1286 gStyle->SetPalette(1);
1287 gStyle->SetOptStat(1111);
1288 gStyle->SetPadBorderMode(0);
1289 gStyle->SetCanvasColor(10);
1290 gStyle->SetPadLeftMargin(0.13);
1291 gStyle->SetPadRightMargin(0.01);
1292 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1295 NbEntries->Draw("");
1297 NbEntriesPerSp->SetStats(0);
1298 NbEntriesPerSp->Draw("");
1299 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1301 NbEntriesPerGroup->SetStats(0);
1302 NbEntriesPerGroup->Draw("");
1308 //____________________________________________________________________________
1309 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1312 // Return a Int_t[4] with:
1313 // 0 Mean number of entries
1314 // 1 median of number of entries
1315 // 2 rms of number of entries
1316 // 3 number of group with entries
1319 Double_t *stat = new Double_t[4];
1322 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
1323 Double_t *weight = new Double_t[nbofgroups];
1324 Int_t *nonul = new Int_t[nbofgroups];
1326 for(Int_t k = 0; k < nbofgroups; k++){
1327 if(fEntriesCH[k] > 0) {
1329 nonul[(Int_t)stat[3]] = fEntriesCH[k];
1332 else weight[k] = 0.0;
1334 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
1335 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
1336 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1341 //____________________________________________________________________________
1342 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1345 // Return a Int_t[4] with:
1346 // 0 Mean number of entries
1347 // 1 median of number of entries
1348 // 2 rms of number of entries
1349 // 3 number of group with entries
1352 Double_t *stat = new Double_t[4];
1355 Int_t nbofgroups = 540;
1356 Double_t *weight = new Double_t[nbofgroups];
1357 Int_t *nonul = new Int_t[nbofgroups];
1359 for(Int_t k = 0; k < nbofgroups; k++){
1360 if(fEntriesLinearFitter[k] > 0) {
1362 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1365 else weight[k] = 0.0;
1367 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1368 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1369 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1374 //_____________________________________________________________________________
1375 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1378 // Should be between 0 and 6
1381 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1382 AliInfo("The number of groups must be between 0 and 6!");
1385 fNgroupprf = numberGroupsPRF;
1389 //_____________________________________________________________________________
1390 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1393 // Set the factor that will divide the deposited charge
1394 // to fit in the histo range [0,300]
1397 if (RelativeScale > 0.0) {
1398 fRelativeScale = RelativeScale;
1401 AliInfo("RelativeScale must be strict positif!");
1406 //_____________________________________________________________________________
1407 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1410 // Set the mode of calibration group in the z direction for the parameter i
1415 fCalibraMode->SetNz(i, Nz);
1418 AliInfo("You have to choose between 0 and 4");
1423 //_____________________________________________________________________________
1424 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1427 // Set the mode of calibration group in the rphi direction for the parameter i
1432 fCalibraMode->SetNrphi(i ,Nrphi);
1435 AliInfo("You have to choose between 0 and 6");
1439 //____________Protected Functions______________________________________________
1440 //____________Create the 2D histo to be filled online__________________________
1442 //_____________________________________________________________________________
1443 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1446 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1447 // If fNgroupprf is zero then no binning in tnp
1451 name += fCalibraMode->GetNz(2);
1453 name += fCalibraMode->GetNrphi(2);
1457 if(fNgroupprf != 0){
1459 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1460 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1461 fPRF2d->SetYTitle("Det/pad groups");
1462 fPRF2d->SetXTitle("Position x/W [pad width units]");
1463 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1464 fPRF2d->SetStats(0);
1467 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1468 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1469 fPRF2d->SetYTitle("Det/pad groups");
1470 fPRF2d->SetXTitle("Position x/W [pad width units]");
1471 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1472 fPRF2d->SetStats(0);
1477 //_____________________________________________________________________________
1478 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1481 // Create the 2D histos
1485 name += fCalibraMode->GetNz(1);
1487 name += fCalibraMode->GetNrphi(1);
1489 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1490 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1492 fPH2d->SetYTitle("Det/pad groups");
1493 fPH2d->SetXTitle("time [#mus]");
1494 fPH2d->SetZTitle("<PH> [a.u.]");
1499 //_____________________________________________________________________________
1500 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1503 // Create the 2D histos
1507 name += fCalibraMode->GetNz(0);
1509 name += fCalibraMode->GetNrphi(0);
1511 fCH2d = new TH2I("CH2d",(const Char_t *) name
1512 ,fNumberBinCharge,0,300,nn,0,nn);
1513 fCH2d->SetYTitle("Det/pad groups");
1514 fCH2d->SetXTitle("charge deposit [a.u]");
1515 fCH2d->SetZTitle("counts");
1520 //____________Offine tracking in the AliTRDtracker_____________________________
1521 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1524 // For the offline tracking or mcm tracklets
1525 // This function will be called in the functions UpdateHistogram...
1526 // to fill the info of a track for the relativ gain calibration
1529 Int_t nb = 0; // Nombre de zones traversees
1530 Int_t fd = -1; // Premiere zone non nulle
1531 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1534 // See if the track goes through different zones
1535 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1536 if (fAmpTotal[k] > 0.0) {
1537 totalcharge += fAmpTotal[k];
1550 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1552 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1553 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1556 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1560 if ((fAmpTotal[fd] > 0.0) &&
1561 (fAmpTotal[fd+1] > 0.0)) {
1562 // One of the two very big
1563 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1565 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1566 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1569 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1572 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1574 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1576 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1577 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1580 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1583 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1586 if (fCalibraMode->GetNfragZ(0) > 1) {
1587 if (fAmpTotal[fd] > 0.0) {
1588 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1589 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1590 // One of the two very big
1591 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1593 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1594 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1597 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1600 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1602 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1604 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1605 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1608 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1610 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1621 //____________Offine tracking in the AliTRDtracker_____________________________
1622 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1625 // For the offline tracking or mcm tracklets
1626 // This function will be called in the functions UpdateHistogram...
1627 // to fill the info of a track for the drift velocity calibration
1630 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1631 Int_t fd1 = -1; // Premiere zone non nulle
1632 Int_t fd2 = -1; // Deuxieme zone non nulle
1633 Int_t k1 = -1; // Debut de la premiere zone
1634 Int_t k2 = -1; // Debut de la seconde zone
1636 // See if the track goes through different zones
1637 for (Int_t k = 0; k < fTimeMax; k++) {
1638 if (fPHValue[k] > 0.0) {
1643 if (fPHPlace[k] != fd1) {
1649 if (fPHPlace[k] != fd2) {
1660 for (Int_t i = 0; i < fTimeMax; i++) {
1662 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1665 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1670 if ((fd1 == fd2+1) ||
1672 // One of the two fast all the think
1673 if (k2 > (k1+fDifference)) {
1674 //we choose to fill the fd1 with all the values
1676 for (Int_t i = 0; i < fTimeMax; i++) {
1678 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1681 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1685 if ((k2+fDifference) < fTimeMax) {
1686 //we choose to fill the fd2 with all the values
1688 for (Int_t i = 0; i < fTimeMax; i++) {
1690 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1693 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1698 // Two zones voisines sinon rien!
1699 if (fCalibraMode->GetNfragZ(1) > 1) {
1701 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1702 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1703 // One of the two fast all the think
1704 if (k2 > (k1+fDifference)) {
1705 //we choose to fill the fd1 with all the values
1707 for (Int_t i = 0; i < fTimeMax; i++) {
1709 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1712 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1716 if ((k2+fDifference) < fTimeMax) {
1717 //we choose to fill the fd2 with all the values
1719 for (Int_t i = 0; i < fTimeMax; i++) {
1721 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1724 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1730 // Two zones voisines sinon rien!
1732 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1733 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1734 // One of the two fast all the think
1735 if (k2 > (k1 + fDifference)) {
1736 //we choose to fill the fd1 with all the values
1738 for (Int_t i = 0; i < fTimeMax; i++) {
1740 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1743 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1747 if ((k2+fDifference) < fTimeMax) {
1748 //we choose to fill the fd2 with all the values
1750 for (Int_t i = 0; i < fTimeMax; i++) {
1752 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1755 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1766 //____________Offine tracking in the AliTRDtracker_____________________________
1767 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
1770 // For the offline tracking
1771 // This function will be called in the functions UpdateHistogram...
1772 // to fill the find the parameter P1 of a track for the drift velocity calibration
1775 // Get the parameter object
1776 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1778 AliInfo("Could not get CommonParam");
1782 //Number of points: if less than 3 return kFALSE
1783 Int_t Npoints = fListClusters->GetEntriesFast();
1784 if(Npoints <= 2) return kFALSE;
1787 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
1788 Double_t snp = 0.0; // sin angle in the plan yx track
1789 Double_t y = 0.0; // y clusters in the middle of the chamber
1790 Double_t z = 0.0; // z cluster in the middle of the chamber
1791 Double_t dydt = 0.0; // dydt tracklet after straight line fit
1792 Double_t tnp = 0.0; // tan angle in the plan xy track
1793 Double_t tgl = 0.0; // dz/dl and not dz/dx!
1794 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
1795 Double_t pointError = 0.0; // error after straight line fit
1796 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
1797 Int_t snpright = 1; // if we took in the middle snp
1798 Int_t crossrow = 0; // if it crosses a pad row
1799 Double_t tiltingangle = 0; // tiltingangle of the pad
1800 Float_t dzdx = 0; // dz/dx now from dz/dl
1801 Int_t nbli = 0; // number linear fitter points
1802 AliTRDpadPlane *padplane = parCom->GetPadPlane(GetPlane(detector),GetChamber(detector));
1804 linearFitterTracklet.StoreData(kFALSE);
1805 linearFitterTracklet.ClearPoints();
1807 //if more than one row
1808 Int_t rowp = -1; // if it crosses a pad row
1811 tiltingangle = padplane->GetTiltingAngle();
1812 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
1815 for(Int_t k = 0; k < Npoints; k++){
1817 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1818 Double_t ycluster = cl->GetY();
1819 Int_t time = cl->GetLocalTimeBin();
1820 Double_t timeis = time/fSf;
1821 //See if cross two pad rows
1822 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
1823 if(k==0) rowp = row;
1824 if(row != rowp) crossrow = 1;
1825 //Take in the middle of the chamber
1826 if(time > (Int_t) 10) {
1832 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1835 if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() <= 10) snpright = 0;
1836 if(nbli <= 2) return kFALSE;
1838 // Do the straight line fit now
1840 linearFitterTracklet.Eval();
1841 linearFitterTracklet.GetParameters(pars);
1842 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
1843 errorpar = linearFitterTracklet.GetParError(1)*pointError;
1846 if( TMath::Abs(snp) < 1.){
1847 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1849 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1851 if(fDebugLevel > 0){
1852 if ( !fDebugStreamer ) {
1854 TDirectory *backup = gDirectory;
1855 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
1856 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1859 (* fDebugStreamer) << "VDRIFT0"<<
1860 "Npoints="<<Npoints<<
1864 (* fDebugStreamer) << "VDRIFT"<<
1865 "snpright="<<snpright<<
1866 "Npoints="<<Npoints<<
1868 "detector="<<detector<<
1877 "crossrow="<<crossrow<<
1878 "errorpar="<<errorpar<<
1879 "pointError="<<pointError<<
1884 if(Npoints < fNumberClusters) return kFALSE;
1885 if(snpright == 0) return kFALSE;
1886 if(pointError >= 0.1) return kFALSE;
1887 if(crossrow == 1) return kFALSE;
1889 if(fLinearFitterOn){
1890 //Add to the linear fitter of the detector
1891 if( TMath::Abs(snp) < 1.){
1892 Double_t x = tnp-dzdx*tnt;
1893 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
1894 if(fLinearFitterDebugOn) {
1895 ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars[1]);
1897 fEntriesLinearFitter[detector]++;
1900 //AliInfo("End of FindP1TrackPH with success!")
1904 //____________Offine tracking in the AliTRDtracker_____________________________
1905 Bool_t AliTRDCalibraFillHisto::HandlePRF()
1908 // For the offline tracking
1909 // Fit the tracklet with a line and take the position as reference for the PRF
1912 // Get the parameter object
1913 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1915 AliInfo("Could not get CommonParam");
1920 Int_t Npoints = fListClusters->GetEntriesFast(); // number of total points
1921 Int_t Nb3pc = 0; // number of three pads clusters used for fit
1922 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
1925 // To see the difference due to the fit
1926 Double_t *padPositions;
1927 padPositions = new Double_t[Npoints];
1928 for(Int_t k = 0; k < Npoints; k++){
1929 padPositions[k] = 0.0;
1933 //Find the position by a fit
1934 TLinearFitter fitter(2,"pol1");
1935 fitter.StoreData(kFALSE);
1936 fitter.ClearPoints();
1937 for(Int_t k = 0; k < Npoints; k++){
1939 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1940 Short_t *signals = cl->GetSignals();
1941 Double_t time = cl->GetLocalTimeBin();
1942 //Calculate x if possible
1943 Float_t xcenter = 0.0;
1944 Bool_t echec = kTRUE;
1945 if((time<=7) || (time>=21)) continue;
1946 // Center 3 balanced: position with the center of the pad
1947 if ((((Float_t) signals[3]) > 0.0) &&
1948 (((Float_t) signals[2]) > 0.0) &&
1949 (((Float_t) signals[4]) > 0.0)) {
1951 // Security if the denomiateur is 0
1952 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1953 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1954 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1955 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1956 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1962 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1964 //if no echec: calculate with the position of the pad
1965 // Position of the cluster
1966 Double_t padPosition = xcenter + cl->GetPad();
1967 padPositions[k] = padPosition;
1969 fitter.AddPoint(&time, padPosition,1);
1972 //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
1973 if(Nb3pc < 3) return kFALSE;
1976 fitter.GetParameters(line);
1977 Float_t pointError = -1.0;
1978 pointError = TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
1982 for(Int_t k = 0; k < Npoints; k++){
1984 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1985 Short_t *signals = cl->GetSignals(); // signal
1986 Double_t time = cl->GetLocalTimeBin(); // time bin
1987 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1988 Float_t padPos = cl->GetPad(); // middle pad
1989 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1990 Float_t ycenter = 0.0; // relative center charge
1991 Float_t ymin = 0.0; // relative left charge
1992 Float_t ymax = 0.0; // relative right charge
1993 Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx
1994 Double_t pt = fPar4[(Int_t)time]; // pt
1995 Float_t dzdx = 0.0; // dzdx
1998 //Requiere simply two pads clusters at least
1999 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2000 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2001 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2002 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2003 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2004 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2008 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2009 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2010 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponidng group
2011 Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan
2012 Float_t xcl = cl->GetY(); // y cluster
2013 Float_t qcl = cl->GetQ(); // charge cluster
2014 Int_t plane = GetPlane(detector); // plane
2015 Int_t chamber = GetChamber(detector); // chamber
2016 Double_t xdiff = dpad; // reconstructed position constant
2017 Double_t x = dpad; // reconstructed position moved
2018 Float_t Ep = pointError; // error of fit
2019 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2020 Float_t signal3 = (Float_t)signals[3]; // signal
2021 Float_t signal2 = (Float_t)signals[2]; // signal
2022 Float_t signal4 = (Float_t)signals[4]; // signal
2023 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2025 if(TMath::Abs(snp) < 1.0){
2026 tnp = snp / (TMath::Sqrt(1-snp*snp));
2027 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2031 if(fDebugLevel > 0){
2032 if ( !fDebugStreamer ) {
2034 TDirectory *backup = gDirectory;
2035 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
2036 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2039 (* fDebugStreamer) << "PRF0"<<
2040 "caligroup="<<caligroup<<
2041 "detector="<<detector<<
2043 "chamber="<<chamber<<
2044 "Npoints="<<Npoints<<
2053 "padPosition="<<padPositions[k]<<
2054 "padPosTracklet="<<padPosTracklet<<
2056 "ycenter="<<ycenter<<
2061 "signal1="<<signal1<<
2062 "signal2="<<signal2<<
2063 "signal3="<<signal3<<
2064 "signal4="<<signal4<<
2065 "signal5="<<signal5<<
2070 Float_t y = ycenter;
2071 (* fDebugStreamer) << "PRFALL"<<
2072 "caligroup="<<caligroup<<
2073 "detector="<<detector<<
2075 "chamber="<<chamber<<
2076 "Npoints="<<Npoints<<
2086 "padPosition="<<padPositions[k]<<
2087 "padPosTracklet="<<padPosTracklet<<
2092 "signal1="<<signal1<<
2093 "signal2="<<signal2<<
2094 "signal3="<<signal3<<
2095 "signal4="<<signal4<<
2096 "signal5="<<signal5<<
2102 (* fDebugStreamer) << "PRFALL"<<
2103 "caligroup="<<caligroup<<
2104 "detector="<<detector<<
2106 "chamber="<<chamber<<
2107 "Npoints="<<Npoints<<
2117 "padPosition="<<padPositions[k]<<
2118 "padPosTracklet="<<padPosTracklet<<
2123 "signal1="<<signal1<<
2124 "signal2="<<signal2<<
2125 "signal3="<<signal3<<
2126 "signal4="<<signal4<<
2127 "signal5="<<signal5<<
2134 (* fDebugStreamer) << "PRFALL"<<
2135 "caligroup="<<caligroup<<
2136 "detector="<<detector<<
2138 "chamber="<<chamber<<
2139 "Npoints="<<Npoints<<
2149 "padPosition="<<padPositions[k]<<
2150 "padPosTracklet="<<padPosTracklet<<
2155 "signal1="<<signal1<<
2156 "signal2="<<signal2<<
2157 "signal3="<<signal3<<
2158 "signal4="<<signal4<<
2159 "signal5="<<signal5<<
2166 if(Npoints < fNumberClusters) continue;
2167 if(Nb3pc <= 5) continue;
2168 if((time >= 21) || (time < 7)) continue;
2169 if(TMath::Abs(snp) >= 1.0) continue;
2170 if(qcl < 80) continue;
2172 Bool_t echec = kFALSE;
2173 Double_t shift = 0.0;
2174 //Calculate the shift in x coresponding to this tnp
2175 if(fNgroupprf != 0.0){
2176 shift = -3.0*(fNgroupprf-1)-1.5;
2177 Double_t limithigh = -0.2*(fNgroupprf-1);
2178 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2180 while(tnp > limithigh){
2186 if (fHisto2d && !echec) {
2187 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2189 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2190 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2193 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2194 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2197 //Not equivalent anymore here!
2198 if (fVector2d && !echec) {
2199 fCalibraVector->UpdateVectorPRF(caligroup,shift+dpad,ycenter);
2201 fCalibraVector->UpdateVectorPRF(caligroup,shift-(dpad+1.0),ymin);
2202 fCalibraVector->UpdateVectorPRF(caligroup,shift+(dpad+1.0),ymin);
2205 fCalibraVector->UpdateVectorPRF(caligroup,shift+1.0-dpad,ymax);
2206 fCalibraVector->UpdateVectorPRF(caligroup,shift-1.0+dpad,ymax);
2214 //____________Some basic geometry function_____________________________________
2216 //_____________________________________________________________________________
2217 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2220 // Reconstruct the plane number from the detector number
2223 return ((Int_t) (d % 6));
2227 //_____________________________________________________________________________
2228 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2231 // Reconstruct the chamber number from the detector number
2235 return ((Int_t) (d % 30) / fgkNplan);
2239 //_____________________________________________________________________________
2240 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2243 // Reconstruct the sector number from the detector number
2247 return ((Int_t) (d / fg));
2250 //_____________________________________________________________________
2251 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency, Bool_t force)
2254 // return pointer to fPH2d TProfile2D
2255 // if force is true create a new TProfile2D if it doesn't exist allready
2257 if ( !force || fPH2d )
2261 fTimeMax = nbtimebin;
2262 fSf = samplefrequency;
2265 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
2266 ,fCalibraMode->GetNz(1)
2267 ,fCalibraMode->GetNrphi(1)));
2269 // Calcul the number of Xbins
2271 fCalibraMode->ModePadCalibration(2,1);
2272 fCalibraMode->ModePadFragmentation(0,2,0,1);
2273 fCalibraMode->SetDetChamb2(1);
2274 Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
2275 fCalibraMode->ModePadCalibration(0,1);
2276 fCalibraMode->ModePadFragmentation(0,0,0,1);
2277 fCalibraMode->SetDetChamb0(1);
2278 Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
2279 AliInfo(Form("Total number of Xbins: %d",Ntotal1));
2281 CreatePH2d(Ntotal1);
2288 //_____________________________________________________________________
2289 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2292 // return pointer to TLinearFitter Calibration
2293 // if force is true create a new TLinearFitter if it doesn't exist allready
2295 if ( !force || fLinearFitterArray.UncheckedAt(detector) )
2296 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2298 // if we are forced and TLinearFitter doesn't yet exist create it
2300 // new TLinearFitter
2301 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2302 fLinearFitterArray.AddAt(linearfitter,detector);
2303 return linearfitter;
2305 //_____________________________________________________________________
2306 TH2F* AliTRDCalibraFillHisto::GetLinearFitterHisto(Int_t detector, Bool_t force)
2309 // return pointer to TH2F histo
2310 // if force is true create a new histo if it doesn't exist allready
2312 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
2313 return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
2315 // if we are forced and TLinearFitter doesn't yes exist create it
2318 TString name("LFDV");
2321 TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
2324 lfdv->SetXTitle("tan(phi_{track})");
2325 lfdv->SetYTitle("dy/dt");
2326 lfdv->SetZTitle("Number of clusters");
2328 lfdv->SetDirectory(0);
2330 fLinearFitterHistoArray.AddAt(lfdv,detector);
2333 //_____________________________________________________________________
2334 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2337 // FillCH2d: Marian style
2340 //skip simply the value out of range
2341 if((y>=300.0) || (y<0.0)) return;
2343 //Calcul the y place
2344 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2345 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2348 fCH2d->GetArray()[place]++;