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 "AliTRDpadPlane.h"
61 #include "AliTRDcluster.h"
62 #include "AliTRDtrack.h"
63 #include "AliTRDtrackV1.h"
64 #include "AliTRDrawStreamBase.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 ,fMcmCorrectAngle(kFALSE)
128 ,fLinearFitterOn(kFALSE)
129 ,fLinearFitterDebugOn(kFALSE)
131 ,fThresholdClusterPRF2(15.0)
132 ,fLimitChargeIntegration(kFALSE)
133 ,fFillWithZero(kFALSE)
134 ,fNormalizeNbOfCluster(kFALSE)
137 ,fCalibraMode(new AliTRDCalibraMode())
140 ,fDetectorPreviousTrack(-1)
144 ,fNumberClustersf(30)
150 ,fNumberBinCharge(100)
156 ,fGoodTracklet(kTRUE)
157 ,fLinearFitterTracklet(0x0)
159 ,fEntriesLinearFitter(0x0)
164 ,fLinearFitterArray(540)
165 ,fLinearVdriftFit(0x0)
170 // Default constructor
174 // Init some default values
177 fNumberUsedCh[0] = 0;
178 fNumberUsedCh[1] = 0;
179 fNumberUsedPh[0] = 0;
180 fNumberUsedPh[1] = 0;
182 fGeo = new AliTRDgeometry();
186 //______________________________________________________________________________________
187 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
191 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
194 ,fPRF2dOn(c.fPRF2dOn)
195 ,fHisto2d(c.fHisto2d)
196 ,fVector2d(c.fVector2d)
197 ,fLinearFitterOn(c.fLinearFitterOn)
198 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
199 ,fRelativeScale(c.fRelativeScale)
200 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
201 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
202 ,fFillWithZero(c.fFillWithZero)
203 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
204 ,fMaxCluster(c.fMaxCluster)
205 ,fNbMaxCluster(c.fNbMaxCluster)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
210 ,fMCMPrevious(c.fMCMPrevious)
211 ,fROBPrevious(c.fROBPrevious)
212 ,fNumberClusters(c.fNumberClusters)
213 ,fNumberClustersf(c.fNumberClustersf)
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)
225 ,fGoodTracklet(c.fGoodTracklet)
226 ,fLinearFitterTracklet(0x0)
228 ,fEntriesLinearFitter(0x0)
233 ,fLinearFitterArray(540)
234 ,fLinearVdriftFit(0x0)
241 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
242 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
244 fPH2d = new TProfile2D(*c.fPH2d);
245 fPH2d->SetDirectory(0);
248 fPRF2d = new TProfile2D(*c.fPRF2d);
249 fPRF2d->SetDirectory(0);
252 fCH2d = new TH2I(*c.fCH2d);
253 fCH2d->SetDirectory(0);
255 if(c.fLinearVdriftFit){
256 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
259 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
260 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
265 fGeo = new AliTRDgeometry();
268 //____________________________________________________________________________________
269 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
272 // AliTRDCalibraFillHisto destructor
276 if ( fDebugStreamer ) delete fDebugStreamer;
278 if ( fCalDetGain ) delete fCalDetGain;
279 if ( fCalROCGain ) delete fCalROCGain;
281 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
285 delete [] fEntriesCH;
286 delete [] fEntriesLinearFitter;
289 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
290 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
298 //_____________________________________________________________________________
299 void AliTRDCalibraFillHisto::Destroy()
310 //_____________________________________________________________________________
311 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
314 // Delete DebugStreamer
317 if ( fDebugStreamer ) delete fDebugStreamer;
318 fDebugStreamer = 0x0;
321 //_____________________________________________________________________________
322 void AliTRDCalibraFillHisto::ClearHistos()
342 //////////////////////////////////////////////////////////////////////////////////
343 // calibration with AliTRDtrackV1: Init, Update
344 //////////////////////////////////////////////////////////////////////////////////
345 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
346 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
349 // Init the histograms and stuff to be filled
354 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
356 AliInfo("Could not get calibDB");
359 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
361 AliInfo("Could not get CommonParam");
366 fTimeMax = cal->GetNumberOfTimeBins();
367 fSf = parCom->GetSamplingFrequency();
368 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
369 else fRelativeScale = 1.18;
370 fNumberClustersf = fTimeMax;
371 fNumberClusters = (Int_t)(0.5*fTimeMax);
373 // Init linear fitter
374 if(!fLinearFitterTracklet) {
375 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
376 fLinearFitterTracklet->StoreData(kTRUE);
379 //calib object from database used for reconstruction
381 fCalDetGain->~AliTRDCalDet();
382 new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
383 }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
385 // Calcul Xbins Chambd0, Chamb2
386 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
387 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
388 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
390 // If vector method On initialised all the stuff
392 fCalibraVector = new AliTRDCalibraVector();
393 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
394 fCalibraVector->SetTimeMax(fTimeMax);
395 if(fNgroupprf != 0) {
396 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
397 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
400 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
401 fCalibraVector->SetPRFRange(1.5);
403 for(Int_t k = 0; k < 3; k++){
404 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
405 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
407 TString namech("Nz");
408 namech += fCalibraMode->GetNz(0);
410 namech += fCalibraMode->GetNrphi(0);
411 fCalibraVector->SetNameCH((const char* ) namech);
412 TString nameph("Nz");
413 nameph += fCalibraMode->GetNz(1);
415 nameph += fCalibraMode->GetNrphi(1);
416 fCalibraVector->SetNamePH((const char* ) nameph);
417 TString nameprf("Nz");
418 nameprf += fCalibraMode->GetNz(2);
420 nameprf += fCalibraMode->GetNrphi(2);
422 nameprf += fNgroupprf;
423 fCalibraVector->SetNamePRF((const char* ) nameprf);
426 // Create the 2D histos corresponding to the pad groupCalibration mode
429 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
430 ,fCalibraMode->GetNz(0)
431 ,fCalibraMode->GetNrphi(0)));
433 // Create the 2D histo
438 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
439 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
443 fEntriesCH = new Int_t[ntotal0];
444 for(Int_t k = 0; k < ntotal0; k++){
451 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
452 ,fCalibraMode->GetNz(1)
453 ,fCalibraMode->GetNrphi(1)));
455 // Create the 2D histo
460 fPHPlace = new Short_t[fTimeMax];
461 for (Int_t k = 0; k < fTimeMax; k++) {
464 fPHValue = new Float_t[fTimeMax];
465 for (Int_t k = 0; k < fTimeMax; k++) {
469 if (fLinearFitterOn) {
470 //fLinearFitterArray.Expand(540);
471 fLinearFitterArray.SetName("ArrayLinearFitters");
472 fEntriesLinearFitter = new Int_t[540];
473 for(Int_t k = 0; k < 540; k++){
474 fEntriesLinearFitter[k] = 0;
476 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
481 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
482 ,fCalibraMode->GetNz(2)
483 ,fCalibraMode->GetNrphi(2)));
484 // Create the 2D histo
486 CreatePRF2d(ntotal2);
493 //____________Offline tracking in the AliTRDtracker____________________________
494 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
497 // Use AliTRDtrack for the calibration
501 AliTRDcluster *cl = 0x0; // pointeur to cluster
502 Int_t index0 = 0; // index of the first cluster in the current chamber
503 Int_t index1 = 0; // index of the current cluster in the current chamber
505 //////////////////////////////
506 // loop over the clusters
507 ///////////////////////////////
508 while((cl = t->GetCluster(index1))){
510 /////////////////////////////////////////////////////////////////////////
511 // Fill the infos for the previous clusters if not the same detector
512 ////////////////////////////////////////////////////////////////////////
513 Int_t detector = cl->GetDetector();
514 if ((detector != fDetectorPreviousTrack) &&
515 (index0 != index1)) {
519 //If the same track, then look if the previous detector is in
520 //the same plane, if yes: not a good track
521 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
525 // Fill only if the track doesn't touch a masked pad or doesn't
528 // drift velocity unables to cut bad tracklets
529 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
533 FillTheInfoOfTheTrackCH(index1-index0);
538 FillTheInfoOfTheTrackPH();
541 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
544 } // if a good tracklet
547 ResetfVariablestracklet();
550 } // Fill at the end the charge
553 /////////////////////////////////////////////////////////////
554 // Localise and take the calib gain object for the detector
555 ////////////////////////////////////////////////////////////
556 if (detector != fDetectorPreviousTrack) {
558 //Localise the detector
559 LocalisationDetectorXbins(detector);
562 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
564 AliInfo("Could not get calibDB");
571 fCalROCGain->~AliTRDCalROC();
572 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
574 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
578 // Reset the detectbjobsor
579 fDetectorPreviousTrack = detector;
581 ///////////////////////////////////////
582 // Store the info of the cluster
583 ///////////////////////////////////////
584 Int_t row = cl->GetPadRow();
585 Int_t col = cl->GetPadCol();
586 CheckGoodTrackletV1(cl);
587 Int_t group[2] = {0,0};
588 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
589 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
590 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
594 } // while on clusters
596 ///////////////////////////
597 // Fill the last plane
598 //////////////////////////
599 if( index0 != index1 ){
605 // drift velocity unables to cut bad tracklets
606 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
610 FillTheInfoOfTheTrackCH(index1-index0);
615 FillTheInfoOfTheTrackPH();
618 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
620 } // if a good tracklet
625 ResetfVariablestracklet();
630 //____________Offline tracking in the AliTRDtracker____________________________
631 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
634 // Use AliTRDtrackV1 for the calibration
638 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
639 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
640 Bool_t newtr = kTRUE; // new track
643 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
645 AliInfo("Could not get calibDB");
649 ///////////////////////////
650 // loop over the tracklet
651 ///////////////////////////
652 for(Int_t itr = 0; itr < 6; itr++){
654 if(!(tracklet = t->GetTracklet(itr))) continue;
655 if(!tracklet->IsOK()) continue;
657 ResetfVariablestracklet();
659 //////////////////////////////////////////
660 // localisation of the tracklet and dqdl
661 //////////////////////////////////////////
662 Int_t layer = tracklet->GetPlane();
664 while(!(cl = tracklet->GetClusters(ic++))) continue;
665 Int_t detector = cl->GetDetector();
666 if (detector != fDetectorPreviousTrack) {
667 // if not a new track
669 // don't use the rest of this track if in the same plane
670 if (layer == GetLayer(fDetectorPreviousTrack)) {
671 //printf("bad tracklet, same layer for detector %d\n",detector);
675 //Localise the detector bin
676 LocalisationDetectorXbins(detector);
680 fCalROCGain->~AliTRDCalROC();
681 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
683 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
686 fDetectorPreviousTrack = detector;
690 ////////////////////////////
691 // loop over the clusters
692 ////////////////////////////
693 Int_t nbclusters = 0;
694 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
695 if(!(cl = tracklet->GetClusters(jc))) continue;
698 // Store the info bis of the tracklet
699 Int_t row = cl->GetPadRow();
700 Int_t col = cl->GetPadCol();
701 CheckGoodTrackletV1(cl);
702 Int_t group[2] = {0,0};
703 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
704 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
705 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col);
708 ////////////////////////////////////////
709 // Fill the stuffs if a good tracklet
710 ////////////////////////////////////////
713 // drift velocity unables to cut bad tracklets
714 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
718 FillTheInfoOfTheTrackCH(nbclusters);
723 FillTheInfoOfTheTrackPH();
726 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
728 } // if a good tracklet
734 ///////////////////////////////////////////////////////////////////////////////////
735 // Routine inside the update with AliTRDtrack
736 ///////////////////////////////////////////////////////////////////////////////////
737 //____________Offine tracking in the AliTRDtracker_____________________________
738 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
741 // Drift velocity calibration:
742 // Fit the clusters with a straight line
743 // From the slope find the drift velocity
746 //Number of points: if less than 3 return kFALSE
747 Int_t npoints = index1-index0;
748 if(npoints <= 2) return kFALSE;
753 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
754 // parameters of the track
755 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
756 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
757 Double_t tnp = 0.0; // tan angle in the plan xy track
758 if( TMath::Abs(snp) < 1.){
759 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
761 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
762 // tilting pad and cross row
763 Int_t crossrow = 0; // if it crosses a pad row
764 Int_t rowp = -1; // if it crosses a pad row
765 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
766 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
767 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
769 fLinearFitterTracklet->ClearPoints();
770 Double_t dydt = 0.0; // dydt tracklet after straight line fit
771 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
772 Double_t pointError = 0.0; // error after straight line fit
773 Int_t nbli = 0; // number linear fitter points
775 //////////////////////////////
776 // loop over clusters
777 ////////////////////////////
778 for(Int_t k = 0; k < npoints; k++){
780 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
781 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
782 Double_t ycluster = cl->GetY();
783 Int_t time = cl->GetPadTime();
784 Double_t timeis = time/fSf;
785 //Double_t q = cl->GetQ();
786 //See if cross two pad rows
787 Int_t row = cl->GetPadRow();
789 if(row != rowp) crossrow = 1;
791 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
796 //////////////////////////////
798 //////////////////////////////
800 fLinearFitterTracklet->ClearPoints();
804 fLinearFitterTracklet->Eval();
805 fLinearFitterTracklet->GetParameters(pars);
806 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
807 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
809 fLinearFitterTracklet->ClearPoints();
811 /////////////////////////////
813 ////////////////////////////
815 if ( !fDebugStreamer ) {
817 TDirectory *backup = gDirectory;
818 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
819 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
823 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
824 //"snpright="<<snpright<<
825 "npoints="<<npoints<<
827 "detector="<<detector<<
834 "crossrow="<<crossrow<<
835 "errorpar="<<errorpar<<
836 "pointError="<<pointError<<
840 Int_t nbclusters = index1-index0;
841 Int_t layer = GetLayer(fDetectorPreviousTrack);
843 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
844 //"snpright="<<snpright<<
845 "nbclusters="<<nbclusters<<
846 "detector="<<fDetectorPreviousTrack<<
852 ///////////////////////////
854 ///////////////////////////
855 if(npoints < fNumberClusters) return kFALSE;
856 if(npoints > fNumberClustersf) return kFALSE;
857 if(pointError >= 0.1) return kFALSE;
858 if(crossrow == 1) return kFALSE;
860 ////////////////////////////
862 ////////////////////////////
864 //Add to the linear fitter of the detector
865 if( TMath::Abs(snp) < 1.){
866 Double_t x = tnp-dzdx*tnt;
867 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
868 if(fLinearFitterDebugOn) {
869 fLinearVdriftFit->Update(detector,x,pars[1]);
871 fEntriesLinearFitter[detector]++;
877 //____________Offine tracking in the AliTRDtracker_____________________________
878 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
881 // Drift velocity calibration:
882 // Fit the clusters with a straight line
883 // From the slope find the drift velocity
886 ////////////////////////////////////////////////
887 //Number of points: if less than 3 return kFALSE
888 /////////////////////////////////////////////////
889 if(nbclusters <= 2) return kFALSE;
894 // results of the linear fit
895 Double_t dydt = 0.0; // dydt tracklet after straight line fit
896 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
897 Double_t pointError = 0.0; // error after straight line fit
898 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
899 Int_t crossrow = 0; // if it crosses a pad row
900 Int_t rowp = -1; // if it crosses a pad row
901 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
902 fLinearFitterTracklet->ClearPoints();
905 ///////////////////////////////////////////
906 // Take the parameters of the track
907 //////////////////////////////////////////
908 // take now the snp, tnp and tgl from the track
909 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
910 Double_t tnp = 0.0; // dy/dx at the end of the chamber
911 if( TMath::Abs(snp) < 1.){
912 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
914 Double_t tgl = tracklet->GetTgl(); // dz/dl
915 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
917 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
918 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
919 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
920 // at the end with correction due to linear fit
921 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
922 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
925 ////////////////////////////
926 // loop over the clusters
927 ////////////////////////////
929 AliTRDcluster *cl = 0x0;
930 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
931 if(!(cl = tracklet->GetClusters(ic))) continue;
932 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
934 Double_t ycluster = cl->GetY();
935 Int_t time = cl->GetPadTime();
936 Double_t timeis = time/fSf;
937 //See if cross two pad rows
938 Int_t row = cl->GetPadRow();
939 if(rowp==-1) rowp = row;
940 if(row != rowp) crossrow = 1;
942 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
945 //////////////////////////////
946 // Check no shared clusters
947 //////////////////////////////
948 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
949 if((cl = tracklet->GetClusters(icc))) crossrow = 1;
953 ////////////////////////////////////
954 // Do the straight line fit now
955 ///////////////////////////////////
957 fLinearFitterTracklet->ClearPoints();
961 fLinearFitterTracklet->Eval();
962 fLinearFitterTracklet->GetParameters(pars);
963 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
964 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
966 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
967 fLinearFitterTracklet->ClearPoints();
969 ////////////////////////////////
971 ///////////////////////////////
975 if ( !fDebugStreamer ) {
977 TDirectory *backup = gDirectory;
978 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
979 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
983 Int_t layer = GetLayer(fDetectorPreviousTrack);
985 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
986 //"snpright="<<snpright<<
988 "nbclusters="<<nbclusters<<
989 "detector="<<fDetectorPreviousTrack<<
997 "crossrow="<<crossrow<<
998 "errorpar="<<errorpar<<
999 "pointError="<<pointError<<
1004 /////////////////////////
1006 ////////////////////////
1008 if(nbclusters < fNumberClusters) return kFALSE;
1009 if(nbclusters > fNumberClustersf) return kFALSE;
1010 if(pointError >= 0.3) return kFALSE;
1011 if(crossrow == 1) return kFALSE;
1013 ///////////////////////
1015 //////////////////////
1017 if(fLinearFitterOn){
1018 //Add to the linear fitter of the detector
1019 if( TMath::Abs(snp) < 1.){
1020 Double_t x = tnp-dzdx*tnt;
1021 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1022 if(fLinearFitterDebugOn) {
1023 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1025 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1031 //____________Offine tracking in the AliTRDtracker_____________________________
1032 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1035 // PRF width calibration
1036 // Assume a Gaussian shape: determinate the position of the three pad clusters
1037 // Fit with a straight line
1038 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1039 // Fill the PRF as function of angle of the track
1044 //////////////////////////
1046 /////////////////////////
1047 Int_t npoints = index1-index0; // number of total points
1048 Int_t nb3pc = 0; // number of three pads clusters used for fit
1049 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1050 // To see the difference due to the fit
1051 Double_t *padPositions;
1052 padPositions = new Double_t[npoints];
1053 for(Int_t k = 0; k < npoints; k++){
1054 padPositions[k] = 0.0;
1056 // Take the tgl and snp with the track t now
1057 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1058 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1059 Float_t dzdx = 0.0; // dzdx
1061 if(TMath::Abs(snp) < 1.0){
1062 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1063 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1066 fLinearFitterTracklet->ClearPoints();
1068 ///////////////////////////
1069 // calcul the tnp group
1070 ///////////////////////////
1071 Bool_t echec = kFALSE;
1072 Double_t shift = 0.0;
1073 //Calculate the shift in x coresponding to this tnp
1074 if(fNgroupprf != 0.0){
1075 shift = -3.0*(fNgroupprf-1)-1.5;
1076 Double_t limithigh = -0.2*(fNgroupprf-1);
1077 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1079 while(tnp > limithigh){
1086 delete [] padPositions;
1090 //////////////////////
1092 /////////////////////
1093 for(Int_t k = 0; k < npoints; k++){
1095 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1096 Short_t *signals = cl->GetSignals();
1097 Double_t time = cl->GetLocalTimeBin();
1098 //Calculate x if possible
1099 Float_t xcenter = 0.0;
1100 Bool_t echec1 = kTRUE;
1101 if((time<=7) || (time>=21)) continue;
1102 // Center 3 balanced: position with the center of the pad
1103 if ((((Float_t) signals[3]) > 0.0) &&
1104 (((Float_t) signals[2]) > 0.0) &&
1105 (((Float_t) signals[4]) > 0.0)) {
1107 // Security if the denomiateur is 0
1108 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1109 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1110 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1111 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1112 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1118 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1120 //if no echec: calculate with the position of the pad
1121 // Position of the cluster
1122 Double_t padPosition = xcenter + cl->GetPadCol();
1123 padPositions[k] = padPosition;
1125 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1129 /////////////////////////////
1131 ////////////////////////////
1133 delete [] padPositions;
1134 fLinearFitterTracklet->ClearPoints();
1137 fLinearFitterTracklet->Eval();
1139 fLinearFitterTracklet->GetParameters(line);
1140 Float_t pointError = -1.0;
1141 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1142 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1144 fLinearFitterTracklet->ClearPoints();
1146 /////////////////////////////////////////////////////
1147 // Now fill the PRF: second loop over clusters
1148 /////////////////////////////////////////////////////
1149 for(Int_t k = 0; k < npoints; k++){
1151 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1152 Short_t *signals = cl->GetSignals(); // signal
1153 Double_t time = cl->GetLocalTimeBin(); // time bin
1154 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1155 Float_t padPos = cl->GetPadCol(); // middle pad
1156 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1157 Float_t ycenter = 0.0; // relative center charge
1158 Float_t ymin = 0.0; // relative left charge
1159 Float_t ymax = 0.0; // relative right charge
1161 //Requiere simply two pads clusters at least
1162 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1163 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1164 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1165 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1166 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1167 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1171 Int_t rowcl = cl->GetPadRow(); // row of cluster
1172 Int_t colcl = cl->GetPadCol(); // col of cluster
1173 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1174 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1175 Float_t xcl = cl->GetY(); // y cluster
1176 Float_t qcl = cl->GetQ(); // charge cluster
1177 Int_t layer = GetLayer(detector); // layer
1178 Int_t stack = GetStack(detector); // stack
1179 Double_t xdiff = dpad; // reconstructed position constant
1180 Double_t x = dpad; // reconstructed position moved
1181 Float_t ep = pointError; // error of fit
1182 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1183 Float_t signal3 = (Float_t)signals[3]; // signal
1184 Float_t signal2 = (Float_t)signals[2]; // signal
1185 Float_t signal4 = (Float_t)signals[4]; // signal
1186 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1188 //////////////////////////////
1190 /////////////////////////////
1191 if(fDebugLevel > 0){
1192 if ( !fDebugStreamer ) {
1194 TDirectory *backup = gDirectory;
1195 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1196 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1202 Float_t y = ycenter;
1203 (* fDebugStreamer) << "HandlePRFtracklet"<<
1204 "caligroup="<<caligroup<<
1205 "detector="<<detector<<
1208 "npoints="<<npoints<<
1217 "padPosition="<<padPositions[k]<<
1218 "padPosTracklet="<<padPosTracklet<<
1223 "signal1="<<signal1<<
1224 "signal2="<<signal2<<
1225 "signal3="<<signal3<<
1226 "signal4="<<signal4<<
1227 "signal5="<<signal5<<
1233 (* fDebugStreamer) << "HandlePRFtracklet"<<
1234 "caligroup="<<caligroup<<
1235 "detector="<<detector<<
1238 "npoints="<<npoints<<
1247 "padPosition="<<padPositions[k]<<
1248 "padPosTracklet="<<padPosTracklet<<
1253 "signal1="<<signal1<<
1254 "signal2="<<signal2<<
1255 "signal3="<<signal3<<
1256 "signal4="<<signal4<<
1257 "signal5="<<signal5<<
1263 (* fDebugStreamer) << "HandlePRFtracklet"<<
1264 "caligroup="<<caligroup<<
1265 "detector="<<detector<<
1268 "npoints="<<npoints<<
1277 "padPosition="<<padPositions[k]<<
1278 "padPosTracklet="<<padPosTracklet<<
1283 "signal1="<<signal1<<
1284 "signal2="<<signal2<<
1285 "signal3="<<signal3<<
1286 "signal4="<<signal4<<
1287 "signal5="<<signal5<<
1293 ////////////////////////////
1295 ///////////////////////////
1296 if(npoints < fNumberClusters) continue;
1297 if(npoints > fNumberClustersf) continue;
1298 if(nb3pc <= 5) continue;
1299 if((time >= 21) || (time < 7)) continue;
1300 if(TMath::Abs(snp) >= 1.0) continue;
1301 if(TMath::Abs(qcl) < 80) continue;
1303 ////////////////////////////
1305 ///////////////////////////
1307 if(TMath::Abs(dpad) < 1.5) {
1308 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1309 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1311 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1312 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1313 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1315 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1316 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1317 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1321 if(TMath::Abs(dpad) < 1.5) {
1322 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1323 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1325 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1326 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1327 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1329 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1330 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1331 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1335 delete [] padPositions;
1339 //____________Offine tracking in the AliTRDtracker_____________________________
1340 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1343 // PRF width calibration
1344 // Assume a Gaussian shape: determinate the position of the three pad clusters
1345 // Fit with a straight line
1346 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1347 // Fill the PRF as function of angle of the track
1351 //printf("begin\n");
1352 ///////////////////////////////////////////
1353 // Take the parameters of the track
1354 //////////////////////////////////////////
1355 // take now the snp, tnp and tgl from the track
1356 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1357 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1358 if( TMath::Abs(snp) < 1.){
1359 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1361 Double_t tgl = tracklet->GetTgl(); // dz/dl
1362 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1364 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1365 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1366 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1367 // at the end with correction due to linear fit
1368 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1369 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1371 ///////////////////////////////
1372 // Calculate tnp group shift
1373 ///////////////////////////////
1374 Bool_t echec = kFALSE;
1375 Double_t shift = 0.0;
1376 //Calculate the shift in x coresponding to this tnp
1377 if(fNgroupprf != 0.0){
1378 shift = -3.0*(fNgroupprf-1)-1.5;
1379 Double_t limithigh = -0.2*(fNgroupprf-1);
1380 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1382 while(tnp > limithigh){
1388 // do nothing if out of tnp range
1389 //printf("echec %d\n",(Int_t)echec);
1390 if(echec) return kFALSE;
1392 ///////////////////////
1394 //////////////////////
1396 Int_t nb3pc = 0; // number of three pads clusters used for fit
1397 // to see the difference between the fit and the 3 pad clusters position
1398 Double_t padPositions[AliTRDseedV1::kNtb];
1399 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1400 fLinearFitterTracklet->ClearPoints();
1402 //printf("loop clusters \n");
1403 ////////////////////////////
1404 // loop over the clusters
1405 ////////////////////////////
1406 AliTRDcluster *cl = 0x0;
1407 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1408 // reject shared clusters on pad row
1409 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1411 Double_t time = cl->GetLocalTimeBin();
1412 if((time<=7) || (time>=21)) continue;
1413 Short_t *signals = cl->GetSignals();
1414 Float_t xcenter = 0.0;
1415 Bool_t echec1 = kTRUE;
1417 /////////////////////////////////////////////////////////////
1418 // Center 3 balanced: position with the center of the pad
1419 /////////////////////////////////////////////////////////////
1420 if ((((Float_t) signals[3]) > 0.0) &&
1421 (((Float_t) signals[2]) > 0.0) &&
1422 (((Float_t) signals[4]) > 0.0)) {
1424 // Security if the denomiateur is 0
1425 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1426 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1427 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1428 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1429 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1435 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1436 if(echec1) continue;
1438 ////////////////////////////////////////////////////////
1439 //if no echec1: calculate with the position of the pad
1440 // Position of the cluster
1441 // fill the linear fitter
1442 ///////////////////////////////////////////////////////
1443 Double_t padPosition = xcenter + cl->GetPadCol();
1444 padPositions[ic] = padPosition;
1446 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1451 //printf("Fin loop clusters \n");
1452 //////////////////////////////
1453 // fit with a straight line
1454 /////////////////////////////
1456 fLinearFitterTracklet->ClearPoints();
1459 fLinearFitterTracklet->Eval();
1461 fLinearFitterTracklet->GetParameters(line);
1462 Float_t pointError = -1.0;
1463 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1464 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1466 fLinearFitterTracklet->ClearPoints();
1468 //printf("PRF second loop \n");
1469 ////////////////////////////////////////////////
1470 // Fill the PRF: Second loop over clusters
1471 //////////////////////////////////////////////
1472 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1473 // reject shared clusters on pad row
1474 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1476 if(!(cl = tracklet->GetClusters(ic))) continue;
1478 Short_t *signals = cl->GetSignals(); // signal
1479 Double_t time = cl->GetLocalTimeBin(); // time bin
1480 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1481 Float_t padPos = cl->GetPadCol(); // middle pad
1482 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1483 Float_t ycenter = 0.0; // relative center charge
1484 Float_t ymin = 0.0; // relative left charge
1485 Float_t ymax = 0.0; // relative right charge
1487 ////////////////////////////////////////////////////////////////
1488 // Calculate ycenter, ymin and ymax even for two pad clusters
1489 ////////////////////////////////////////////////////////////////
1490 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1491 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1492 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1493 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1494 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1495 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1498 /////////////////////////
1499 // Calibration group
1500 ////////////////////////
1501 Int_t rowcl = cl->GetPadRow(); // row of cluster
1502 Int_t colcl = cl->GetPadCol(); // col of cluster
1503 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1504 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1505 Float_t xcl = cl->GetY(); // y cluster
1506 Float_t qcl = cl->GetQ(); // charge cluster
1507 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1508 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1509 Double_t xdiff = dpad; // reconstructed position constant
1510 Double_t x = dpad; // reconstructed position moved
1511 Float_t ep = pointError; // error of fit
1512 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1513 Float_t signal3 = (Float_t)signals[3]; // signal
1514 Float_t signal2 = (Float_t)signals[2]; // signal
1515 Float_t signal4 = (Float_t)signals[4]; // signal
1516 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1520 /////////////////////
1522 ////////////////////
1524 if(fDebugLevel > 0){
1525 if ( !fDebugStreamer ) {
1527 TDirectory *backup = gDirectory;
1528 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1529 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1534 Float_t y = ycenter;
1535 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1536 "caligroup="<<caligroup<<
1537 "detector="<<fDetectorPreviousTrack<<
1540 "npoints="<<nbclusters<<
1549 "padPosition="<<padPositions[ic]<<
1550 "padPosTracklet="<<padPosTracklet<<
1555 "signal1="<<signal1<<
1556 "signal2="<<signal2<<
1557 "signal3="<<signal3<<
1558 "signal4="<<signal4<<
1559 "signal5="<<signal5<<
1565 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1566 "caligroup="<<caligroup<<
1567 "detector="<<fDetectorPreviousTrack<<
1570 "npoints="<<nbclusters<<
1579 "padPosition="<<padPositions[ic]<<
1580 "padPosTracklet="<<padPosTracklet<<
1585 "signal1="<<signal1<<
1586 "signal2="<<signal2<<
1587 "signal3="<<signal3<<
1588 "signal4="<<signal4<<
1589 "signal5="<<signal5<<
1595 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1596 "caligroup="<<caligroup<<
1597 "detector="<<fDetectorPreviousTrack<<
1600 "npoints="<<nbclusters<<
1609 "padPosition="<<padPositions[ic]<<
1610 "padPosTracklet="<<padPosTracklet<<
1615 "signal1="<<signal1<<
1616 "signal2="<<signal2<<
1617 "signal3="<<signal3<<
1618 "signal4="<<signal4<<
1619 "signal5="<<signal5<<
1625 /////////////////////
1627 /////////////////////
1628 if(nbclusters < fNumberClusters) continue;
1629 if(nbclusters > fNumberClustersf) continue;
1630 if(nb3pc <= 5) continue;
1631 if((time >= 21) || (time < 7)) continue;
1632 if(TMath::Abs(qcl) < 80) continue;
1633 if( TMath::Abs(snp) > 1.) continue;
1636 ////////////////////////
1638 ///////////////////////
1640 if(TMath::Abs(dpad) < 1.5) {
1641 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1642 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1643 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1645 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1646 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1647 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1649 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1650 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1651 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1656 if(TMath::Abs(dpad) < 1.5) {
1657 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1658 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1660 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1661 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1662 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1664 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1665 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1666 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1669 } // second loop over clusters
1674 ///////////////////////////////////////////////////////////////////////////////////////
1675 // Pad row col stuff: see if masked or not
1676 ///////////////////////////////////////////////////////////////////////////////////////
1677 //_____________________________________________________________________________
1678 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1681 // See if we are not near a masked pad
1684 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1688 //_____________________________________________________________________________
1689 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1692 // See if we are not near a masked pad
1695 if (!IsPadOn(detector, col, row)) {
1696 fGoodTracklet = kFALSE;
1700 if (!IsPadOn(detector, col-1, row)) {
1701 fGoodTracklet = kFALSE;
1706 if (!IsPadOn(detector, col+1, row)) {
1707 fGoodTracklet = kFALSE;
1712 //_____________________________________________________________________________
1713 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1716 // Look in the choosen database if the pad is On.
1717 // If no the track will be "not good"
1720 // Get the parameter object
1721 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1723 AliInfo("Could not get calibDB");
1727 if (!cal->IsChamberInstalled(detector) ||
1728 cal->IsChamberMasked(detector) ||
1729 cal->IsPadMasked(detector,col,row)) {
1737 ///////////////////////////////////////////////////////////////////////////////////////
1738 // Calibration groups: calculate the number of groups, localise...
1739 ////////////////////////////////////////////////////////////////////////////////////////
1740 //_____________________________________________________________________________
1741 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1744 // Calculate the calibration group number for i
1747 // Row of the cluster and position in the pad groups
1749 if (fCalibraMode->GetNnZ(i) != 0) {
1750 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1754 // Col of the cluster and position in the pad groups
1756 if (fCalibraMode->GetNnRphi(i) != 0) {
1757 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1760 return posc*fCalibraMode->GetNfragZ(i)+posr;
1763 //____________________________________________________________________________________
1764 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1767 // Calculate the total number of calibration groups
1773 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1775 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1780 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1782 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1787 fCalibraMode->ModePadCalibration(2,i);
1788 fCalibraMode->ModePadFragmentation(0,2,0,i);
1789 fCalibraMode->SetDetChamb2(i);
1790 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1791 fCalibraMode->ModePadCalibration(0,i);
1792 fCalibraMode->ModePadFragmentation(0,0,0,i);
1793 fCalibraMode->SetDetChamb0(i);
1794 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1795 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1799 //_____________________________________________________________________________
1800 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1803 // Set the mode of calibration group in the z direction for the parameter i
1808 fCalibraMode->SetNz(i, Nz);
1811 AliInfo("You have to choose between 0 and 4");
1815 //_____________________________________________________________________________
1816 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1819 // Set the mode of calibration group in the rphi direction for the parameter i
1824 fCalibraMode->SetNrphi(i ,Nrphi);
1827 AliInfo("You have to choose between 0 and 6");
1832 //_____________________________________________________________________________
1833 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1836 // Set the mode of calibration group all together
1838 if(fVector2d == kTRUE) {
1839 AliInfo("Can not work with the vector method");
1842 fCalibraMode->SetAllTogether(i);
1846 //_____________________________________________________________________________
1847 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1850 // Set the mode of calibration group per supermodule
1852 if(fVector2d == kTRUE) {
1853 AliInfo("Can not work with the vector method");
1856 fCalibraMode->SetPerSuperModule(i);
1860 //____________Set the pad calibration variables for the detector_______________
1861 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1864 // For the detector calcul the first Xbins and set the number of row
1865 // and col pads per calibration groups, the number of calibration
1866 // groups in the detector.
1869 // first Xbins of the detector
1871 fCalibraMode->CalculXBins(detector,0);
1874 fCalibraMode->CalculXBins(detector,1);
1877 fCalibraMode->CalculXBins(detector,2);
1880 // fragmentation of idect
1881 for (Int_t i = 0; i < 3; i++) {
1882 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1883 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1884 , (Int_t) GetStack(detector)
1885 , (Int_t) GetSector(detector),i);
1891 //_____________________________________________________________________________
1892 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1895 // Should be between 0 and 6
1898 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1899 AliInfo("The number of groups must be between 0 and 6!");
1902 fNgroupprf = numberGroupsPRF;
1906 ///////////////////////////////////////////////////////////////////////////////////////////
1907 // Per tracklet: store or reset the info, fill the histos with the info
1908 //////////////////////////////////////////////////////////////////////////////////////////
1909 //_____________________________________________________________________________
1910 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1913 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1914 // Correct from the gain correction before
1917 // time bin of the cluster not corrected
1918 Int_t time = cl->GetPadTime();
1920 //Correct for the gain coefficient used in the database for reconstruction
1921 Float_t correctthegain = 1.0;
1922 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1923 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1924 Float_t correction = 1.0;
1925 Float_t normalisation = 6.67;
1926 // we divide with gain in AliTRDclusterizer::Transform...
1927 if( correctthegain > 0 ) normalisation /= correctthegain;
1930 // take dd/dl corrected from the angle of the track
1931 correction = dqdl / (normalisation);
1934 // Fill the fAmpTotal with the charge
1936 if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1939 // Fill the fPHPlace and value
1941 fPHPlace[time] = group[1];
1942 fPHValue[time] = correction;
1946 //____________Offine tracking in the AliTRDtracker_____________________________
1947 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1950 // Reset values per tracklet
1953 //Reset good tracklet
1954 fGoodTracklet = kTRUE;
1956 // Reset the fPHValue
1958 //Reset the fPHValue and fPHPlace
1959 for (Int_t k = 0; k < fTimeMax; k++) {
1965 // Reset the fAmpTotal where we put value
1967 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1972 //____________Offine tracking in the AliTRDtracker_____________________________
1973 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1976 // For the offline tracking or mcm tracklets
1977 // This function will be called in the functions UpdateHistogram...
1978 // to fill the info of a track for the relativ gain calibration
1981 Int_t nb = 0; // Nombre de zones traversees
1982 Int_t fd = -1; // Premiere zone non nulle
1983 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1985 if(nbclusters < fNumberClusters) return;
1986 if(nbclusters > fNumberClustersf) return;
1989 // Normalize with the number of clusters
1990 Double_t normalizeCst = fRelativeScale;
1991 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1993 // See if the track goes through different zones
1994 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1995 if (fAmpTotal[k] > 0.0) {
1996 totalcharge += fAmpTotal[k];
2009 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2011 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2012 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2015 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2019 if ((fAmpTotal[fd] > 0.0) &&
2020 (fAmpTotal[fd+1] > 0.0)) {
2021 // One of the two very big
2022 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2024 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2025 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2028 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2031 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2033 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2035 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2036 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2039 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2042 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2045 if (fCalibraMode->GetNfragZ(0) > 1) {
2046 if (fAmpTotal[fd] > 0.0) {
2047 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2048 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2049 // One of the two very big
2050 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2052 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2053 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2056 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2059 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2061 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2063 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2064 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2067 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2069 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2080 //____________Offine tracking in the AliTRDtracker_____________________________
2081 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2084 // For the offline tracking or mcm tracklets
2085 // This function will be called in the functions UpdateHistogram...
2086 // to fill the info of a track for the drift velocity calibration
2089 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2090 Int_t fd1 = -1; // Premiere zone non nulle
2091 Int_t fd2 = -1; // Deuxieme zone non nulle
2092 Int_t k1 = -1; // Debut de la premiere zone
2093 Int_t k2 = -1; // Debut de la seconde zone
2094 Int_t nbclusters = 0; // number of clusters
2098 // See if the track goes through different zones
2099 for (Int_t k = 0; k < fTimeMax; k++) {
2100 if (fPHValue[k] > 0.0) {
2106 if (fPHPlace[k] != fd1) {
2112 if (fPHPlace[k] != fd2) {
2119 // See if noise before and after
2120 if(fMaxCluster > 0) {
2121 if(fPHValue[0] > fMaxCluster) return;
2122 if(fTimeMax > fNbMaxCluster) {
2123 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2124 if(fPHValue[k] > fMaxCluster) return;
2129 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2131 if(nbclusters < fNumberClusters) return;
2132 if(nbclusters > fNumberClustersf) return;
2138 for (Int_t i = 0; i < fTimeMax; i++) {
2140 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2142 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2144 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2147 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2149 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2155 if ((fd1 == fd2+1) ||
2157 // One of the two fast all the think
2158 if (k2 > (k1+fDifference)) {
2159 //we choose to fill the fd1 with all the values
2161 for (Int_t i = 0; i < fTimeMax; i++) {
2163 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2165 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2169 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2171 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2176 if ((k2+fDifference) < fTimeMax) {
2177 //we choose to fill the fd2 with all the values
2179 for (Int_t i = 0; i < fTimeMax; i++) {
2181 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2183 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2187 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2189 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2195 // Two zones voisines sinon rien!
2196 if (fCalibraMode->GetNfragZ(1) > 1) {
2198 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2199 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2200 // One of the two fast all the think
2201 if (k2 > (k1+fDifference)) {
2202 //we choose to fill the fd1 with all the values
2204 for (Int_t i = 0; i < fTimeMax; i++) {
2206 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2208 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2212 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2214 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2219 if ((k2+fDifference) < fTimeMax) {
2220 //we choose to fill the fd2 with all the values
2222 for (Int_t i = 0; i < fTimeMax; i++) {
2224 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2226 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2230 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2232 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2239 // Two zones voisines sinon rien!
2241 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2242 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2243 // One of the two fast all the think
2244 if (k2 > (k1 + fDifference)) {
2245 //we choose to fill the fd1 with all the values
2247 for (Int_t i = 0; i < fTimeMax; i++) {
2249 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2251 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2255 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2257 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2262 if ((k2+fDifference) < fTimeMax) {
2263 //we choose to fill the fd2 with all the values
2265 for (Int_t i = 0; i < fTimeMax; i++) {
2267 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2269 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2273 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2275 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2287 //////////////////////////////////////////////////////////////////////////////////////////
2288 // DAQ process functions
2289 /////////////////////////////////////////////////////////////////////////////////////////
2290 //_____________________________________________________________________
2291 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2294 // Event Processing loop - AliTRDrawStreamBase
2295 // TestBeam 2007 version
2296 // 0 timebin problem
2300 // Same algorithm as TestBeam but different reader
2303 Int_t withInput = 1;
2305 Double_t phvalue[16][144][36];
2306 for(Int_t k = 0; k < 36; k++){
2307 for(Int_t j = 0; j < 16; j++){
2308 for(Int_t c = 0; c < 144; c++){
2309 phvalue[j][c][k] = 0.0;
2314 fDetectorPreviousTrack = -1;
2317 Int_t nbtimebin = 0;
2318 Int_t baseline = 10;
2325 while (rawStream->Next()) {
2327 Int_t idetector = rawStream->GetDet(); // current detector
2328 Int_t imcm = rawStream->GetMCM(); // current MCM
2329 Int_t irob = rawStream->GetROB(); // current ROB
2331 //printf("Detector %d\n",idetector);
2333 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2336 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2340 for(Int_t k = 0; k < 36; k++){
2341 for(Int_t j = 0; j < 16; j++){
2342 for(Int_t c = 0; c < 144; c++){
2343 phvalue[j][c][k] = 0.0;
2349 fDetectorPreviousTrack = idetector;
2350 fMCMPrevious = imcm;
2351 fROBPrevious = irob;
2353 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2354 if(nbtimebin == 0) return 0;
2355 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2356 fTimeMax = nbtimebin;
2358 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2359 fNumberClustersf = fTimeMax;
2360 fNumberClusters = (Int_t)(0.6*fTimeMax);
2363 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2364 Int_t col = rawStream->GetCol();
2365 Int_t row = rawStream->GetRow();
2368 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2371 for(Int_t itime = 0; itime < nbtimebin; itime++){
2372 phvalue[row][col][itime] = signal[itime]-baseline;
2376 // fill the last one
2377 if(fDetectorPreviousTrack != -1){
2380 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2383 for(Int_t k = 0; k < 36; k++){
2384 for(Int_t j = 0; j < 16; j++){
2385 for(Int_t c = 0; c < 144; c++){
2386 phvalue[j][c][k] = 0.0;
2395 while (rawStream->Next()) {
2397 Int_t idetector = rawStream->GetDet(); // current detector
2398 Int_t imcm = rawStream->GetMCM(); // current MCM
2399 Int_t irob = rawStream->GetROB(); // current ROB
2401 //printf("Detector %d\n",idetector);
2403 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2406 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2409 for(Int_t k = 0; k < 36; k++){
2410 for(Int_t j = 0; j < 16; j++){
2411 for(Int_t c = 0; c < 144; c++){
2412 phvalue[j][c][k] = 0.0;
2418 fDetectorPreviousTrack = idetector;
2419 fMCMPrevious = imcm;
2420 fROBPrevious = irob;
2422 //baseline = rawStream->GetCommonAdditive(); // common baseline
2424 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2425 fNumberClustersf = fTimeMax;
2426 fNumberClusters = (Int_t)(0.6*fTimeMax);
2427 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2428 Int_t col = rawStream->GetCol();
2429 Int_t row = rawStream->GetRow();
2432 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2434 for(Int_t itime = 0; itime < fTimeMax; itime++){
2435 phvalue[row][col][itime] = signal[itime]-baseline;
2439 // fill the last one
2440 if(fDetectorPreviousTrack != -1){
2443 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2446 for(Int_t k = 0; k < 36; k++){
2447 for(Int_t j = 0; j < 16; j++){
2448 for(Int_t c = 0; c < 144; c++){
2449 phvalue[j][c][k] = 0.0;
2459 //_____________________________________________________________________
2460 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2463 // Event processing loop - AliRawReader
2464 // Testbeam 2007 version
2467 AliTRDrawStreamBase rawStream(rawReader);
2469 rawReader->Select("TRD");
2471 return ProcessEventDAQ(&rawStream, nocheck);
2474 //_________________________________________________________________________
2475 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2477 eventHeaderStruct *event,
2480 eventHeaderStruct* /*event*/,
2487 // process date event
2488 // Testbeam 2007 version
2491 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2492 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2496 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2501 //////////////////////////////////////////////////////////////////////////////
2502 // Routine inside the DAQ process
2503 /////////////////////////////////////////////////////////////////////////////
2504 //_______________________________________________________________________
2505 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2508 // Look for the maximum by collapsing over the time
2509 // Sum over four pad col and two pad row
2515 Int_t idect = fDetectorPreviousTrack;
2516 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2518 for(Int_t tb = 0; tb < 36; tb++){
2522 //fGoodTracklet = kTRUE;
2523 //fDetectorPreviousTrack = 0;
2526 ///////////////////////////
2528 /////////////////////////
2532 Double_t integralMax = -1;
2534 for (Int_t ir = 1; ir <= 15; ir++)
2536 for (Int_t ic = 2; ic <= 142; ic++)
2538 Double_t integral = 0;
2539 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2541 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2543 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2544 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2547 for(Int_t tb = 0; tb< fTimeMax; tb++){
2548 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2554 if (integralMax < integral)
2558 integralMax = integral;
2559 } // check max integral
2563 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2565 if((imaxRow == 0) || (imaxCol == 0)) {
2569 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2570 //if(!fGoodTracklet) used = 1;;
2572 // /////////////////////////////////////////////////////
2573 // sum ober 2 row and 4 pad cols for each time bins
2574 // ////////////////////////////////////////////////////
2577 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2579 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2581 for(Int_t it = 0; it < fTimeMax; it++){
2582 sum[it] += phvalue[ir][ic][it];
2588 Double_t sumcharge = 0.0;
2589 for(Int_t it = 0; it < fTimeMax; it++){
2590 sumcharge += sum[it];
2591 if(sum[it] > 20.0) nbcl++;
2595 /////////////////////////////////////////////////////////
2597 ////////////////////////////////////////////////////////
2598 if(fDebugLevel > 0){
2599 if ( !fDebugStreamer ) {
2601 TDirectory *backup = gDirectory;
2602 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2603 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2606 Double_t amph0 = sum[0];
2607 Double_t amphlast = sum[fTimeMax-1];
2608 Double_t rms = TMath::RMS(fTimeMax,sum);
2609 Int_t goodtracklet = (Int_t) fGoodTracklet;
2610 for(Int_t it = 0; it < fTimeMax; it++){
2611 Double_t clustera = sum[it];
2613 (* fDebugStreamer) << "FillDAQa"<<
2614 "ampTotal="<<sumcharge<<
2617 "detector="<<idect<<
2619 "amphlast="<<amphlast<<
2620 "goodtracklet="<<goodtracklet<<
2621 "clustera="<<clustera<<
2628 ////////////////////////////////////////////////////////
2630 ///////////////////////////////////////////////////////
2631 if(sum[0] > 100.0) used = 1;
2632 if(nbcl < fNumberClusters) used = 1;
2633 if(nbcl > fNumberClustersf) used = 1;
2635 //if(fDetectorPreviousTrack == 15){
2636 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2638 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2640 for(Int_t it = 0; it < fTimeMax; it++){
2641 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2643 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2648 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2650 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2657 //____________Online trackling in AliTRDtrigger________________________________
2658 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2662 // Fill a simple average pulse height
2666 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2672 //____________Write_____________________________________________________
2673 //_____________________________________________________________________
2674 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2677 // Write infos to file
2681 if ( fDebugStreamer ) {
2682 delete fDebugStreamer;
2683 fDebugStreamer = 0x0;
2686 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2691 ,fNumberUsedPh[1]));
2693 TDirectory *backup = gDirectory;
2699 option = "recreate";
2701 TFile f(filename,option.Data());
2703 TStopwatch stopwatch;
2706 f.WriteTObject(fCalibraVector);
2711 f.WriteTObject(fCH2d);
2716 f.WriteTObject(fPH2d);
2721 f.WriteTObject(fPRF2d);
2724 if(fLinearFitterOn){
2725 AnalyseLinearFitter();
2726 f.WriteTObject(fLinearVdriftFit);
2731 if ( backup ) backup->cd();
2733 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2734 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2736 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2738 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2739 //___________________________________________probe the histos__________________________________________________
2740 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2743 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2744 // debug mode with 2 for TH2I and 3 for TProfile2D
2745 // It gives a pointer to a Double_t[7] with the info following...
2746 // [0] : number of calibration groups with entries
2747 // [1] : minimal number of entries found
2748 // [2] : calibration group number of the min
2749 // [3] : maximal number of entries found
2750 // [4] : calibration group number of the max
2751 // [5] : mean number of entries found
2752 // [6] : mean relative error
2755 Double_t *info = new Double_t[7];
2757 // Number of Xbins (detectors or groups of pads)
2758 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2759 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2762 Double_t nbwe = 0; //number of calibration groups with entries
2763 Double_t minentries = 0; //minimal number of entries found
2764 Double_t maxentries = 0; //maximal number of entries found
2765 Double_t placemin = 0; //calibration group number of the min
2766 Double_t placemax = -1; //calibration group number of the max
2767 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2768 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2770 Double_t counter = 0;
2773 TH1F *nbEntries = 0x0;//distribution of the number of entries
2774 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2775 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2777 // Beginning of the loop over the calibration groups
2778 for (Int_t idect = 0; idect < nbins; idect++) {
2780 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2781 projch->SetDirectory(0);
2783 // Number of entries for this calibration group
2784 Double_t nentries = 0.0;
2786 for (Int_t k = 0; k < nxbins; k++) {
2787 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2791 for (Int_t k = 0; k < nxbins; k++) {
2792 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2793 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2794 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2802 if((!((Bool_t)nbEntries)) && (nentries > 0)){
2803 nbEntries = new TH1F("Number of entries","Number of entries"
2804 ,100,(Int_t)nentries/2,nentries*2);
2805 nbEntries->SetDirectory(0);
2806 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
2808 nbEntriesPerGroup->SetDirectory(0);
2809 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
2810 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2811 nbEntriesPerSp->SetDirectory(0);
2814 if(nentries > 0) nbEntries->Fill(nentries);
2815 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2816 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2821 if(nentries > maxentries){
2822 maxentries = nentries;
2826 minentries = nentries;
2828 if(nentries < minentries){
2829 minentries = nentries;
2835 meanstats += nentries;
2837 }//calibration groups loop
2839 if(nbwe > 0) meanstats /= nbwe;
2840 if(counter > 0) meanrelativerror /= counter;
2842 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2843 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2844 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2845 AliInfo(Form("The mean number of entries is %f",meanstats));
2846 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2849 info[1] = minentries;
2851 info[3] = maxentries;
2853 info[5] = meanstats;
2854 info[6] = meanrelativerror;
2857 gStyle->SetPalette(1);
2858 gStyle->SetOptStat(1111);
2859 gStyle->SetPadBorderMode(0);
2860 gStyle->SetCanvasColor(10);
2861 gStyle->SetPadLeftMargin(0.13);
2862 gStyle->SetPadRightMargin(0.01);
2863 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2866 nbEntries->Draw("");
2868 nbEntriesPerSp->SetStats(0);
2869 nbEntriesPerSp->Draw("");
2870 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2872 nbEntriesPerGroup->SetStats(0);
2873 nbEntriesPerGroup->Draw("");
2879 //____________________________________________________________________________
2880 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2883 // Return a Int_t[4] with:
2884 // 0 Mean number of entries
2885 // 1 median of number of entries
2886 // 2 rms of number of entries
2887 // 3 number of group with entries
2890 Double_t *stat = new Double_t[4];
2893 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2894 Double_t *weight = new Double_t[nbofgroups];
2895 Int_t *nonul = new Int_t[nbofgroups];
2897 for(Int_t k = 0; k < nbofgroups; k++){
2898 if(fEntriesCH[k] > 0) {
2900 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2903 else weight[k] = 0.0;
2905 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2906 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2907 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2912 //____________________________________________________________________________
2913 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2916 // Return a Int_t[4] with:
2917 // 0 Mean number of entries
2918 // 1 median of number of entries
2919 // 2 rms of number of entries
2920 // 3 number of group with entries
2923 Double_t *stat = new Double_t[4];
2926 Int_t nbofgroups = 540;
2927 Double_t *weight = new Double_t[nbofgroups];
2928 Int_t *nonul = new Int_t[nbofgroups];
2930 for(Int_t k = 0; k < nbofgroups; k++){
2931 if(fEntriesLinearFitter[k] > 0) {
2933 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2936 else weight[k] = 0.0;
2938 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2939 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2940 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2945 //////////////////////////////////////////////////////////////////////////////////////
2947 //////////////////////////////////////////////////////////////////////////////////////
2948 //_____________________________________________________________________________
2949 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2952 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2953 // If fNgroupprf is zero then no binning in tnp
2957 name += fCalibraMode->GetNz(2);
2959 name += fCalibraMode->GetNrphi(2);
2963 if(fNgroupprf != 0){
2965 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2966 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2967 fPRF2d->SetYTitle("Det/pad groups");
2968 fPRF2d->SetXTitle("Position x/W [pad width units]");
2969 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2970 fPRF2d->SetStats(0);
2973 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2974 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2975 fPRF2d->SetYTitle("Det/pad groups");
2976 fPRF2d->SetXTitle("Position x/W [pad width units]");
2977 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2978 fPRF2d->SetStats(0);
2983 //_____________________________________________________________________________
2984 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2987 // Create the 2D histos
2991 name += fCalibraMode->GetNz(1);
2993 name += fCalibraMode->GetNrphi(1);
2995 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2996 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2998 fPH2d->SetYTitle("Det/pad groups");
2999 fPH2d->SetXTitle("time [#mus]");
3000 fPH2d->SetZTitle("<PH> [a.u.]");
3004 //_____________________________________________________________________________
3005 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3008 // Create the 2D histos
3012 name += fCalibraMode->GetNz(0);
3014 name += fCalibraMode->GetNrphi(0);
3016 fCH2d = new TH2I("CH2d",(const Char_t *) name
3017 ,fNumberBinCharge,0,300,nn,0,nn);
3018 fCH2d->SetYTitle("Det/pad groups");
3019 fCH2d->SetXTitle("charge deposit [a.u]");
3020 fCH2d->SetZTitle("counts");
3025 //////////////////////////////////////////////////////////////////////////////////
3026 // Set relative scale
3027 /////////////////////////////////////////////////////////////////////////////////
3028 //_____________________________________________________________________________
3029 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3032 // Set the factor that will divide the deposited charge
3033 // to fit in the histo range [0,300]
3036 if (RelativeScale > 0.0) {
3037 fRelativeScale = RelativeScale;
3040 AliInfo("RelativeScale must be strict positif!");
3044 //////////////////////////////////////////////////////////////////////////////////
3045 // Quick way to fill a histo
3046 //////////////////////////////////////////////////////////////////////////////////
3047 //_____________________________________________________________________
3048 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3051 // FillCH2d: Marian style
3054 //skip simply the value out of range
3055 if((y>=300.0) || (y<0.0)) return;
3057 //Calcul the y place
3058 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3059 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3062 fCH2d->GetArray()[place]++;
3066 //////////////////////////////////////////////////////////////////////////////////
3067 // Geometrical functions
3068 ///////////////////////////////////////////////////////////////////////////////////
3069 //_____________________________________________________________________________
3070 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3073 // Reconstruct the layer number from the detector number
3076 return ((Int_t) (d % 6));
3080 //_____________________________________________________________________________
3081 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3084 // Reconstruct the stack number from the detector number
3086 const Int_t kNlayer = 6;
3088 return ((Int_t) (d % 30) / kNlayer);
3092 //_____________________________________________________________________________
3093 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3096 // Reconstruct the sector number from the detector number
3100 return ((Int_t) (d / fg));
3103 ///////////////////////////////////////////////////////////////////////////////////
3104 // Getter functions for DAQ of the CH2d and the PH2d
3105 //////////////////////////////////////////////////////////////////////////////////
3106 //_____________________________________________________________________
3107 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3110 // return pointer to fPH2d TProfile2D
3111 // create a new TProfile2D if it doesn't exist allready
3117 fTimeMax = nbtimebin;
3118 fSf = samplefrequency;
3124 //_____________________________________________________________________
3125 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3128 // return pointer to fCH2d TH2I
3129 // create a new TH2I if it doesn't exist allready
3138 ////////////////////////////////////////////////////////////////////////////////////////////
3139 // Drift velocity calibration
3140 ///////////////////////////////////////////////////////////////////////////////////////////
3141 //_____________________________________________________________________
3142 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3145 // return pointer to TLinearFitter Calibration
3146 // if force is true create a new TLinearFitter if it doesn't exist allready
3149 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3150 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3153 // if we are forced and TLinearFitter doesn't yet exist create it
3155 // new TLinearFitter
3156 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3157 fLinearFitterArray.AddAt(linearfitter,detector);
3158 return linearfitter;
3161 //____________________________________________________________________________
3162 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3165 // Analyse array of linear fitter because can not be written
3166 // Store two arrays: one with the param the other one with the error param + number of entries
3169 for(Int_t k = 0; k < 540; k++){
3170 TLinearFitter *linearfitter = GetLinearFitter(k);
3171 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3172 TVectorD *par = new TVectorD(2);
3173 TVectorD pare = TVectorD(2);
3174 TVectorD *parE = new TVectorD(3);
3175 linearfitter->Eval();
3176 linearfitter->GetParameters(*par);
3177 linearfitter->GetErrors(pare);
3178 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3179 (*parE)[0] = pare[0]*ppointError;
3180 (*parE)[1] = pare[1]*ppointError;
3181 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3182 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3183 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);