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 //////////////////////////////////////////////////////////////////////////////////////
36 #include <TProfile2D.h>
41 #include <TObjArray.h>
46 #include <TStopwatch.h>
48 #include <TDirectory.h>
49 #include <TTreeStream.h>
54 #include "AliTRDCalibraFillHisto.h"
55 #include "AliTRDCalibraMode.h"
56 #include "AliTRDCalibraVector.h"
57 #include "AliTRDCalibraVdriftLinearFit.h"
58 #include "AliTRDcalibDB.h"
59 #include "AliTRDCommonParam.h"
60 #include "AliTRDmcmTracklet.h"
61 #include "AliTRDpadPlane.h"
62 #include "AliTRDcluster.h"
63 #include "AliTRDtrack.h"
64 #include "AliTRDrawStreamTB.h"
65 #include "AliRawReader.h"
66 #include "AliRawReaderDate.h"
67 #include "AliTRDgeometry.h"
68 #include "./Cal/AliTRDCalROC.h"
69 #include "./Cal/AliTRDCalDet.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) {
117 //______________________________________________________________________________________
118 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
122 ,fMcmTracking(kFALSE)
123 ,fMcmCorrectAngle(kFALSE)
129 ,fLinearFitterOn(kFALSE)
130 ,fLinearFitterDebugOn(kFALSE)
132 ,fThresholdClusterPRF2(15.0)
133 ,fCalibraMode(new AliTRDCalibraMode())
136 ,fDetectorAliTRDtrack(kFALSE)
137 ,fDetectorPreviousTrack(-1)
146 ,fNumberBinCharge(100)
149 ,fListClusters(new TObjArray())
158 ,fGoodTracklet(kTRUE)
161 ,fEntriesLinearFitter(0x0)
166 ,fLinearFitterArray(540)
167 ,fLinearVdriftFit(0x0)
174 // Default constructor
178 // Init some default values
181 fNumberUsedCh[0] = 0;
182 fNumberUsedCh[1] = 0;
183 fNumberUsedPh[0] = 0;
184 fNumberUsedPh[1] = 0;
186 fGeo = new AliTRDgeometry();
190 //______________________________________________________________________________________
191 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
194 ,fMITracking(c.fMITracking)
195 ,fMcmTracking(c.fMcmTracking)
196 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
199 ,fPRF2dOn(c.fPRF2dOn)
200 ,fHisto2d(c.fHisto2d)
201 ,fVector2d(c.fVector2d)
202 ,fLinearFitterOn(c.fLinearFitterOn)
203 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
204 ,fRelativeScale(c.fRelativeScale)
205 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
210 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
211 ,fMCMPrevious(c.fMCMPrevious)
212 ,fROBPrevious(c.fROBPrevious)
213 ,fNumberClusters(c.fNumberClusters)
214 ,fProcent(c.fProcent)
215 ,fDifference(c.fDifference)
216 ,fNumberTrack(c.fNumberTrack)
217 ,fTimeMax(c.fTimeMax)
219 ,fNumberBinCharge(c.fNumberBinCharge)
220 ,fNumberBinPRF(c.fNumberBinPRF)
221 ,fNgroupprf(c.fNgroupprf)
222 ,fListClusters(new TObjArray())
231 ,fGoodTracklet(c.fGoodTracklet)
232 ,fGoodTrack(c.fGoodTrack)
234 ,fEntriesLinearFitter(0x0)
239 ,fLinearFitterArray(540)
240 ,fLinearVdriftFit(0x0)
249 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
250 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
252 fPH2d = new TProfile2D(*c.fPH2d);
253 fPH2d->SetDirectory(0);
256 fPRF2d = new TProfile2D(*c.fPRF2d);
257 fPRF2d->SetDirectory(0);
260 fCH2d = new TH2I(*c.fCH2d);
261 fCH2d->SetDirectory(0);
263 if(c.fLinearVdriftFit){
264 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
267 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
268 if(c.fCalDetT0) fCalDetT0 = new AliTRDCalDet(*c.fCalDetT0);
269 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
270 if(c.fCalROCT0) fCalROCT0 = new AliTRDCalROC(*c.fCalROCT0);
275 fGeo = new AliTRDgeometry();
278 //____________________________________________________________________________________
279 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
282 // AliTRDCalibraFillHisto destructor
286 if ( fDebugStreamer ) delete fDebugStreamer;
288 if ( fCalDetGain ) delete fCalDetGain;
289 if ( fCalDetT0 ) delete fCalDetT0;
290 if ( fCalROCGain ) delete fCalROCGain;
291 if ( fCalROCT0 ) delete fCalROCT0;
299 //_____________________________________________________________________________
300 void AliTRDCalibraFillHisto::Destroy()
313 //_____________________________________________________________________________
314 void AliTRDCalibraFillHisto::ClearHistos()
335 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
336 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
339 // For the offline tracking
340 // This function will be called in the function AliReconstruction::Run()
341 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
346 //Init the tracklet parameters
347 fPar0 = new Double_t[fTimeMax];
348 fPar1 = new Double_t[fTimeMax];
349 fPar2 = new Double_t[fTimeMax];
350 fPar3 = new Double_t[fTimeMax];
351 fPar4 = new Double_t[fTimeMax];
353 for(Int_t k = 0; k < fTimeMax; k++){
363 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
364 Bool_t AliTRDCalibraFillHisto::Init2Dhistostrack()
367 // For the offline tracking
368 // This function will be called in the function AliReconstruction::Run()
369 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
374 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
376 AliInfo("Could not get calibDB");
379 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
381 AliInfo("Could not get CommonParam");
386 fTimeMax = cal->GetNumberOfTimeBins();
387 fSf = parCom->GetSamplingFrequency();
390 //calib object from database used for reconstruction
391 if(fCalDetGain) delete fCalDetGain;
392 fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
393 if(fCalDetT0) delete fCalDetT0;
394 fCalDetT0 = new AliTRDCalDet(*(cal->GetT0Det()));
396 // Calcul Xbins Chambd0, Chamb2
397 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
398 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
399 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
401 // If vector method On initialised all the stuff
403 fCalibraVector = new AliTRDCalibraVector();
404 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
405 fCalibraVector->SetTimeMax(fTimeMax);
406 if(fNgroupprf != 0) {
407 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
408 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
411 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
412 fCalibraVector->SetPRFRange(1.5);
414 for(Int_t k = 0; k < 3; k++){
415 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
416 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
418 TString namech("Nz");
419 namech += fCalibraMode->GetNz(0);
421 namech += fCalibraMode->GetNrphi(0);
422 fCalibraVector->SetNameCH((const char* ) namech);
423 TString nameph("Nz");
424 nameph += fCalibraMode->GetNz(1);
426 nameph += fCalibraMode->GetNrphi(1);
427 fCalibraVector->SetNamePH((const char* ) nameph);
428 TString nameprf("Nz");
429 nameprf += fCalibraMode->GetNz(2);
431 nameprf += fCalibraMode->GetNrphi(2);
433 nameprf += fNgroupprf;
434 fCalibraVector->SetNamePRF((const char* ) nameprf);
437 // Create the 2D histos corresponding to the pad groupCalibration mode
440 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
441 ,fCalibraMode->GetNz(0)
442 ,fCalibraMode->GetNrphi(0)));
444 // Create the 2D histo
449 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
450 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
454 fEntriesCH = new Int_t[ntotal0];
455 for(Int_t k = 0; k < ntotal0; k++){
462 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
463 ,fCalibraMode->GetNz(1)
464 ,fCalibraMode->GetNrphi(1)));
466 // Create the 2D histo
471 fPHPlace = new Short_t[fTimeMax];
472 for (Int_t k = 0; k < fTimeMax; k++) {
475 fPHValue = new Float_t[fTimeMax];
476 for (Int_t k = 0; k < fTimeMax; k++) {
480 if (fLinearFitterOn) {
481 //fLinearFitterArray.Expand(540);
482 fLinearFitterArray.SetName("ArrayLinearFitters");
483 fEntriesLinearFitter = new Int_t[540];
484 for(Int_t k = 0; k < 540; k++){
485 fEntriesLinearFitter[k] = 0;
487 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
492 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
493 ,fCalibraMode->GetNz(2)
494 ,fCalibraMode->GetNrphi(2)));
495 // Create the 2D histo
497 CreatePRF2d(ntotal2);
505 //____________Functions for filling the histos in the code_____________________
507 //____________Offine tracking in the AliTRDtracker_____________________________
508 Bool_t AliTRDCalibraFillHisto::ResetTrack()
511 // For the offline tracking
512 // This function will be called in the function
513 // AliTRDtracker::FollowBackPropagation() at the beginning
514 // Reset the parameter to know we have a new TRD track
517 fDetectorAliTRDtrack = kFALSE;
518 //fGoodTrack = kTRUE;
522 //____________Offine tracking in the AliTRDtracker_____________________________
523 void AliTRDCalibraFillHisto::ResetfVariables()
526 // Reset values per tracklet
529 // Reset the list of clusters
530 fListClusters->Clear();
532 //Reset the tracklet parameters
533 for(Int_t k = 0; k < fTimeMax; k++){
541 ResetfVariablestrack();
544 //____________Offine tracking in the AliTRDtracker_____________________________
545 void AliTRDCalibraFillHisto::ResetfVariablestrack()
548 // Reset values per tracklet
551 //Reset good tracklet
552 fGoodTracklet = kTRUE;
554 // Reset the fPHValue
556 //Reset the fPHValue and fPHPlace
557 for (Int_t k = 0; k < fTimeMax; k++) {
563 // Reset the fAmpTotal where we put value
565 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
570 //____________Offline tracking in the AliTRDtracker____________________________
571 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
574 // For the offline tracking
575 // This function will be called in the function
576 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
578 // Fill the 2D histos or the vectors with the info of the clusters at
579 // the end of a detectors if the track is "good"
584 AliTRDcluster *cl = 0x0;
588 // reset if good track
592 // loop over the clusters
593 while((cl = t->GetCluster(index1))){
595 // Localisation of the detector
596 Int_t detector = cl->GetDetector();
599 // Fill the infos for the previous clusters if not the same
600 // detector anymore but this time it should be the same track
601 if ((detector != fDetectorPreviousTrack) &&
602 (index0 != index1)) {
606 //printf("detector %d, fPreviousdetector %d, plane %d, planeprevious %d, index0 %d, index1 %d la\n",detector,fDetectorPreviousTrack,GetPlane(detector),GetPlane(fDetectorPreviousTrack),index0,index1);
608 //If the same track, then look if the previous detector is in
609 //the same plane, if yes: not a good track
611 //if (fDetectorAliTRDtrack &&
612 // (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
614 if ((GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) {
618 // Fill only if the track doesn't touch a masked pad or doesn't
619 // appear in the middle (fGoodTrack)
620 if (fGoodTrack && fGoodTracklet) {
622 // drift velocity unables to cut bad tracklets
623 Bool_t pass = FindP1TrackPHtrack(t,index0,index1);
627 FillTheInfoOfTheTrackCH();
632 FillTheInfoOfTheTrackPH();
635 if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1);
641 ResetfVariablestrack();
644 } // Fill at the end the charge
646 // Calcul the position of the detector and take the calib objects
647 if (detector != fDetectorPreviousTrack) {
649 //Localise the detector
650 LocalisationDetectorXbins(detector);
653 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
655 AliInfo("Could not get calibDB");
660 if( fCalROCGain ) delete fCalROCGain;
661 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
662 if( fCalROCT0 ) delete fCalROCT0;
663 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector)));
667 // Reset the detectbjobsor
668 fDetectorPreviousTrack = detector;
670 // Store the info bis of the tracklet
671 Int_t *rowcol = CalculateRowCol(cl);
672 CheckGoodTracklet(detector,rowcol);
673 Int_t group[2] = {0,0};
674 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
675 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
676 StoreInfoCHPHtrack(cl,t,index1,group,rowcol);
680 } // while on clusters
682 // Fill the last plane
683 if( index0 != index1 ){
685 //printf("fPreviousdetector %d, planeprevious %d, index0 %d, index1 %d li\n",fDetectorPreviousTrack,GetPlane(fDetectorPreviousTrack),index0,index1);
689 if (fGoodTrack && fGoodTracklet) {
691 // drift velocity unables to cut bad tracklets
692 Bool_t pass = FindP1TrackPHtrack(t,index0,index1);
696 FillTheInfoOfTheTrackCH();
701 FillTheInfoOfTheTrackPH();
704 if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1);
711 ResetfVariablestrack();
716 //____________Offline tracking in the AliTRDtracker____________________________
717 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
720 // For the offline tracking
721 // This function will be called in the function
722 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
724 // Fill the 2D histos or the vectors with the info of the clusters at
725 // the end of a detectors if the track is "good"
728 // Localisation of the detector
729 Int_t detector = cl->GetDetector();
731 // Fill the infos for the previous clusters if not the same
732 // detector anymore or if not the same track
733 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
734 (fDetectorPreviousTrack != -1)) {
738 // If the same track, then look if the previous detector is in
739 // the same plane, if yes: not a good track
741 if (fDetectorAliTRDtrack &&
742 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
744 //if (fDetectorAliTRDtrack &&
745 // (GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) {
749 // Fill only if the track doesn't touch a masked pad or doesn't
750 // appear in the middle (fGoodTrack)
751 if (fGoodTrack && fGoodTracklet) {
753 // drift velocity unables to cut bad tracklets
754 Bool_t pass = FindP1TrackPH();
758 FillTheInfoOfTheTrackCH();
763 FillTheInfoOfTheTrackPH();
766 if(pass && fPRF2dOn) HandlePRF();
772 if(!fDetectorAliTRDtrack) fGoodTrack = kTRUE;
774 } // Fill at the end the charge
776 // Calcul the position of the detector and take the calib objects
777 if (detector != fDetectorPreviousTrack) {
778 //Localise the detector
779 LocalisationDetectorXbins(detector);
782 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
784 AliInfo("Could not get calibDB");
789 if( fCalROCGain ) delete fCalROCGain;
790 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
791 if( fCalROCT0 ) delete fCalROCT0;
792 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector)));
795 // Reset the detector
796 fDetectorPreviousTrack = detector;
797 fDetectorAliTRDtrack = kTRUE;
799 // Store the infos of the tracklets
800 AliTRDcluster *kcl = new AliTRDcluster(*cl);
801 fListClusters->Add((TObject *)kcl);
802 Int_t time = cl->GetLocalTimeBin();
803 fPar0[time] = t->GetY();
804 fPar1[time] = t->GetZ();
805 fPar2[time] = t->GetSnp();
806 fPar3[time] = t->GetTgl();
807 fPar4[time] = t->GetSigned1Pt();
809 // Store the info bis of the tracklet
810 Int_t *rowcol = CalculateRowCol(cl);
811 CheckGoodTracklet(detector,rowcol);
812 Int_t group[2] = {0,0};
813 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
814 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
815 StoreInfoCHPH(cl,t,group,rowcol);
820 //____________Online trackling in AliTRDtrigger________________________________
821 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
825 // This function will be called in the function AliTRDtrigger::TestTracklet
826 // before applying the pt cut on the tracklets
827 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
830 // Localisation of the Xbins involved
831 Int_t idect = trk->GetDetector();
832 fDetectorPreviousTrack = idect;
833 LocalisationDetectorXbins(idect);
835 // Get the parameter object
836 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
838 AliInfo("Could not get calibDB");
846 if( fCalROCGain ) delete fCalROCGain;
847 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(idect)));
848 if( fCalROCT0 ) delete fCalROCT0;
849 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(idect)));
851 // Row of the tracklet and position in the pad groups
852 Int_t *rowcol = new Int_t[2];
853 rowcol[0] = trk->GetRow();
854 Int_t group[3] = {-1,-1,-1};
856 // Eventuelle correction due to track angle in z direction
857 Float_t correction = 1.0;
858 if (fMcmCorrectAngle) {
859 Float_t z = trk->GetRowz();
860 Float_t r = trk->GetTime0();
861 correction = r / TMath::Sqrt((r*r+z*z));
864 // Boucle sur les clusters
865 // Condition on number of cluster: don't come from the middle of the detector
866 if (trk->GetNclusters() >= fNumberClusters) {
868 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
870 Float_t amp[3] = { 0.0, 0.0, 0.0 };
871 Int_t time = trk->GetClusterTime(icl);
872 rowcol[1] = trk->GetClusterCol(icl);
874 amp[0] = trk->GetClusterADC(icl)[0] * correction;
875 amp[1] = trk->GetClusterADC(icl)[1] * correction;
876 amp[2] = trk->GetClusterADC(icl)[2] * correction;
879 if ((amp[0] < 0.0) ||
885 // Col of cluster and position in the pad groups
887 group[0] = CalculateCalibrationGroup(0,rowcol);
888 fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
891 group[1] = CalculateCalibrationGroup(1,rowcol);
892 fPHPlace[time] = group[1];
893 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
895 if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
897 // See if we are not near a masked pad fGoodTracklet
898 CheckGoodTracklet(idect,rowcol);
900 // Fill PRF direct without tnp bins...only for monitoring...
901 if (fPRF2dOn && fGoodTracklet) {
903 if ((amp[0] > fThresholdClusterPRF2) &&
904 (amp[1] > fThresholdClusterPRF2) &&
905 (amp[2] > fThresholdClusterPRF2) &&
906 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
908 // Security of the denomiateur is 0
909 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
910 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
911 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
912 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
913 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
915 if (TMath::Abs(xcenter) < 0.5) {
916 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
917 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
918 // Fill only if it is in the drift region!
919 if (((Float_t) time / fSf) > 0.3) {
921 fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
922 fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
923 fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
926 fCalibraVector->UpdateVectorPRF(idect,group[2],xcenter,ycenter);
927 fCalibraVector->UpdateVectorPRF(idect,group[2],-(xcenter+1.0),yminus);
928 fCalibraVector->UpdateVectorPRF(idect,group[2],(1.0-xcenter),ymax);
930 }//in the drift region
932 }//denominateur security
933 }//cluster shape and thresholds
940 if (fCH2dOn) FillTheInfoOfTheTrackCH();
941 if (fPH2dOn) FillTheInfoOfTheTrackPH();
946 } // Condition on number of clusters
951 //_____________________________________________________________________________
952 Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
955 // Calculate the row and col number of the cluster
959 Int_t *rowcol = new Int_t[2];
963 // Localisation of the detector
964 Int_t detector = cl->GetDetector();
965 Int_t chamber = GetChamber(detector);
966 Int_t plane = GetPlane(detector);
968 // Localisation of the cluster
969 Double_t pos[3] = { 0.0, 0.0, 0.0 };
970 pos[0] = ((AliCluster *)cl)->GetX();
974 // Position of the cluster
975 AliTRDpadPlane *padplane = fGeo->GetPadPlane(plane,chamber);
976 Int_t row = padplane->GetPadRowNumber(pos[2]);
977 //Do not take from here because it was corrected from ExB already....
978 //Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
979 //Double_t offsettilt = padplane->GetTiltOffset(offsetz);
980 //Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
981 //Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt);
982 Int_t col = cl->GetPadCol();
990 //_____________________________________________________________________________
991 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
994 // See if we are not near a masked pad
997 Int_t row = rowcol[0];
998 Int_t col = rowcol[1];
1000 if (!IsPadOn(detector, col, row)) {
1001 fGoodTracklet = kFALSE;
1005 if (!IsPadOn(detector, col-1, row)) {
1006 fGoodTracklet = kFALSE;
1011 if (!IsPadOn(detector, col+1, row)) {
1012 fGoodTracklet = kFALSE;
1017 //_____________________________________________________________________________
1018 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
1021 // Look in the choosen database if the pad is On.
1022 // If no the track will be "not good"
1025 // Get the parameter object
1026 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1028 AliInfo("Could not get calibDB");
1032 if (!cal->IsChamberInstalled(detector) ||
1033 cal->IsChamberMasked(detector) ||
1034 cal->IsPadMasked(detector,col,row)) {
1042 //_____________________________________________________________________________
1043 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
1046 // Calculate the calibration group number for i
1049 // Row of the cluster and position in the pad groups
1051 if (fCalibraMode->GetNnZ(i) != 0) {
1052 posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
1056 // Col of the cluster and position in the pad groups
1058 if (fCalibraMode->GetNnRphi(i) != 0) {
1059 posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
1062 return posc*fCalibraMode->GetNfragZ(i)+posr;
1065 //____________________________________________________________________________________
1066 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1069 // Calculate the total number of calibration groups
1073 fCalibraMode->ModePadCalibration(2,i);
1074 fCalibraMode->ModePadFragmentation(0,2,0,i);
1075 fCalibraMode->SetDetChamb2(i);
1076 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1077 fCalibraMode->ModePadCalibration(0,i);
1078 fCalibraMode->ModePadFragmentation(0,0,0,i);
1079 fCalibraMode->SetDetChamb0(i);
1080 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1081 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1085 //____________Set the pad calibration variables for the detector_______________
1086 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1089 // For the detector calcul the first Xbins and set the number of row
1090 // and col pads per calibration groups, the number of calibration
1091 // groups in the detector.
1094 // first Xbins of the detector
1096 fCalibraMode->CalculXBins(detector,0);
1099 fCalibraMode->CalculXBins(detector,1);
1102 fCalibraMode->CalculXBins(detector,2);
1105 // fragmentation of idect
1106 for (Int_t i = 0; i < 3; i++) {
1107 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1108 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1109 , (Int_t) GetChamber(detector)
1110 , (Int_t) GetSector(detector),i);
1116 //_____________________________________________________________________________
1117 void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group, Int_t *rowcol)
1120 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1123 // Charge in the cluster
1124 Float_t q = TMath::Abs(cl->GetQ());
1125 Int_t time = cl->GetLocalTimeBin();
1127 //Correct for the gain coefficient used in the database for reconstruction
1128 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]);
1129 Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]);
1131 // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458)
1132 Int_t timec = Arrondi((Double_t)(time+correcttheT0));
1133 if((correcttheT0+0.5)==(int(correcttheT0+0.5))) {
1136 if( timec < 0 ) return;
1139 // Correction due to the track angle
1140 Float_t correction = 1.0;
1141 Float_t normalisation = 6.67;
1142 // we divide with gain in AliTRDclusterizerV1::Transform...
1143 if( correctthegain > 0 ) normalisation /= correctthegain;
1144 if ((q >0) && (t->GetNdedx() > 0)) {
1145 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
1148 // Fill the fAmpTotal with the charge
1150 fAmpTotal[(Int_t) group[0]] += correction;
1153 // Fill the fPHPlace and value
1155 fPHPlace[timec] = group[1];
1156 fPHValue[timec] = correction;
1160 //_____________________________________________________________________________
1161 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, AliTRDtrack *t, Int_t index, Int_t *group, Int_t *rowcol)
1164 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1167 // Charge in the cluster
1168 Float_t q = TMath::Abs(cl->GetQ());
1169 Int_t time = cl->GetLocalTimeBin();
1171 //Correct for the gain coefficient used in the database for reconstruction
1172 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]);
1173 Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]);
1175 // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458)
1176 Int_t timec = Arrondi((Double_t)(time+correcttheT0));
1177 if((correcttheT0+0.5)==(int(correcttheT0+0.5))) {
1180 if( timec < 0 ) return;
1182 // Correction due to the track angle
1183 Float_t correction = 1.0;
1184 Float_t normalisation = 6.67;
1185 // we divide with gain in AliTRDclusterizerV1::Transform...
1186 if( correctthegain > 0 ) normalisation /= correctthegain;
1188 correction = t->GetClusterdQdl(index) / (normalisation);
1191 // Fill the fAmpTotal with the charge
1193 fAmpTotal[(Int_t) group[0]] += correction;
1196 // Fill the fPHPlace and value
1198 fPHPlace[timec] = group[1];
1199 fPHValue[timec] = correction;
1203 //_____________________________________________________________________
1204 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamTB *rawStream, Bool_t nocheck)
1207 // Event Processing loop - AliTRDrawStreamTB
1208 // 0 timebin problem
1213 Int_t withInput = 1;
1215 Int_t phvalue[21][36];
1216 for(Int_t k = 0; k < 36; k++){
1217 for(Int_t j = 0; j < 21; j++){
1222 fDetectorPreviousTrack = -1;
1225 Int_t nbtimebin = 0;
1226 Int_t baseline = 10;
1228 // For selecting the signal
1229 Double_t mean[21] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1230 Int_t first[21] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1231 Int_t select[21] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1237 while (rawStream->Next()) {
1239 Int_t idetector = rawStream->GetDet(); // current detector
1240 Int_t imcm = rawStream->GetMCM(); // current MCM
1241 Int_t irob = rawStream->GetROB(); // current ROB
1243 if(((fMCMPrevious != imcm) || (fDetectorPreviousTrack != idetector) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
1245 // take the mean values and check the first time bin
1246 for(Int_t j = 0; j < 21; j++){
1247 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1249 if(phvalue[j][0] > 200.0) first[j] = 1;
1254 for(Int_t j = 1; j < 20; j++){
1255 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1262 for(Int_t j = 1; j < 20; j++){
1265 for(Int_t k = 0; k < fTimeMax; k++){
1266 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1267 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1270 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1271 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1273 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,(3*baseline),fTimeMax);
1280 for(Int_t j = 0; j < 21; j++){
1284 for(Int_t k = 0; k < 36; k++){
1285 phvalue[j][k] = baseline;
1290 fDetectorPreviousTrack = idetector;
1291 fMCMPrevious = imcm;
1292 fROBPrevious = irob;
1294 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1295 if(nbtimebin == 0) return 0;
1296 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
1297 fTimeMax = nbtimebin;
1299 baseline = rawStream->GetCommonAdditive(); // common additive baseline
1301 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
1302 Int_t *signal = rawStream->GetSignals(); // current ADC signal
1303 Int_t col = (rawStream->GetCol())%18; // current COL MCM
1305 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
1307 if((col < 0) || (col >= 21)) return 0;
1308 if((imcm>=16) || (imcm < 0)) return 0;
1310 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
1312 for(Int_t itime = iTimeBin; itime < fin; itime++){
1313 phvalue[col][itime] = signal[n];
1318 // fill the last one
1319 if(fDetectorPreviousTrack != -1){
1321 // take the mean values and check the first time bin
1322 for(Int_t j = 0; j < 21; j++){
1323 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1325 if(phvalue[j][0] > 200.0) first[j] = 1;
1330 for(Int_t j = 1; j < 20; j++){
1331 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1338 for(Int_t j = 1; j < 20; j++){
1341 for(Int_t k = 0; k < fTimeMax; k++){
1342 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1343 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1346 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1347 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1349 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,(3*baseline),fTimeMax);
1356 for(Int_t j = 0; j < 21; j++){
1360 for(Int_t k = 0; k < 36; k++){
1361 phvalue[j][k] = baseline;
1369 while (rawStream->Next()) {
1371 Int_t idetector = rawStream->GetDet(); // current detector
1372 Int_t imcm = rawStream->GetMCM(); // current MCM
1373 Int_t irob = rawStream->GetROB(); // current ROB
1375 if(((fMCMPrevious != imcm) || (fDetectorPreviousTrack != idetector) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
1377 // take the mean values and check the first time bin
1378 for(Int_t j = 0; j < 21; j++){
1379 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1381 if(phvalue[j][0] > 200.0) first[j] = 1;
1383 //printf("meanvalue %f, first %d\n",mean[j],first[j]);
1387 for(Int_t j = 1; j < 20; j++){
1388 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1395 for(Int_t j = 1; j < 20; j++){
1396 //printf("select %d\n",select[j]);
1399 for(Int_t k = 0; k < fTimeMax; k++){
1400 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1402 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1405 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1407 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1409 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1416 for(Int_t j = 0; j < 21; j++){
1420 for(Int_t k = 0; k < 36; k++){
1421 phvalue[j][k] = baseline;
1426 fDetectorPreviousTrack = idetector;
1427 fMCMPrevious = imcm;
1428 fROBPrevious = irob;
1432 baseline = rawStream->GetCommonAdditive(); // common baseline
1434 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1435 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
1436 Int_t *signal = rawStream->GetSignals(); // current ADC signal
1437 Int_t col = (rawStream->GetCol())%18; // current COL MCM
1439 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
1442 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
1445 if((col < 0) || (col >= 21)) return 0;
1446 if((imcm>=16) || (imcm < 0)) return 0;
1448 for(Int_t itime = iTimeBin; itime < fin; itime++){
1449 phvalue[col][itime] = signal[n];
1454 // fill the last one
1455 if(fDetectorPreviousTrack != -1){
1457 // take the mean values and check the first time bin
1458 for(Int_t j = 0; j < 21; j++){
1459 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1461 if(phvalue[j][0] > 200.0) first[j] = 1;
1466 for(Int_t j = 1; j < 20; j++){
1467 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1474 for(Int_t j = 1; j < 20; j++){
1477 for(Int_t k = 0; k < fTimeMax; k++){
1478 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1479 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1482 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1483 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1485 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1492 for(Int_t j = 0; j < 21; j++){
1496 for(Int_t k = 0; k < 36; k++){
1497 phvalue[j][k] = baseline;
1506 //_____________________________________________________________________
1507 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
1510 // Event processing loop - AliRawReader
1514 AliTRDrawStreamTB rawStream(rawReader);
1516 rawReader->Select("TRD");
1518 return ProcessEventDAQ(&rawStream, nocheck);
1520 //_________________________________________________________________________
1521 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
1523 eventHeaderStruct *event,
1526 eventHeaderStruct* /*event*/,
1533 // process date event
1536 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1537 Int_t result=ProcessEventDAQ(rawReader, nocheck);
1541 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1546 //____________Online trackling in AliTRDtrigger________________________________
1547 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1551 // Fill a simple average pulse height
1554 // Localisation of the Xbins involved
1555 //LocalisationDetectorXbins(det);
1557 // Row and position in the pad groups
1559 //if (fCalibraMode->GetNnZ(1) != 0) {
1560 // posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1563 // Col of cluster and position in the pad groups
1565 //if (fCalibraMode->GetNnRphi(1) != 0) {
1566 // posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1569 //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1571 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1576 //____________Write_____________________________________________________
1577 //_____________________________________________________________________
1578 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1581 // Write infos to file
1585 if ( fDebugStreamer ) {
1586 delete fDebugStreamer;
1587 fDebugStreamer = 0x0;
1590 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1595 ,fNumberUsedPh[1]));
1597 TDirectory *backup = gDirectory;
1603 option = "recreate";
1605 TFile f(filename,option.Data());
1607 TStopwatch stopwatch;
1610 f.WriteTObject(fCalibraVector);
1615 f.WriteTObject(fCH2d);
1620 f.WriteTObject(fPH2d);
1625 f.WriteTObject(fPRF2d);
1628 if(fLinearFitterOn){
1629 AnalyseLinearFitter();
1630 f.WriteTObject(fLinearVdriftFit);
1635 if ( backup ) backup->cd();
1637 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1638 ,stopwatch.RealTime(),stopwatch.CpuTime()));
1640 //___________________________________________probe the histos__________________________________________________
1641 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1644 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1645 // debug mode with 2 for TH2I and 3 for TProfile2D
1646 // It gives a pointer to a Double_t[7] with the info following...
1647 // [0] : number of calibration groups with entries
1648 // [1] : minimal number of entries found
1649 // [2] : calibration group number of the min
1650 // [3] : maximal number of entries found
1651 // [4] : calibration group number of the max
1652 // [5] : mean number of entries found
1653 // [6] : mean relative error
1656 Double_t *info = new Double_t[7];
1658 // Number of Xbins (detectors or groups of pads)
1659 Int_t nbins = h->GetNbinsY(); //number of calibration groups
1660 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
1663 Double_t nbwe = 0; //number of calibration groups with entries
1664 Double_t minentries = 0; //minimal number of entries found
1665 Double_t maxentries = 0; //maximal number of entries found
1666 Double_t placemin = 0; //calibration group number of the min
1667 Double_t placemax = -1; //calibration group number of the max
1668 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1669 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1671 Double_t counter = 0;
1674 TH1F *nbEntries = 0x0;//distribution of the number of entries
1675 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
1676 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
1678 // Beginning of the loop over the calibration groups
1679 for (Int_t idect = 0; idect < nbins; idect++) {
1681 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1682 projch->SetDirectory(0);
1684 // Number of entries for this calibration group
1685 Double_t nentries = 0.0;
1687 for (Int_t k = 0; k < nxbins; k++) {
1688 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1692 for (Int_t k = 0; k < nxbins; k++) {
1693 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1694 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1695 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1703 if((!((Bool_t)nbEntries)) && (nentries > 0)){
1704 nbEntries = new TH1F("Number of entries","Number of entries"
1705 ,100,(Int_t)nentries/2,nentries*2);
1706 nbEntries->SetDirectory(0);
1707 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1709 nbEntriesPerGroup->SetDirectory(0);
1710 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1711 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1712 nbEntriesPerSp->SetDirectory(0);
1715 if(nentries > 0) nbEntries->Fill(nentries);
1716 nbEntriesPerGroup->Fill(idect+0.5,nentries);
1717 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1722 if(nentries > maxentries){
1723 maxentries = nentries;
1727 minentries = nentries;
1729 if(nentries < minentries){
1730 minentries = nentries;
1736 meanstats += nentries;
1738 }//calibration groups loop
1740 if(nbwe > 0) meanstats /= nbwe;
1741 if(counter > 0) meanrelativerror /= counter;
1743 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1744 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1745 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1746 AliInfo(Form("The mean number of entries is %f",meanstats));
1747 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1750 info[1] = minentries;
1752 info[3] = maxentries;
1754 info[5] = meanstats;
1755 info[6] = meanrelativerror;
1758 gStyle->SetPalette(1);
1759 gStyle->SetOptStat(1111);
1760 gStyle->SetPadBorderMode(0);
1761 gStyle->SetCanvasColor(10);
1762 gStyle->SetPadLeftMargin(0.13);
1763 gStyle->SetPadRightMargin(0.01);
1764 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1767 nbEntries->Draw("");
1769 nbEntriesPerSp->SetStats(0);
1770 nbEntriesPerSp->Draw("");
1771 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1773 nbEntriesPerGroup->SetStats(0);
1774 nbEntriesPerGroup->Draw("");
1780 //____________________________________________________________________________
1781 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1784 // Return a Int_t[4] with:
1785 // 0 Mean number of entries
1786 // 1 median of number of entries
1787 // 2 rms of number of entries
1788 // 3 number of group with entries
1791 Double_t *stat = new Double_t[4];
1794 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
1795 Double_t *weight = new Double_t[nbofgroups];
1796 Int_t *nonul = new Int_t[nbofgroups];
1798 for(Int_t k = 0; k < nbofgroups; k++){
1799 if(fEntriesCH[k] > 0) {
1801 nonul[(Int_t)stat[3]] = fEntriesCH[k];
1804 else weight[k] = 0.0;
1806 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
1807 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
1808 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1813 //____________________________________________________________________________
1814 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1817 // Return a Int_t[4] with:
1818 // 0 Mean number of entries
1819 // 1 median of number of entries
1820 // 2 rms of number of entries
1821 // 3 number of group with entries
1824 Double_t *stat = new Double_t[4];
1827 Int_t nbofgroups = 540;
1828 Double_t *weight = new Double_t[nbofgroups];
1829 Int_t *nonul = new Int_t[nbofgroups];
1831 for(Int_t k = 0; k < nbofgroups; k++){
1832 if(fEntriesLinearFitter[k] > 0) {
1834 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1837 else weight[k] = 0.0;
1839 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1840 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1841 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1846 //_____________________________________________________________________________
1847 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1850 // Should be between 0 and 6
1853 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1854 AliInfo("The number of groups must be between 0 and 6!");
1857 fNgroupprf = numberGroupsPRF;
1861 //_____________________________________________________________________________
1862 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1865 // Set the factor that will divide the deposited charge
1866 // to fit in the histo range [0,300]
1869 if (RelativeScale > 0.0) {
1870 fRelativeScale = RelativeScale;
1873 AliInfo("RelativeScale must be strict positif!");
1878 //_____________________________________________________________________________
1879 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1882 // Set the mode of calibration group in the z direction for the parameter i
1887 fCalibraMode->SetNz(i, Nz);
1890 AliInfo("You have to choose between 0 and 4");
1895 //_____________________________________________________________________________
1896 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1899 // Set the mode of calibration group in the rphi direction for the parameter i
1904 fCalibraMode->SetNrphi(i ,Nrphi);
1907 AliInfo("You have to choose between 0 and 6");
1911 //____________Protected Functions______________________________________________
1912 //____________Create the 2D histo to be filled online__________________________
1914 //_____________________________________________________________________________
1915 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1918 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1919 // If fNgroupprf is zero then no binning in tnp
1923 name += fCalibraMode->GetNz(2);
1925 name += fCalibraMode->GetNrphi(2);
1929 if(fNgroupprf != 0){
1931 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1932 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1933 fPRF2d->SetYTitle("Det/pad groups");
1934 fPRF2d->SetXTitle("Position x/W [pad width units]");
1935 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1936 fPRF2d->SetStats(0);
1939 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1940 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1941 fPRF2d->SetYTitle("Det/pad groups");
1942 fPRF2d->SetXTitle("Position x/W [pad width units]");
1943 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1944 fPRF2d->SetStats(0);
1949 //_____________________________________________________________________________
1950 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1953 // Create the 2D histos
1957 name += fCalibraMode->GetNz(1);
1959 name += fCalibraMode->GetNrphi(1);
1961 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1962 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1964 fPH2d->SetYTitle("Det/pad groups");
1965 fPH2d->SetXTitle("time [#mus]");
1966 fPH2d->SetZTitle("<PH> [a.u.]");
1970 //_____________________________________________________________________________
1971 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1974 // Create the 2D histos
1978 name += fCalibraMode->GetNz(0);
1980 name += fCalibraMode->GetNrphi(0);
1982 fCH2d = new TH2I("CH2d",(const Char_t *) name
1983 ,fNumberBinCharge,0,300,nn,0,nn);
1984 fCH2d->SetYTitle("Det/pad groups");
1985 fCH2d->SetXTitle("charge deposit [a.u]");
1986 fCH2d->SetZTitle("counts");
1992 //____________Offine tracking in the AliTRDtracker_____________________________
1993 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1996 // For the offline tracking or mcm tracklets
1997 // This function will be called in the functions UpdateHistogram...
1998 // to fill the info of a track for the relativ gain calibration
2001 Int_t nb = 0; // Nombre de zones traversees
2002 Int_t fd = -1; // Premiere zone non nulle
2003 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
2008 // See if the track goes through different zones
2009 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2010 if (fAmpTotal[k] > 0.0) {
2011 totalcharge += fAmpTotal[k];
2024 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2026 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2027 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
2030 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
2034 if ((fAmpTotal[fd] > 0.0) &&
2035 (fAmpTotal[fd+1] > 0.0)) {
2036 // One of the two very big
2037 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2039 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2040 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
2043 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
2046 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2048 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2050 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
2051 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
2054 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
2057 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2060 if (fCalibraMode->GetNfragZ(0) > 1) {
2061 if (fAmpTotal[fd] > 0.0) {
2062 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2063 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2064 // One of the two very big
2065 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2067 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2068 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
2071 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
2074 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2076 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2078 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
2079 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2082 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2084 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
2095 //____________Offine tracking in the AliTRDtracker_____________________________
2096 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2099 // For the offline tracking or mcm tracklets
2100 // This function will be called in the functions UpdateHistogram...
2101 // to fill the info of a track for the drift velocity calibration
2104 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2105 Int_t fd1 = -1; // Premiere zone non nulle
2106 Int_t fd2 = -1; // Deuxieme zone non nulle
2107 Int_t k1 = -1; // Debut de la premiere zone
2108 Int_t k2 = -1; // Debut de la seconde zone
2110 // See if the track goes through different zones
2111 for (Int_t k = 0; k < fTimeMax; k++) {
2112 if (fPHValue[k] > 0.0) {
2117 if (fPHPlace[k] != fd1) {
2123 if (fPHPlace[k] != fd2) {
2135 for (Int_t i = 0; i < fTimeMax; i++) {
2137 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2140 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2145 if ((fd1 == fd2+1) ||
2147 // One of the two fast all the think
2148 if (k2 > (k1+fDifference)) {
2149 //we choose to fill the fd1 with all the values
2151 for (Int_t i = 0; i < fTimeMax; i++) {
2153 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2156 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2160 if ((k2+fDifference) < fTimeMax) {
2161 //we choose to fill the fd2 with all the values
2163 for (Int_t i = 0; i < fTimeMax; i++) {
2165 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2168 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2173 // Two zones voisines sinon rien!
2174 if (fCalibraMode->GetNfragZ(1) > 1) {
2176 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2177 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2178 // One of the two fast all the think
2179 if (k2 > (k1+fDifference)) {
2180 //we choose to fill the fd1 with all the values
2182 for (Int_t i = 0; i < fTimeMax; i++) {
2184 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2187 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2191 if ((k2+fDifference) < fTimeMax) {
2192 //we choose to fill the fd2 with all the values
2194 for (Int_t i = 0; i < fTimeMax; i++) {
2196 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2199 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2205 // Two zones voisines sinon rien!
2207 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2208 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2209 // One of the two fast all the think
2210 if (k2 > (k1 + fDifference)) {
2211 //we choose to fill the fd1 with all the values
2213 for (Int_t i = 0; i < fTimeMax; i++) {
2215 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2218 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2222 if ((k2+fDifference) < fTimeMax) {
2223 //we choose to fill the fd2 with all the values
2225 for (Int_t i = 0; i < fTimeMax; i++) {
2227 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2230 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2241 //____________Offine tracking in the AliTRDtracker_____________________________
2242 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
2245 // For the offline tracking
2246 // This function will be called in the functions UpdateHistogram...
2247 // to fill the find the parameter P1 of a track for the drift velocity calibration
2251 //Number of points: if less than 3 return kFALSE
2252 Int_t npoints = fListClusters->GetEntriesFast();
2253 if(npoints <= 2) return kFALSE;
2256 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2257 Double_t snp = 0.0; // sin angle in the plan yx track
2258 Double_t y = 0.0; // y clusters in the middle of the chamber
2259 Double_t z = 0.0; // z cluster in the middle of the chamber
2260 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2261 Double_t tnp = 0.0; // tan angle in the plan xy track
2262 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2263 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2264 Double_t pointError = 0.0; // error after straight line fit
2265 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
2266 Int_t snpright = 1; // if we took in the middle snp
2267 Int_t crossrow = 0; // if it crosses a pad row
2268 Double_t tiltingangle = 0; // tiltingangle of the pad
2269 Float_t dzdx = 0; // dz/dx now from dz/dl
2270 Int_t nbli = 0; // number linear fitter points
2271 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2273 linearFitterTracklet.StoreData(kFALSE);
2274 linearFitterTracklet.ClearPoints();
2276 //if more than one row
2277 Int_t rowp = -1; // if it crosses a pad row
2280 tiltingangle = padplane->GetTiltingAngle();
2281 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2284 for(Int_t k = 0; k < npoints; k++){
2286 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2287 Double_t ycluster = cl->GetY();
2288 Int_t time = cl->GetLocalTimeBin();
2289 Double_t timeis = time/fSf;
2290 //See if cross two pad rows
2291 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2292 if(k==0) rowp = row;
2293 if(row != rowp) crossrow = 1;
2294 //Take in the middle of the chamber
2296 if(time > (Int_t) 10) {
2298 //if(time < (Int_t) 11) {
2304 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2308 if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
2310 //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
2311 if(nbli <= 2) return kFALSE;
2313 // Do the straight line fit now
2315 linearFitterTracklet.Eval();
2316 linearFitterTracklet.GetParameters(pars);
2317 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2318 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2321 if( TMath::Abs(snp) < 1.){
2322 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2324 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2326 if(fDebugLevel > 0){
2327 if ( !fDebugStreamer ) {
2329 TDirectory *backup = gDirectory;
2330 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2331 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2334 (* fDebugStreamer) << "VDRIFT0"<<
2335 "npoints="<<npoints<<
2339 (* fDebugStreamer) << "VDRIFT"<<
2340 "snpright="<<snpright<<
2341 "npoints="<<npoints<<
2343 "detector="<<detector<<
2352 "crossrow="<<crossrow<<
2353 "errorpar="<<errorpar<<
2354 "pointError="<<pointError<<
2359 if(npoints < fNumberClusters) return kFALSE;
2360 if(snpright == 0) return kFALSE;
2361 if(pointError >= 0.1) return kFALSE;
2362 if(crossrow == 1) return kFALSE;
2364 if(fLinearFitterOn){
2365 //Add to the linear fitter of the detector
2366 if( TMath::Abs(snp) < 1.){
2367 Double_t x = tnp-dzdx*tnt;
2368 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2369 if(fLinearFitterDebugOn) {
2370 fLinearVdriftFit->Update(detector,x,pars[1]);
2372 fEntriesLinearFitter[detector]++;
2375 //AliInfo("End of FindP1TrackPH with success!")
2379 //____________Offine tracking in the AliTRDtracker_____________________________
2380 Bool_t AliTRDCalibraFillHisto::HandlePRF()
2383 // For the offline tracking
2384 // Fit the tracklet with a line and take the position as reference for the PRF
2388 Int_t npoints = fListClusters->GetEntriesFast(); // number of total points
2389 Int_t nb3pc = 0; // number of three pads clusters used for fit
2390 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
2393 // To see the difference due to the fit
2394 Double_t *padPositions;
2395 padPositions = new Double_t[npoints];
2396 for(Int_t k = 0; k < npoints; k++){
2397 padPositions[k] = 0.0;
2401 //Find the position by a fit
2402 TLinearFitter fitter(2,"pol1");
2403 fitter.StoreData(kFALSE);
2404 fitter.ClearPoints();
2405 for(Int_t k = 0; k < npoints; k++){
2407 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2408 Short_t *signals = cl->GetSignals();
2409 Double_t time = cl->GetLocalTimeBin();
2410 //Calculate x if possible
2411 Float_t xcenter = 0.0;
2412 Bool_t echec = kTRUE;
2413 if((time<=7) || (time>=21)) continue;
2414 // Center 3 balanced: position with the center of the pad
2415 if ((((Float_t) signals[3]) > 0.0) &&
2416 (((Float_t) signals[2]) > 0.0) &&
2417 (((Float_t) signals[4]) > 0.0)) {
2419 // Security if the denomiateur is 0
2420 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2421 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2422 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2423 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2424 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2430 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2432 //if no echec: calculate with the position of the pad
2433 // Position of the cluster
2434 Double_t padPosition = xcenter + cl->GetPadCol();
2435 padPositions[k] = padPosition;
2437 fitter.AddPoint(&time, padPosition,1);
2440 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2441 if(nb3pc < 3) return kFALSE;
2444 fitter.GetParameters(line);
2445 Float_t pointError = -1.0;
2446 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
2450 for(Int_t k = 0; k < npoints; k++){
2452 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2453 Short_t *signals = cl->GetSignals(); // signal
2454 Double_t time = cl->GetLocalTimeBin(); // time bin
2455 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
2456 Float_t padPos = cl->GetPadCol(); // middle pad
2457 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2458 Float_t ycenter = 0.0; // relative center charge
2459 Float_t ymin = 0.0; // relative left charge
2460 Float_t ymax = 0.0; // relative right charge
2461 Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx
2462 Double_t pt = fPar4[(Int_t)time]; // pt
2463 Float_t dzdx = 0.0; // dzdx
2466 //Requiere simply two pads clusters at least
2467 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2468 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2469 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2470 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2471 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2472 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2476 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2477 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2478 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
2479 Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan
2480 Float_t xcl = cl->GetY(); // y cluster
2481 Float_t qcl = cl->GetQ(); // charge cluster
2482 Int_t plane = GetPlane(detector); // plane
2483 Int_t chamber = GetChamber(detector); // chamber
2484 Double_t xdiff = dpad; // reconstructed position constant
2485 Double_t x = dpad; // reconstructed position moved
2486 Float_t ep = pointError; // error of fit
2487 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2488 Float_t signal3 = (Float_t)signals[3]; // signal
2489 Float_t signal2 = (Float_t)signals[2]; // signal
2490 Float_t signal4 = (Float_t)signals[4]; // signal
2491 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2493 if(TMath::Abs(snp) < 1.0){
2494 tnp = snp / (TMath::Sqrt(1-snp*snp));
2495 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2499 if(fDebugLevel > 0){
2500 if ( !fDebugStreamer ) {
2502 TDirectory *backup = gDirectory;
2503 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2504 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2507 (* fDebugStreamer) << "PRF0"<<
2508 "caligroup="<<caligroup<<
2509 "detector="<<detector<<
2511 "chamber="<<chamber<<
2512 "npoints="<<npoints<<
2521 "padPosition="<<padPositions[k]<<
2522 "padPosTracklet="<<padPosTracklet<<
2524 "ycenter="<<ycenter<<
2529 "signal1="<<signal1<<
2530 "signal2="<<signal2<<
2531 "signal3="<<signal3<<
2532 "signal4="<<signal4<<
2533 "signal5="<<signal5<<
2538 Float_t y = ycenter;
2539 (* fDebugStreamer) << "PRFALL"<<
2540 "caligroup="<<caligroup<<
2541 "detector="<<detector<<
2543 "chamber="<<chamber<<
2544 "npoints="<<npoints<<
2554 "padPosition="<<padPositions[k]<<
2555 "padPosTracklet="<<padPosTracklet<<
2560 "signal1="<<signal1<<
2561 "signal2="<<signal2<<
2562 "signal3="<<signal3<<
2563 "signal4="<<signal4<<
2564 "signal5="<<signal5<<
2570 (* fDebugStreamer) << "PRFALL"<<
2571 "caligroup="<<caligroup<<
2572 "detector="<<detector<<
2574 "chamber="<<chamber<<
2575 "npoints="<<npoints<<
2585 "padPosition="<<padPositions[k]<<
2586 "padPosTracklet="<<padPosTracklet<<
2591 "signal1="<<signal1<<
2592 "signal2="<<signal2<<
2593 "signal3="<<signal3<<
2594 "signal4="<<signal4<<
2595 "signal5="<<signal5<<
2602 (* fDebugStreamer) << "PRFALL"<<
2603 "caligroup="<<caligroup<<
2604 "detector="<<detector<<
2606 "chamber="<<chamber<<
2607 "npoints="<<npoints<<
2617 "padPosition="<<padPositions[k]<<
2618 "padPosTracklet="<<padPosTracklet<<
2623 "signal1="<<signal1<<
2624 "signal2="<<signal2<<
2625 "signal3="<<signal3<<
2626 "signal4="<<signal4<<
2627 "signal5="<<signal5<<
2634 if(npoints < fNumberClusters) continue;
2635 if(nb3pc <= 5) continue;
2636 if((time >= 21) || (time < 7)) continue;
2637 if(TMath::Abs(snp) >= 1.0) continue;
2638 if(qcl < 80) continue;
2640 Bool_t echec = kFALSE;
2641 Double_t shift = 0.0;
2642 //Calculate the shift in x coresponding to this tnp
2643 if(fNgroupprf != 0.0){
2644 shift = -3.0*(fNgroupprf-1)-1.5;
2645 Double_t limithigh = -0.2*(fNgroupprf-1);
2646 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2648 while(tnp > limithigh){
2654 if (fHisto2d && !echec) {
2655 if(TMath::Abs(dpad) < 1.5) {
2656 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2657 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2659 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2660 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2661 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2663 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2664 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2665 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2668 //Not equivalent anymore here!
2669 if (fVector2d && !echec) {
2670 if(TMath::Abs(dpad) < 1.5) {
2671 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2672 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2674 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2675 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2676 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2678 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2679 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2680 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2687 //____________Offine tracking in the AliTRDtracker_____________________________
2688 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2691 // For the offline tracking
2692 // This function will be called in the functions UpdateHistogram...
2693 // to fill the find the parameter P1 of a track for the drift velocity calibration
2697 //Number of points: if less than 3 return kFALSE
2698 Int_t npoints = index1-index0;
2699 if(npoints <= 2) return kFALSE;
2702 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2703 Double_t snp = 0.0; // sin angle in the plan yx track
2704 Double_t y = 0.0; // y clusters in the middle of the chamber
2705 Double_t z = 0.0; // z cluster in the middle of the chamber
2706 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2707 Double_t tnp = 0.0; // tan angle in the plan xy track
2708 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2709 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2710 Double_t pointError = 0.0; // error after straight line fit
2711 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
2712 //Int_t snpright = 1; // if we took in the middle snp
2713 Int_t crossrow = 0; // if it crosses a pad row
2714 Double_t tiltingangle = 0; // tiltingangle of the pad
2715 Float_t dzdx = 0; // dz/dx now from dz/dl
2716 Int_t nbli = 0; // number linear fitter points
2717 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2719 linearFitterTracklet.StoreData(kFALSE);
2720 linearFitterTracklet.ClearPoints();
2722 //if more than one row
2723 Int_t rowp = -1; // if it crosses a pad row
2726 tiltingangle = padplane->GetTiltingAngle();
2727 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2730 for(Int_t k = 0; k < npoints; k++){
2732 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2733 Double_t ycluster = cl->GetY();
2734 Int_t time = cl->GetLocalTimeBin();
2735 Double_t timeis = time/fSf;
2736 //See if cross two pad rows
2737 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2738 if(k==0) rowp = row;
2739 if(row != rowp) crossrow = 1;
2740 //Take in the middle of the chamber
2742 //if(time > (Int_t) 10) {
2744 if(time < (Int_t) 11) {
2748 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2752 // take now the snp, tnp and tgl from the track
2753 snp = t->GetSnpPlane(GetPlane(detector));
2754 tgl = t->GetTglPlane(GetPlane(detector));
2757 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0;
2759 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0;
2760 if(nbli <= 2) return kFALSE;
2762 // Do the straight line fit now
2764 linearFitterTracklet.Eval();
2765 linearFitterTracklet.GetParameters(pars);
2766 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2767 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2770 if( TMath::Abs(snp) < 1.){
2771 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2773 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2775 if(fDebugLevel > 0){
2776 if ( !fDebugStreamer ) {
2778 TDirectory *backup = gDirectory;
2779 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2780 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2784 (* fDebugStreamer) << "VDRIFT"<<
2785 //"snpright="<<snpright<<
2786 "npoints="<<npoints<<
2788 "detector="<<detector<<
2797 "crossrow="<<crossrow<<
2798 "errorpar="<<errorpar<<
2799 "pointError="<<pointError<<
2804 if(npoints < fNumberClusters) return kFALSE;
2805 //if(snpright == 0) return kFALSE;
2806 if(pointError >= 0.1) return kFALSE;
2807 if(crossrow == 1) return kFALSE;
2809 if(fLinearFitterOn){
2810 //Add to the linear fitter of the detector
2811 if( TMath::Abs(snp) < 1.){
2812 Double_t x = tnp-dzdx*tnt;
2813 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2814 if(fLinearFitterDebugOn) {
2815 fLinearVdriftFit->Update(detector,x,pars[1]);
2817 fEntriesLinearFitter[detector]++;
2820 //AliInfo("End of FindP1TrackPH with success!")
2824 //____________Offine tracking in the AliTRDtracker_____________________________
2825 Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2828 // For the offline tracking
2829 // Fit the tracklet with a line and take the position as reference for the PRF
2833 Int_t npoints = index1-index0; // number of total points
2834 Int_t nb3pc = 0; // number of three pads clusters used for fit
2835 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
2838 // To see the difference due to the fit
2839 Double_t *padPositions;
2840 padPositions = new Double_t[npoints];
2841 for(Int_t k = 0; k < npoints; k++){
2842 padPositions[k] = 0.0;
2846 //Find the position by a fit
2847 TLinearFitter fitter(2,"pol1");
2848 fitter.StoreData(kFALSE);
2849 fitter.ClearPoints();
2850 for(Int_t k = 0; k < npoints; k++){
2852 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2853 Short_t *signals = cl->GetSignals();
2854 Double_t time = cl->GetLocalTimeBin();
2855 //Calculate x if possible
2856 Float_t xcenter = 0.0;
2857 Bool_t echec = kTRUE;
2858 if((time<=7) || (time>=21)) continue;
2859 // Center 3 balanced: position with the center of the pad
2860 if ((((Float_t) signals[3]) > 0.0) &&
2861 (((Float_t) signals[2]) > 0.0) &&
2862 (((Float_t) signals[4]) > 0.0)) {
2864 // Security if the denomiateur is 0
2865 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2866 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2867 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2868 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2869 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2875 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2877 //if no echec: calculate with the position of the pad
2878 // Position of the cluster
2879 Double_t padPosition = xcenter + cl->GetPadCol();
2880 padPositions[k] = padPosition;
2882 fitter.AddPoint(&time, padPosition,1);
2885 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2886 if(nb3pc < 3) return kFALSE;
2889 fitter.GetParameters(line);
2890 Float_t pointError = -1.0;
2891 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
2893 // Take the tgl and snp with the track t now
2894 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
2895 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
2896 Float_t dzdx = 0.0; // dzdx
2898 if(TMath::Abs(snp) < 1.0){
2899 tnp = snp / (TMath::Sqrt(1-snp*snp));
2900 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2905 for(Int_t k = 0; k < npoints; k++){
2907 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2908 Short_t *signals = cl->GetSignals(); // signal
2909 Double_t time = cl->GetLocalTimeBin(); // time bin
2910 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
2911 Float_t padPos = cl->GetPadCol(); // middle pad
2912 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2913 Float_t ycenter = 0.0; // relative center charge
2914 Float_t ymin = 0.0; // relative left charge
2915 Float_t ymax = 0.0; // relative right charge
2919 //Requiere simply two pads clusters at least
2920 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2921 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2922 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2923 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2924 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2925 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2929 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2930 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2931 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
2932 Float_t xcl = cl->GetY(); // y cluster
2933 Float_t qcl = cl->GetQ(); // charge cluster
2934 Int_t plane = GetPlane(detector); // plane
2935 Int_t chamber = GetChamber(detector); // chamber
2936 Double_t xdiff = dpad; // reconstructed position constant
2937 Double_t x = dpad; // reconstructed position moved
2938 Float_t ep = pointError; // error of fit
2939 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2940 Float_t signal3 = (Float_t)signals[3]; // signal
2941 Float_t signal2 = (Float_t)signals[2]; // signal
2942 Float_t signal4 = (Float_t)signals[4]; // signal
2943 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2947 if(fDebugLevel > 0){
2948 if ( !fDebugStreamer ) {
2950 TDirectory *backup = gDirectory;
2951 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2952 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2955 (* fDebugStreamer) << "PRF0"<<
2956 "caligroup="<<caligroup<<
2957 "detector="<<detector<<
2959 "chamber="<<chamber<<
2960 "npoints="<<npoints<<
2968 "padPosition="<<padPositions[k]<<
2969 "padPosTracklet="<<padPosTracklet<<
2971 "ycenter="<<ycenter<<
2976 "signal1="<<signal1<<
2977 "signal2="<<signal2<<
2978 "signal3="<<signal3<<
2979 "signal4="<<signal4<<
2980 "signal5="<<signal5<<
2985 Float_t y = ycenter;
2986 (* fDebugStreamer) << "PRFALL"<<
2987 "caligroup="<<caligroup<<
2988 "detector="<<detector<<
2990 "chamber="<<chamber<<
2991 "npoints="<<npoints<<
3000 "padPosition="<<padPositions[k]<<
3001 "padPosTracklet="<<padPosTracklet<<
3006 "signal1="<<signal1<<
3007 "signal2="<<signal2<<
3008 "signal3="<<signal3<<
3009 "signal4="<<signal4<<
3010 "signal5="<<signal5<<
3016 (* fDebugStreamer) << "PRFALL"<<
3017 "caligroup="<<caligroup<<
3018 "detector="<<detector<<
3020 "chamber="<<chamber<<
3021 "npoints="<<npoints<<
3030 "padPosition="<<padPositions[k]<<
3031 "padPosTracklet="<<padPosTracklet<<
3036 "signal1="<<signal1<<
3037 "signal2="<<signal2<<
3038 "signal3="<<signal3<<
3039 "signal4="<<signal4<<
3040 "signal5="<<signal5<<
3047 (* fDebugStreamer) << "PRFALL"<<
3048 "caligroup="<<caligroup<<
3049 "detector="<<detector<<
3051 "chamber="<<chamber<<
3052 "npoints="<<npoints<<
3061 "padPosition="<<padPositions[k]<<
3062 "padPosTracklet="<<padPosTracklet<<
3067 "signal1="<<signal1<<
3068 "signal2="<<signal2<<
3069 "signal3="<<signal3<<
3070 "signal4="<<signal4<<
3071 "signal5="<<signal5<<
3078 if(npoints < fNumberClusters) continue;
3079 if(nb3pc <= 5) continue;
3080 if((time >= 21) || (time < 7)) continue;
3081 if(TMath::Abs(snp) >= 1.0) continue;
3082 if(qcl < 80) continue;
3084 Bool_t echec = kFALSE;
3085 Double_t shift = 0.0;
3086 //Calculate the shift in x coresponding to this tnp
3087 if(fNgroupprf != 0.0){
3088 shift = -3.0*(fNgroupprf-1)-1.5;
3089 Double_t limithigh = -0.2*(fNgroupprf-1);
3090 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
3092 while(tnp > limithigh){
3098 if (fHisto2d && !echec) {
3099 if(TMath::Abs(dpad) < 1.5) {
3100 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
3101 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
3103 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3104 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
3105 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
3107 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3108 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
3109 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
3112 //Not equivalent anymore here!
3113 if (fVector2d && !echec) {
3114 if(TMath::Abs(dpad) < 1.5) {
3115 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
3116 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
3118 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3119 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
3120 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
3122 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3123 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
3124 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
3132 //____________Some basic geometry function_____________________________________
3134 //_____________________________________________________________________________
3135 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3138 // Reconstruct the plane number from the detector number
3141 return ((Int_t) (d % 6));
3145 //_____________________________________________________________________________
3146 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3149 // Reconstruct the chamber number from the detector number
3153 return ((Int_t) (d % 30) / fgkNplan);
3157 //_____________________________________________________________________________
3158 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3161 // Reconstruct the sector number from the detector number
3165 return ((Int_t) (d / fg));
3168 //_____________________________________________________________________
3169 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3172 // return pointer to fPH2d TProfile2D
3173 // create a new TProfile2D if it doesn't exist allready
3179 fTimeMax = nbtimebin;
3180 fSf = samplefrequency;
3183 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
3184 ,fCalibraMode->GetNz(1)
3185 ,fCalibraMode->GetNrphi(1)));
3187 // Calcul the number of Xbins
3189 fCalibraMode->ModePadCalibration(2,1);
3190 fCalibraMode->ModePadFragmentation(0,2,0,1);
3191 fCalibraMode->SetDetChamb2(1);
3192 ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
3193 fCalibraMode->ModePadCalibration(0,1);
3194 fCalibraMode->ModePadFragmentation(0,0,0,1);
3195 fCalibraMode->SetDetChamb0(1);
3196 ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
3197 AliInfo(Form("Total number of Xbins: %d",ntotal1));
3199 CreatePH2d(ntotal1);
3206 //_____________________________________________________________________
3207 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3210 // return pointer to TLinearFitter Calibration
3211 // if force is true create a new TLinearFitter if it doesn't exist allready
3214 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3215 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3218 // if we are forced and TLinearFitter doesn't yet exist create it
3220 // new TLinearFitter
3221 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3222 fLinearFitterArray.AddAt(linearfitter,detector);
3223 return linearfitter;
3226 //_____________________________________________________________________
3227 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3230 // FillCH2d: Marian style
3233 //skip simply the value out of range
3234 if((y>=300.0) || (y<0.0)) return;
3236 //Calcul the y place
3237 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3238 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3241 fCH2d->GetArray()[place]++;
3245 //____________________________________________________________________________
3246 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3249 // Analyse array of linear fitter because can not be written
3250 // Store two arrays: one with the param the other one with the error param + number of entries
3253 for(Int_t k = 0; k < 540; k++){
3254 TLinearFitter *linearfitter = GetLinearFitter(k);
3255 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3256 TVectorD *par = new TVectorD(2);
3257 TVectorD pare = TVectorD(2);
3258 TVectorD *parE = new TVectorD(3);
3259 linearfitter->Eval();
3260 linearfitter->GetParameters(*par);
3261 linearfitter->GetErrors(pare);
3262 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3263 (*parE)[0] = pare[0]*ppointError;
3264 (*parE)[1] = pare[1]*ppointError;
3265 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3266 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3267 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3271 //________________________________________________________________________________
3272 Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const
3274 // Partie entiere of the (x+0.5)
3279 //if (x + 0.5 == Float_t(i) && i & 1) i--;
3282 //if (x - 0.5 == Float_t(i) && i & 1) i++;