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()
127 ,fLinearFitterOn(kFALSE)
128 ,fLinearFitterDebugOn(kFALSE)
130 ,fThresholdClusterPRF2(15.0)
131 ,fLimitChargeIntegration(kFALSE)
132 ,fFillWithZero(kFALSE)
133 ,fNormalizeNbOfCluster(kFALSE)
136 ,fCalibraMode(new AliTRDCalibraMode())
139 ,fDetectorPreviousTrack(-1)
143 ,fNumberClustersf(30)
149 ,fNumberBinCharge(50)
155 ,fGoodTracklet(kTRUE)
156 ,fLinearFitterTracklet(0x0)
158 ,fEntriesLinearFitter(0x0)
163 ,fLinearFitterArray(540)
164 ,fLinearVdriftFit(0x0)
169 // Default constructor
173 // Init some default values
176 fNumberUsedCh[0] = 0;
177 fNumberUsedCh[1] = 0;
178 fNumberUsedPh[0] = 0;
179 fNumberUsedPh[1] = 0;
181 fGeo = new AliTRDgeometry();
185 //______________________________________________________________________________________
186 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
192 ,fPRF2dOn(c.fPRF2dOn)
193 ,fHisto2d(c.fHisto2d)
194 ,fVector2d(c.fVector2d)
195 ,fLinearFitterOn(c.fLinearFitterOn)
196 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
197 ,fRelativeScale(c.fRelativeScale)
198 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
199 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
200 ,fFillWithZero(c.fFillWithZero)
201 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
202 ,fMaxCluster(c.fMaxCluster)
203 ,fNbMaxCluster(c.fNbMaxCluster)
206 ,fDebugLevel(c.fDebugLevel)
207 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
208 ,fMCMPrevious(c.fMCMPrevious)
209 ,fROBPrevious(c.fROBPrevious)
210 ,fNumberClusters(c.fNumberClusters)
211 ,fNumberClustersf(c.fNumberClustersf)
212 ,fProcent(c.fProcent)
213 ,fDifference(c.fDifference)
214 ,fNumberTrack(c.fNumberTrack)
215 ,fTimeMax(c.fTimeMax)
217 ,fNumberBinCharge(c.fNumberBinCharge)
218 ,fNumberBinPRF(c.fNumberBinPRF)
219 ,fNgroupprf(c.fNgroupprf)
223 ,fGoodTracklet(c.fGoodTracklet)
224 ,fLinearFitterTracklet(0x0)
226 ,fEntriesLinearFitter(0x0)
231 ,fLinearFitterArray(540)
232 ,fLinearVdriftFit(0x0)
239 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
240 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
242 fPH2d = new TProfile2D(*c.fPH2d);
243 fPH2d->SetDirectory(0);
246 fPRF2d = new TProfile2D(*c.fPRF2d);
247 fPRF2d->SetDirectory(0);
250 fCH2d = new TH2I(*c.fCH2d);
251 fCH2d->SetDirectory(0);
253 if(c.fLinearVdriftFit){
254 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
257 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
258 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
263 fGeo = new AliTRDgeometry();
266 //____________________________________________________________________________________
267 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
270 // AliTRDCalibraFillHisto destructor
274 if ( fDebugStreamer ) delete fDebugStreamer;
276 if ( fCalDetGain ) delete fCalDetGain;
277 if ( fCalROCGain ) delete fCalROCGain;
279 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
283 delete [] fEntriesCH;
284 delete [] fEntriesLinearFitter;
287 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
288 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
296 //_____________________________________________________________________________
297 void AliTRDCalibraFillHisto::Destroy()
308 //_____________________________________________________________________________
309 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
312 // Delete DebugStreamer
315 if ( fDebugStreamer ) delete fDebugStreamer;
316 fDebugStreamer = 0x0;
319 //_____________________________________________________________________________
320 void AliTRDCalibraFillHisto::ClearHistos()
340 //////////////////////////////////////////////////////////////////////////////////
341 // calibration with AliTRDtrackV1: Init, Update
342 //////////////////////////////////////////////////////////////////////////////////
343 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
344 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
347 // Init the histograms and stuff to be filled
352 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
354 AliInfo("Could not get calibDB");
357 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
359 AliInfo("Could not get CommonParam");
364 fTimeMax = cal->GetNumberOfTimeBins();
365 fSf = parCom->GetSamplingFrequency();
366 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
367 else fRelativeScale = 1.18;
368 fNumberClustersf = fTimeMax;
369 fNumberClusters = (Int_t)(0.5*fTimeMax);
371 // Init linear fitter
372 if(!fLinearFitterTracklet) {
373 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
374 fLinearFitterTracklet->StoreData(kTRUE);
377 //calib object from database used for reconstruction
379 fCalDetGain->~AliTRDCalDet();
380 new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
381 }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
383 // Calcul Xbins Chambd0, Chamb2
384 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
385 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
386 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
388 // If vector method On initialised all the stuff
390 fCalibraVector = new AliTRDCalibraVector();
391 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
392 fCalibraVector->SetTimeMax(fTimeMax);
393 if(fNgroupprf != 0) {
394 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
395 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
398 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
399 fCalibraVector->SetPRFRange(1.5);
401 for(Int_t k = 0; k < 3; k++){
402 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
403 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
405 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
406 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
407 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
408 fCalibraVector->SetNbGroupPRF(fNgroupprf);
411 // Create the 2D histos corresponding to the pad groupCalibration mode
414 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
415 ,fCalibraMode->GetNz(0)
416 ,fCalibraMode->GetNrphi(0)));
418 // Create the 2D histo
423 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
424 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
428 fEntriesCH = new Int_t[ntotal0];
429 for(Int_t k = 0; k < ntotal0; k++){
436 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
437 ,fCalibraMode->GetNz(1)
438 ,fCalibraMode->GetNrphi(1)));
440 // Create the 2D histo
445 fPHPlace = new Short_t[fTimeMax];
446 for (Int_t k = 0; k < fTimeMax; k++) {
449 fPHValue = new Float_t[fTimeMax];
450 for (Int_t k = 0; k < fTimeMax; k++) {
454 if (fLinearFitterOn) {
455 //fLinearFitterArray.Expand(540);
456 fLinearFitterArray.SetName("ArrayLinearFitters");
457 fEntriesLinearFitter = new Int_t[540];
458 for(Int_t k = 0; k < 540; k++){
459 fEntriesLinearFitter[k] = 0;
461 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
466 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
467 ,fCalibraMode->GetNz(2)
468 ,fCalibraMode->GetNrphi(2)));
469 // Create the 2D histo
471 CreatePRF2d(ntotal2);
478 //____________Offline tracking in the AliTRDtracker____________________________
479 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
482 // Use AliTRDtrack for the calibration
486 AliTRDcluster *cl = 0x0; // pointeur to cluster
487 Int_t index0 = 0; // index of the first cluster in the current chamber
488 Int_t index1 = 0; // index of the current cluster in the current chamber
490 //////////////////////////////
491 // loop over the clusters
492 ///////////////////////////////
493 while((cl = t->GetCluster(index1))){
495 /////////////////////////////////////////////////////////////////////////
496 // Fill the infos for the previous clusters if not the same detector
497 ////////////////////////////////////////////////////////////////////////
498 Int_t detector = cl->GetDetector();
499 if ((detector != fDetectorPreviousTrack) &&
500 (index0 != index1)) {
504 //If the same track, then look if the previous detector is in
505 //the same plane, if yes: not a good track
506 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
510 // Fill only if the track doesn't touch a masked pad or doesn't
513 // drift velocity unables to cut bad tracklets
514 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
518 FillTheInfoOfTheTrackCH(index1-index0);
523 FillTheInfoOfTheTrackPH();
526 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
529 } // if a good tracklet
532 ResetfVariablestracklet();
535 } // Fill at the end the charge
538 /////////////////////////////////////////////////////////////
539 // Localise and take the calib gain object for the detector
540 ////////////////////////////////////////////////////////////
541 if (detector != fDetectorPreviousTrack) {
543 //Localise the detector
544 LocalisationDetectorXbins(detector);
547 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
549 AliInfo("Could not get calibDB");
556 fCalROCGain->~AliTRDCalROC();
557 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
559 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
563 // Reset the detectbjobsor
564 fDetectorPreviousTrack = detector;
566 ///////////////////////////////////////
567 // Store the info of the cluster
568 ///////////////////////////////////////
569 Int_t row = cl->GetPadRow();
570 Int_t col = cl->GetPadCol();
571 CheckGoodTrackletV1(cl);
572 Int_t group[2] = {0,0};
573 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
574 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
575 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
579 } // while on clusters
581 ///////////////////////////
582 // Fill the last plane
583 //////////////////////////
584 if( index0 != index1 ){
590 // drift velocity unables to cut bad tracklets
591 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
595 FillTheInfoOfTheTrackCH(index1-index0);
600 FillTheInfoOfTheTrackPH();
603 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
605 } // if a good tracklet
610 ResetfVariablestracklet();
615 //____________Offline tracking in the AliTRDtracker____________________________
616 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
619 // Use AliTRDtrackV1 for the calibration
623 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
624 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
625 Bool_t newtr = kTRUE; // new track
628 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
630 AliInfo("Could not get calibDB");
634 ///////////////////////////
635 // loop over the tracklet
636 ///////////////////////////
637 for(Int_t itr = 0; itr < 6; itr++){
639 if(!(tracklet = t->GetTracklet(itr))) continue;
640 if(!tracklet->IsOK()) continue;
642 ResetfVariablestracklet();
644 //////////////////////////////////////////
645 // localisation of the tracklet and dqdl
646 //////////////////////////////////////////
647 Int_t layer = tracklet->GetPlane();
649 while(!(cl = tracklet->GetClusters(ic++))) continue;
650 Int_t detector = cl->GetDetector();
651 if (detector != fDetectorPreviousTrack) {
652 // if not a new track
654 // don't use the rest of this track if in the same plane
655 if (layer == GetLayer(fDetectorPreviousTrack)) {
656 //printf("bad tracklet, same layer for detector %d\n",detector);
660 //Localise the detector bin
661 LocalisationDetectorXbins(detector);
665 fCalROCGain->~AliTRDCalROC();
666 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
668 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
671 fDetectorPreviousTrack = detector;
675 ////////////////////////////
676 // loop over the clusters
677 ////////////////////////////
678 Int_t nbclusters = 0;
679 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
680 if(!(cl = tracklet->GetClusters(jc))) continue;
683 // Store the info bis of the tracklet
684 Int_t row = cl->GetPadRow();
685 Int_t col = cl->GetPadCol();
686 CheckGoodTrackletV1(cl);
687 Int_t group[2] = {0,0};
688 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
689 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
690 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col);
693 ////////////////////////////////////////
694 // Fill the stuffs if a good tracklet
695 ////////////////////////////////////////
698 // drift velocity unables to cut bad tracklets
699 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
701 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
705 FillTheInfoOfTheTrackCH(nbclusters);
710 FillTheInfoOfTheTrackPH();
713 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
715 } // if a good tracklet
721 ///////////////////////////////////////////////////////////////////////////////////
722 // Routine inside the update with AliTRDtrack
723 ///////////////////////////////////////////////////////////////////////////////////
724 //____________Offine tracking in the AliTRDtracker_____________________________
725 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
728 // Drift velocity calibration:
729 // Fit the clusters with a straight line
730 // From the slope find the drift velocity
733 //Number of points: if less than 3 return kFALSE
734 Int_t npoints = index1-index0;
735 if(npoints <= 2) return kFALSE;
740 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
741 // parameters of the track
742 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
743 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
744 Double_t tnp = 0.0; // tan angle in the plan xy track
745 if( TMath::Abs(snp) < 1.){
746 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
748 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
749 // tilting pad and cross row
750 Int_t crossrow = 0; // if it crosses a pad row
751 Int_t rowp = -1; // if it crosses a pad row
752 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
753 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
754 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
756 fLinearFitterTracklet->ClearPoints();
757 Double_t dydt = 0.0; // dydt tracklet after straight line fit
758 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
759 Double_t pointError = 0.0; // error after straight line fit
760 Int_t nbli = 0; // number linear fitter points
762 //////////////////////////////
763 // loop over clusters
764 ////////////////////////////
765 for(Int_t k = 0; k < npoints; k++){
767 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
768 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
769 Double_t ycluster = cl->GetY();
770 Int_t time = cl->GetPadTime();
771 Double_t timeis = time/fSf;
772 //Double_t q = cl->GetQ();
773 //See if cross two pad rows
774 Int_t row = cl->GetPadRow();
776 if(row != rowp) crossrow = 1;
778 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
783 //////////////////////////////
785 //////////////////////////////
787 fLinearFitterTracklet->ClearPoints();
791 fLinearFitterTracklet->Eval();
792 fLinearFitterTracklet->GetParameters(pars);
793 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
794 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
796 fLinearFitterTracklet->ClearPoints();
798 /////////////////////////////
800 ////////////////////////////
802 if ( !fDebugStreamer ) {
804 TDirectory *backup = gDirectory;
805 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
806 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
810 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
811 //"snpright="<<snpright<<
812 "npoints="<<npoints<<
814 "detector="<<detector<<
821 "crossrow="<<crossrow<<
822 "errorpar="<<errorpar<<
823 "pointError="<<pointError<<
827 Int_t nbclusters = index1-index0;
828 Int_t layer = GetLayer(fDetectorPreviousTrack);
830 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
831 //"snpright="<<snpright<<
832 "nbclusters="<<nbclusters<<
833 "detector="<<fDetectorPreviousTrack<<
839 ///////////////////////////
841 ///////////////////////////
842 if(npoints < fNumberClusters) return kFALSE;
843 if(npoints > fNumberClustersf) return kFALSE;
844 if(pointError >= 0.1) return kFALSE;
845 if(crossrow == 1) return kFALSE;
847 ////////////////////////////
849 ////////////////////////////
851 //Add to the linear fitter of the detector
852 if( TMath::Abs(snp) < 1.){
853 Double_t x = tnp-dzdx*tnt;
854 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
855 if(fLinearFitterDebugOn) {
856 fLinearVdriftFit->Update(detector,x,pars[1]);
858 fEntriesLinearFitter[detector]++;
864 //____________Offine tracking in the AliTRDtracker_____________________________
865 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
868 // Drift velocity calibration:
869 // Fit the clusters with a straight line
870 // From the slope find the drift velocity
873 ////////////////////////////////////////////////
874 //Number of points: if less than 3 return kFALSE
875 /////////////////////////////////////////////////
876 if(nbclusters <= 2) return kFALSE;
881 // results of the linear fit
882 Double_t dydt = 0.0; // dydt tracklet after straight line fit
883 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
884 Double_t pointError = 0.0; // error after straight line fit
885 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
886 Int_t crossrow = 0; // if it crosses a pad row
887 Int_t rowp = -1; // if it crosses a pad row
888 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
889 fLinearFitterTracklet->ClearPoints();
892 ///////////////////////////////////////////
893 // Take the parameters of the track
894 //////////////////////////////////////////
895 // take now the snp, tnp and tgl from the track
896 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
897 Double_t tnp = 0.0; // dy/dx at the end of the chamber
898 if( TMath::Abs(snp) < 1.){
899 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
901 Double_t tgl = tracklet->GetTgl(); // dz/dl
902 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
904 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
905 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
906 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
907 // at the end with correction due to linear fit
908 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
909 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
912 ////////////////////////////
913 // loop over the clusters
914 ////////////////////////////
916 AliTRDcluster *cl = 0x0;
917 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
918 if(!(cl = tracklet->GetClusters(ic))) continue;
919 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
921 Double_t ycluster = cl->GetY();
922 Int_t time = cl->GetPadTime();
923 Double_t timeis = time/fSf;
924 //See if cross two pad rows
925 Int_t row = cl->GetPadRow();
926 if(rowp==-1) rowp = row;
927 if(row != rowp) crossrow = 1;
929 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
932 //////////////////////////////
933 // Check no shared clusters
934 //////////////////////////////
935 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
936 if((cl = tracklet->GetClusters(icc))) crossrow = 1;
940 ////////////////////////////////////
941 // Do the straight line fit now
942 ///////////////////////////////////
944 fLinearFitterTracklet->ClearPoints();
948 fLinearFitterTracklet->Eval();
949 fLinearFitterTracklet->GetParameters(pars);
950 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
951 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
953 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
954 fLinearFitterTracklet->ClearPoints();
956 ////////////////////////////////
958 ///////////////////////////////
962 if ( !fDebugStreamer ) {
964 TDirectory *backup = gDirectory;
965 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
966 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
970 Int_t layer = GetLayer(fDetectorPreviousTrack);
972 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
973 //"snpright="<<snpright<<
975 "nbclusters="<<nbclusters<<
976 "detector="<<fDetectorPreviousTrack<<
984 "crossrow="<<crossrow<<
985 "errorpar="<<errorpar<<
986 "pointError="<<pointError<<
991 /////////////////////////
993 ////////////////////////
995 if(nbclusters < fNumberClusters) return kFALSE;
996 if(nbclusters > fNumberClustersf) return kFALSE;
997 if(pointError >= 0.3) return kFALSE;
998 if(crossrow == 1) return kFALSE;
1000 ///////////////////////
1002 //////////////////////
1004 if(fLinearFitterOn){
1005 //Add to the linear fitter of the detector
1006 if( TMath::Abs(snp) < 1.){
1007 Double_t x = tnp-dzdx*tnt;
1008 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1009 if(fLinearFitterDebugOn) {
1010 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1012 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1018 //____________Offine tracking in the AliTRDtracker_____________________________
1019 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1022 // PRF width calibration
1023 // Assume a Gaussian shape: determinate the position of the three pad clusters
1024 // Fit with a straight line
1025 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1026 // Fill the PRF as function of angle of the track
1031 //////////////////////////
1033 /////////////////////////
1034 Int_t npoints = index1-index0; // number of total points
1035 Int_t nb3pc = 0; // number of three pads clusters used for fit
1036 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1037 // To see the difference due to the fit
1038 Double_t *padPositions;
1039 padPositions = new Double_t[npoints];
1040 for(Int_t k = 0; k < npoints; k++){
1041 padPositions[k] = 0.0;
1043 // Take the tgl and snp with the track t now
1044 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1045 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1046 Float_t dzdx = 0.0; // dzdx
1048 if(TMath::Abs(snp) < 1.0){
1049 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1050 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1053 fLinearFitterTracklet->ClearPoints();
1055 ///////////////////////////
1056 // calcul the tnp group
1057 ///////////////////////////
1058 Bool_t echec = kFALSE;
1059 Double_t shift = 0.0;
1060 //Calculate the shift in x coresponding to this tnp
1061 if(fNgroupprf != 0.0){
1062 shift = -3.0*(fNgroupprf-1)-1.5;
1063 Double_t limithigh = -0.2*(fNgroupprf-1);
1064 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1066 while(tnp > limithigh){
1073 delete [] padPositions;
1077 //////////////////////
1079 /////////////////////
1080 for(Int_t k = 0; k < npoints; k++){
1082 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1084 Short_t *signals = cl->GetSignals();
1085 Double_t time = cl->GetPadTime();
1086 //Calculate x if possible
1087 Float_t xcenter = 0.0;
1088 Bool_t echec1 = kTRUE;
1089 if((time<=7) || (time>=21)) continue;
1090 // Center 3 balanced: position with the center of the pad
1091 if ((((Float_t) signals[3]) > 0.0) &&
1092 (((Float_t) signals[2]) > 0.0) &&
1093 (((Float_t) signals[4]) > 0.0)) {
1095 // Security if the denomiateur is 0
1096 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1097 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1098 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1099 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1100 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1106 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1108 //if no echec: calculate with the position of the pad
1109 // Position of the cluster
1110 Double_t padPosition = xcenter + cl->GetPadCol();
1111 padPositions[k] = padPosition;
1113 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1117 /////////////////////////////
1119 ////////////////////////////
1121 delete [] padPositions;
1122 fLinearFitterTracklet->ClearPoints();
1125 fLinearFitterTracklet->Eval();
1127 fLinearFitterTracklet->GetParameters(line);
1128 Float_t pointError = -1.0;
1129 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1130 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1132 fLinearFitterTracklet->ClearPoints();
1134 /////////////////////////////////////////////////////
1135 // Now fill the PRF: second loop over clusters
1136 /////////////////////////////////////////////////////
1137 for(Int_t k = 0; k < npoints; k++){
1139 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1140 Short_t *signals = cl->GetSignals(); // signal
1141 Double_t time = cl->GetPadTime(); // time bin
1142 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1143 Float_t padPos = cl->GetPadCol(); // middle pad
1144 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1145 Float_t ycenter = 0.0; // relative center charge
1146 Float_t ymin = 0.0; // relative left charge
1147 Float_t ymax = 0.0; // relative right charge
1149 //Requiere simply two pads clusters at least
1150 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1151 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1152 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1153 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1154 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1155 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1159 Int_t rowcl = cl->GetPadRow(); // row of cluster
1160 Int_t colcl = cl->GetPadCol(); // col of cluster
1161 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1162 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1163 Float_t xcl = cl->GetY(); // y cluster
1164 Float_t qcl = cl->GetQ(); // charge cluster
1165 Int_t layer = GetLayer(detector); // layer
1166 Int_t stack = GetStack(detector); // stack
1167 Double_t xdiff = dpad; // reconstructed position constant
1168 Double_t x = dpad; // reconstructed position moved
1169 Float_t ep = pointError; // error of fit
1170 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1171 Float_t signal3 = (Float_t)signals[3]; // signal
1172 Float_t signal2 = (Float_t)signals[2]; // signal
1173 Float_t signal4 = (Float_t)signals[4]; // signal
1174 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1176 //////////////////////////////
1178 /////////////////////////////
1179 if(fDebugLevel > 0){
1180 if ( !fDebugStreamer ) {
1182 TDirectory *backup = gDirectory;
1183 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1184 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1190 Float_t y = ycenter;
1191 (* fDebugStreamer) << "HandlePRFtracklet"<<
1192 "caligroup="<<caligroup<<
1193 "detector="<<detector<<
1196 "npoints="<<npoints<<
1205 "padPosition="<<padPositions[k]<<
1206 "padPosTracklet="<<padPosTracklet<<
1211 "signal1="<<signal1<<
1212 "signal2="<<signal2<<
1213 "signal3="<<signal3<<
1214 "signal4="<<signal4<<
1215 "signal5="<<signal5<<
1221 (* fDebugStreamer) << "HandlePRFtracklet"<<
1222 "caligroup="<<caligroup<<
1223 "detector="<<detector<<
1226 "npoints="<<npoints<<
1235 "padPosition="<<padPositions[k]<<
1236 "padPosTracklet="<<padPosTracklet<<
1241 "signal1="<<signal1<<
1242 "signal2="<<signal2<<
1243 "signal3="<<signal3<<
1244 "signal4="<<signal4<<
1245 "signal5="<<signal5<<
1251 (* fDebugStreamer) << "HandlePRFtracklet"<<
1252 "caligroup="<<caligroup<<
1253 "detector="<<detector<<
1256 "npoints="<<npoints<<
1265 "padPosition="<<padPositions[k]<<
1266 "padPosTracklet="<<padPosTracklet<<
1271 "signal1="<<signal1<<
1272 "signal2="<<signal2<<
1273 "signal3="<<signal3<<
1274 "signal4="<<signal4<<
1275 "signal5="<<signal5<<
1281 ////////////////////////////
1283 ///////////////////////////
1284 if(npoints < fNumberClusters) continue;
1285 if(npoints > fNumberClustersf) continue;
1286 if(nb3pc <= 5) continue;
1287 if((time >= 21) || (time < 7)) continue;
1288 if(TMath::Abs(snp) >= 1.0) continue;
1289 if(TMath::Abs(qcl) < 80) continue;
1291 ////////////////////////////
1293 ///////////////////////////
1295 if(TMath::Abs(dpad) < 1.5) {
1296 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1297 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1299 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1300 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1301 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1303 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1304 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1305 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1309 if(TMath::Abs(dpad) < 1.5) {
1310 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1311 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1313 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1314 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1315 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1317 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1318 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1319 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1323 delete [] padPositions;
1327 //____________Offine tracking in the AliTRDtracker_____________________________
1328 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1331 // PRF width calibration
1332 // Assume a Gaussian shape: determinate the position of the three pad clusters
1333 // Fit with a straight line
1334 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1335 // Fill the PRF as function of angle of the track
1339 //printf("begin\n");
1340 ///////////////////////////////////////////
1341 // Take the parameters of the track
1342 //////////////////////////////////////////
1343 // take now the snp, tnp and tgl from the track
1344 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1345 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1346 if( TMath::Abs(snp) < 1.){
1347 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1349 Double_t tgl = tracklet->GetTgl(); // dz/dl
1350 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1352 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1353 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1354 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1355 // at the end with correction due to linear fit
1356 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1357 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1359 ///////////////////////////////
1360 // Calculate tnp group shift
1361 ///////////////////////////////
1362 Bool_t echec = kFALSE;
1363 Double_t shift = 0.0;
1364 //Calculate the shift in x coresponding to this tnp
1365 if(fNgroupprf != 0.0){
1366 shift = -3.0*(fNgroupprf-1)-1.5;
1367 Double_t limithigh = -0.2*(fNgroupprf-1);
1368 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1370 while(tnp > limithigh){
1376 // do nothing if out of tnp range
1377 //printf("echec %d\n",(Int_t)echec);
1378 if(echec) return kFALSE;
1380 ///////////////////////
1382 //////////////////////
1384 Int_t nb3pc = 0; // number of three pads clusters used for fit
1385 // to see the difference between the fit and the 3 pad clusters position
1386 Double_t padPositions[AliTRDseedV1::kNtb];
1387 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1388 fLinearFitterTracklet->ClearPoints();
1390 //printf("loop clusters \n");
1391 ////////////////////////////
1392 // loop over the clusters
1393 ////////////////////////////
1394 AliTRDcluster *cl = 0x0;
1395 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1396 // reject shared clusters on pad row
1397 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1398 if((cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1400 if(!(cl = tracklet->GetClusters(ic))) continue;
1401 Double_t time = cl->GetPadTime();
1402 if((time<=7) || (time>=21)) continue;
1403 Short_t *signals = cl->GetSignals();
1404 Float_t xcenter = 0.0;
1405 Bool_t echec1 = kTRUE;
1407 /////////////////////////////////////////////////////////////
1408 // Center 3 balanced: position with the center of the pad
1409 /////////////////////////////////////////////////////////////
1410 if ((((Float_t) signals[3]) > 0.0) &&
1411 (((Float_t) signals[2]) > 0.0) &&
1412 (((Float_t) signals[4]) > 0.0)) {
1414 // Security if the denomiateur is 0
1415 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1416 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1417 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1418 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1419 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1425 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1426 if(echec1) continue;
1428 ////////////////////////////////////////////////////////
1429 //if no echec1: calculate with the position of the pad
1430 // Position of the cluster
1431 // fill the linear fitter
1432 ///////////////////////////////////////////////////////
1433 Double_t padPosition = xcenter + cl->GetPadCol();
1434 padPositions[ic] = padPosition;
1436 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1441 //printf("Fin loop clusters \n");
1442 //////////////////////////////
1443 // fit with a straight line
1444 /////////////////////////////
1446 fLinearFitterTracklet->ClearPoints();
1449 fLinearFitterTracklet->Eval();
1451 fLinearFitterTracklet->GetParameters(line);
1452 Float_t pointError = -1.0;
1453 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1454 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1456 fLinearFitterTracklet->ClearPoints();
1458 //printf("PRF second loop \n");
1459 ////////////////////////////////////////////////
1460 // Fill the PRF: Second loop over clusters
1461 //////////////////////////////////////////////
1462 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1463 // reject shared clusters on pad row
1464 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1466 if(!(cl = tracklet->GetClusters(ic))) continue;
1468 Short_t *signals = cl->GetSignals(); // signal
1469 Double_t time = cl->GetPadTime(); // time bin
1470 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1471 Float_t padPos = cl->GetPadCol(); // middle pad
1472 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1473 Float_t ycenter = 0.0; // relative center charge
1474 Float_t ymin = 0.0; // relative left charge
1475 Float_t ymax = 0.0; // relative right charge
1477 ////////////////////////////////////////////////////////////////
1478 // Calculate ycenter, ymin and ymax even for two pad clusters
1479 ////////////////////////////////////////////////////////////////
1480 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1481 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1482 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1483 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1484 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1485 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1488 /////////////////////////
1489 // Calibration group
1490 ////////////////////////
1491 Int_t rowcl = cl->GetPadRow(); // row of cluster
1492 Int_t colcl = cl->GetPadCol(); // col of cluster
1493 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1494 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1495 Float_t xcl = cl->GetY(); // y cluster
1496 Float_t qcl = cl->GetQ(); // charge cluster
1497 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1498 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1499 Double_t xdiff = dpad; // reconstructed position constant
1500 Double_t x = dpad; // reconstructed position moved
1501 Float_t ep = pointError; // error of fit
1502 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1503 Float_t signal3 = (Float_t)signals[3]; // signal
1504 Float_t signal2 = (Float_t)signals[2]; // signal
1505 Float_t signal4 = (Float_t)signals[4]; // signal
1506 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1510 /////////////////////
1512 ////////////////////
1514 if(fDebugLevel > 0){
1515 if ( !fDebugStreamer ) {
1517 TDirectory *backup = gDirectory;
1518 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1519 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1524 Float_t y = ycenter;
1525 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1526 "caligroup="<<caligroup<<
1527 "detector="<<fDetectorPreviousTrack<<
1530 "npoints="<<nbclusters<<
1539 "padPosition="<<padPositions[ic]<<
1540 "padPosTracklet="<<padPosTracklet<<
1545 "signal1="<<signal1<<
1546 "signal2="<<signal2<<
1547 "signal3="<<signal3<<
1548 "signal4="<<signal4<<
1549 "signal5="<<signal5<<
1555 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1556 "caligroup="<<caligroup<<
1557 "detector="<<fDetectorPreviousTrack<<
1560 "npoints="<<nbclusters<<
1569 "padPosition="<<padPositions[ic]<<
1570 "padPosTracklet="<<padPosTracklet<<
1575 "signal1="<<signal1<<
1576 "signal2="<<signal2<<
1577 "signal3="<<signal3<<
1578 "signal4="<<signal4<<
1579 "signal5="<<signal5<<
1585 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1586 "caligroup="<<caligroup<<
1587 "detector="<<fDetectorPreviousTrack<<
1590 "npoints="<<nbclusters<<
1599 "padPosition="<<padPositions[ic]<<
1600 "padPosTracklet="<<padPosTracklet<<
1605 "signal1="<<signal1<<
1606 "signal2="<<signal2<<
1607 "signal3="<<signal3<<
1608 "signal4="<<signal4<<
1609 "signal5="<<signal5<<
1615 /////////////////////
1617 /////////////////////
1618 if(nbclusters < fNumberClusters) continue;
1619 if(nbclusters > fNumberClustersf) continue;
1620 if(nb3pc <= 5) continue;
1621 if((time >= 21) || (time < 7)) continue;
1622 if(TMath::Abs(qcl) < 80) continue;
1623 if( TMath::Abs(snp) > 1.) continue;
1626 ////////////////////////
1628 ///////////////////////
1630 if(TMath::Abs(dpad) < 1.5) {
1631 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1632 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1633 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1635 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1636 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1637 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1639 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1640 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1641 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1646 if(TMath::Abs(dpad) < 1.5) {
1647 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1648 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1650 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1651 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1652 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1654 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1655 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1656 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1659 } // second loop over clusters
1664 ///////////////////////////////////////////////////////////////////////////////////////
1665 // Pad row col stuff: see if masked or not
1666 ///////////////////////////////////////////////////////////////////////////////////////
1667 //_____________________________________________________________________________
1668 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1671 // See if we are not near a masked pad
1674 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1678 //_____________________________________________________________________________
1679 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1682 // See if we are not near a masked pad
1685 if (!IsPadOn(detector, col, row)) {
1686 fGoodTracklet = kFALSE;
1690 if (!IsPadOn(detector, col-1, row)) {
1691 fGoodTracklet = kFALSE;
1696 if (!IsPadOn(detector, col+1, row)) {
1697 fGoodTracklet = kFALSE;
1702 //_____________________________________________________________________________
1703 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1706 // Look in the choosen database if the pad is On.
1707 // If no the track will be "not good"
1710 // Get the parameter object
1711 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1713 AliInfo("Could not get calibDB");
1717 if (!cal->IsChamberInstalled(detector) ||
1718 cal->IsChamberMasked(detector) ||
1719 cal->IsPadMasked(detector,col,row)) {
1727 ///////////////////////////////////////////////////////////////////////////////////////
1728 // Calibration groups: calculate the number of groups, localise...
1729 ////////////////////////////////////////////////////////////////////////////////////////
1730 //_____________________________________________________________________________
1731 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1734 // Calculate the calibration group number for i
1737 // Row of the cluster and position in the pad groups
1739 if (fCalibraMode->GetNnZ(i) != 0) {
1740 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1744 // Col of the cluster and position in the pad groups
1746 if (fCalibraMode->GetNnRphi(i) != 0) {
1747 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1750 return posc*fCalibraMode->GetNfragZ(i)+posr;
1753 //____________________________________________________________________________________
1754 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1757 // Calculate the total number of calibration groups
1763 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1765 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1770 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1772 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1777 fCalibraMode->ModePadCalibration(2,i);
1778 fCalibraMode->ModePadFragmentation(0,2,0,i);
1779 fCalibraMode->SetDetChamb2(i);
1780 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1781 fCalibraMode->ModePadCalibration(0,i);
1782 fCalibraMode->ModePadFragmentation(0,0,0,i);
1783 fCalibraMode->SetDetChamb0(i);
1784 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1785 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1789 //_____________________________________________________________________________
1790 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1793 // Set the mode of calibration group in the z direction for the parameter i
1798 fCalibraMode->SetNz(i, Nz);
1801 AliInfo("You have to choose between 0 and 4");
1805 //_____________________________________________________________________________
1806 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1809 // Set the mode of calibration group in the rphi direction for the parameter i
1814 fCalibraMode->SetNrphi(i ,Nrphi);
1817 AliInfo("You have to choose between 0 and 6");
1822 //_____________________________________________________________________________
1823 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1826 // Set the mode of calibration group all together
1828 if(fVector2d == kTRUE) {
1829 AliInfo("Can not work with the vector method");
1832 fCalibraMode->SetAllTogether(i);
1836 //_____________________________________________________________________________
1837 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1840 // Set the mode of calibration group per supermodule
1842 if(fVector2d == kTRUE) {
1843 AliInfo("Can not work with the vector method");
1846 fCalibraMode->SetPerSuperModule(i);
1850 //____________Set the pad calibration variables for the detector_______________
1851 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1854 // For the detector calcul the first Xbins and set the number of row
1855 // and col pads per calibration groups, the number of calibration
1856 // groups in the detector.
1859 // first Xbins of the detector
1861 fCalibraMode->CalculXBins(detector,0);
1864 fCalibraMode->CalculXBins(detector,1);
1867 fCalibraMode->CalculXBins(detector,2);
1870 // fragmentation of idect
1871 for (Int_t i = 0; i < 3; i++) {
1872 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1873 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1874 , (Int_t) GetStack(detector)
1875 , (Int_t) GetSector(detector),i);
1881 //_____________________________________________________________________________
1882 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1885 // Should be between 0 and 6
1888 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1889 AliInfo("The number of groups must be between 0 and 6!");
1892 fNgroupprf = numberGroupsPRF;
1896 ///////////////////////////////////////////////////////////////////////////////////////////
1897 // Per tracklet: store or reset the info, fill the histos with the info
1898 //////////////////////////////////////////////////////////////////////////////////////////
1899 //_____________________________________________________________________________
1900 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1903 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1904 // Correct from the gain correction before
1907 // time bin of the cluster not corrected
1908 Int_t time = cl->GetPadTime();
1909 Float_t charge = TMath::Abs(cl->GetQ());
1911 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1913 //Correct for the gain coefficient used in the database for reconstruction
1914 Float_t correctthegain = 1.0;
1915 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1916 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1917 Float_t correction = 1.0;
1918 Float_t normalisation = 6.67;
1919 // we divide with gain in AliTRDclusterizer::Transform...
1920 if( correctthegain > 0 ) normalisation /= correctthegain;
1923 // take dd/dl corrected from the angle of the track
1924 correction = dqdl / (normalisation);
1927 // Fill the fAmpTotal with the charge
1929 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1930 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1931 fAmpTotal[(Int_t) group[0]] += correction;
1935 // Fill the fPHPlace and value
1937 fPHPlace[time] = group[1];
1938 fPHValue[time] = charge;
1942 //____________Offine tracking in the AliTRDtracker_____________________________
1943 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1946 // Reset values per tracklet
1949 //Reset good tracklet
1950 fGoodTracklet = kTRUE;
1952 // Reset the fPHValue
1954 //Reset the fPHValue and fPHPlace
1955 for (Int_t k = 0; k < fTimeMax; k++) {
1961 // Reset the fAmpTotal where we put value
1963 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1968 //____________Offine tracking in the AliTRDtracker_____________________________
1969 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1972 // For the offline tracking or mcm tracklets
1973 // This function will be called in the functions UpdateHistogram...
1974 // to fill the info of a track for the relativ gain calibration
1977 Int_t nb = 0; // Nombre de zones traversees
1978 Int_t fd = -1; // Premiere zone non nulle
1979 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1981 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1983 if(nbclusters < fNumberClusters) return;
1984 if(nbclusters > fNumberClustersf) return;
1987 // Normalize with the number of clusters
1988 Double_t normalizeCst = fRelativeScale;
1989 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1991 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1993 // See if the track goes through different zones
1994 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1995 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1996 if (fAmpTotal[k] > 0.0) {
1997 totalcharge += fAmpTotal[k];
2005 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
2011 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2013 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2014 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2017 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2021 if ((fAmpTotal[fd] > 0.0) &&
2022 (fAmpTotal[fd+1] > 0.0)) {
2023 // One of the two very big
2024 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2026 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2027 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2030 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2033 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2035 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2037 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2038 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2041 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2044 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2047 if (fCalibraMode->GetNfragZ(0) > 1) {
2048 if (fAmpTotal[fd] > 0.0) {
2049 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2050 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2051 // One of the two very big
2052 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2054 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2055 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2058 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2061 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2063 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2065 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2066 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2069 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2071 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2082 //____________Offine tracking in the AliTRDtracker_____________________________
2083 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2086 // For the offline tracking or mcm tracklets
2087 // This function will be called in the functions UpdateHistogram...
2088 // to fill the info of a track for the drift velocity calibration
2091 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2092 Int_t fd1 = -1; // Premiere zone non nulle
2093 Int_t fd2 = -1; // Deuxieme zone non nulle
2094 Int_t k1 = -1; // Debut de la premiere zone
2095 Int_t k2 = -1; // Debut de la seconde zone
2096 Int_t nbclusters = 0; // number of clusters
2100 // See if the track goes through different zones
2101 for (Int_t k = 0; k < fTimeMax; k++) {
2102 if (fPHValue[k] > 0.0) {
2108 if (fPHPlace[k] != fd1) {
2114 if (fPHPlace[k] != fd2) {
2121 // See if noise before and after
2122 if(fMaxCluster > 0) {
2123 if(fPHValue[0] > fMaxCluster) return;
2124 if(fTimeMax > fNbMaxCluster) {
2125 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2126 if(fPHValue[k] > fMaxCluster) return;
2131 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2133 if(nbclusters < fNumberClusters) return;
2134 if(nbclusters > fNumberClustersf) return;
2140 for (Int_t i = 0; i < fTimeMax; i++) {
2142 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2144 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2146 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2149 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2151 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2157 if ((fd1 == fd2+1) ||
2159 // One of the two fast all the think
2160 if (k2 > (k1+fDifference)) {
2161 //we choose to fill the fd1 with all the values
2163 for (Int_t i = 0; i < fTimeMax; i++) {
2165 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2167 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2171 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2173 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2178 if ((k2+fDifference) < fTimeMax) {
2179 //we choose to fill the fd2 with all the values
2181 for (Int_t i = 0; i < fTimeMax; i++) {
2183 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2185 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2189 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2191 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2197 // Two zones voisines sinon rien!
2198 if (fCalibraMode->GetNfragZ(1) > 1) {
2200 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2201 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2202 // One of the two fast all the think
2203 if (k2 > (k1+fDifference)) {
2204 //we choose to fill the fd1 with all the values
2206 for (Int_t i = 0; i < fTimeMax; i++) {
2208 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2210 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2214 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2216 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2221 if ((k2+fDifference) < fTimeMax) {
2222 //we choose to fill the fd2 with all the values
2224 for (Int_t i = 0; i < fTimeMax; i++) {
2226 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2228 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2232 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2234 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2241 // Two zones voisines sinon rien!
2243 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2244 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2245 // One of the two fast all the think
2246 if (k2 > (k1 + fDifference)) {
2247 //we choose to fill the fd1 with all the values
2249 for (Int_t i = 0; i < fTimeMax; i++) {
2251 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2253 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2257 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2259 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2264 if ((k2+fDifference) < fTimeMax) {
2265 //we choose to fill the fd2 with all the values
2267 for (Int_t i = 0; i < fTimeMax; i++) {
2269 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2271 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2275 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2277 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2289 //////////////////////////////////////////////////////////////////////////////////////////
2290 // DAQ process functions
2291 /////////////////////////////////////////////////////////////////////////////////////////
2292 //_____________________________________________________________________
2293 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2296 // Event Processing loop - AliTRDrawStreamBase
2297 // TestBeam 2007 version
2298 // 0 timebin problem
2302 // Same algorithm as TestBeam but different reader
2305 Int_t withInput = 1;
2307 Double_t phvalue[16][144][36];
2308 for(Int_t k = 0; k < 36; k++){
2309 for(Int_t j = 0; j < 16; j++){
2310 for(Int_t c = 0; c < 144; c++){
2311 phvalue[j][c][k] = 0.0;
2316 fDetectorPreviousTrack = -1;
2319 Int_t nbtimebin = 0;
2320 Int_t baseline = 10;
2327 while (rawStream->Next()) {
2329 Int_t idetector = rawStream->GetDet(); // current detector
2330 Int_t imcm = rawStream->GetMCM(); // current MCM
2331 Int_t irob = rawStream->GetROB(); // current ROB
2333 //printf("Detector %d\n",idetector);
2335 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2338 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2342 for(Int_t k = 0; k < 36; k++){
2343 for(Int_t j = 0; j < 16; j++){
2344 for(Int_t c = 0; c < 144; c++){
2345 phvalue[j][c][k] = 0.0;
2351 fDetectorPreviousTrack = idetector;
2352 fMCMPrevious = imcm;
2353 fROBPrevious = irob;
2355 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2356 if(nbtimebin == 0) return 0;
2357 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2358 fTimeMax = nbtimebin;
2360 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2361 fNumberClustersf = fTimeMax;
2362 fNumberClusters = (Int_t)(0.6*fTimeMax);
2365 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2366 Int_t col = rawStream->GetCol();
2367 Int_t row = rawStream->GetRow();
2370 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2373 for(Int_t itime = 0; itime < nbtimebin; itime++){
2374 phvalue[row][col][itime] = signal[itime]-baseline;
2378 // fill the last one
2379 if(fDetectorPreviousTrack != -1){
2382 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2385 for(Int_t k = 0; k < 36; k++){
2386 for(Int_t j = 0; j < 16; j++){
2387 for(Int_t c = 0; c < 144; c++){
2388 phvalue[j][c][k] = 0.0;
2397 while (rawStream->Next()) {
2399 Int_t idetector = rawStream->GetDet(); // current detector
2400 Int_t imcm = rawStream->GetMCM(); // current MCM
2401 Int_t irob = rawStream->GetROB(); // current ROB
2403 //printf("Detector %d\n",idetector);
2405 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2408 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2411 for(Int_t k = 0; k < 36; k++){
2412 for(Int_t j = 0; j < 16; j++){
2413 for(Int_t c = 0; c < 144; c++){
2414 phvalue[j][c][k] = 0.0;
2420 fDetectorPreviousTrack = idetector;
2421 fMCMPrevious = imcm;
2422 fROBPrevious = irob;
2424 //baseline = rawStream->GetCommonAdditive(); // common baseline
2426 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2427 fNumberClustersf = fTimeMax;
2428 fNumberClusters = (Int_t)(0.6*fTimeMax);
2429 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2430 Int_t col = rawStream->GetCol();
2431 Int_t row = rawStream->GetRow();
2434 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2436 for(Int_t itime = 0; itime < fTimeMax; itime++){
2437 phvalue[row][col][itime] = signal[itime]-baseline;
2441 // fill the last one
2442 if(fDetectorPreviousTrack != -1){
2445 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2448 for(Int_t k = 0; k < 36; k++){
2449 for(Int_t j = 0; j < 16; j++){
2450 for(Int_t c = 0; c < 144; c++){
2451 phvalue[j][c][k] = 0.0;
2461 //_____________________________________________________________________
2462 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2465 // Event processing loop - AliRawReader
2466 // Testbeam 2007 version
2469 AliTRDrawStreamBase rawStream(rawReader);
2471 rawReader->Select("TRD");
2473 return ProcessEventDAQ(&rawStream, nocheck);
2476 //_________________________________________________________________________
2477 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2479 eventHeaderStruct *event,
2482 eventHeaderStruct* /*event*/,
2489 // process date event
2490 // Testbeam 2007 version
2493 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2494 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2498 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2503 //////////////////////////////////////////////////////////////////////////////
2504 // Routine inside the DAQ process
2505 /////////////////////////////////////////////////////////////////////////////
2506 //_______________________________________________________________________
2507 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2510 // Look for the maximum by collapsing over the time
2511 // Sum over four pad col and two pad row
2517 Int_t idect = fDetectorPreviousTrack;
2518 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2520 for(Int_t tb = 0; tb < 36; tb++){
2524 //fGoodTracklet = kTRUE;
2525 //fDetectorPreviousTrack = 0;
2528 ///////////////////////////
2530 /////////////////////////
2534 Double_t integralMax = -1;
2536 for (Int_t ir = 1; ir <= 15; ir++)
2538 for (Int_t ic = 2; ic <= 142; ic++)
2540 Double_t integral = 0;
2541 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2543 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2545 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2546 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2549 for(Int_t tb = 0; tb< fTimeMax; tb++){
2550 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2556 if (integralMax < integral)
2560 integralMax = integral;
2561 } // check max integral
2565 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2567 if((imaxRow == 0) || (imaxCol == 0)) {
2571 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2572 //if(!fGoodTracklet) used = 1;;
2574 // /////////////////////////////////////////////////////
2575 // sum ober 2 row and 4 pad cols for each time bins
2576 // ////////////////////////////////////////////////////
2579 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2581 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2583 for(Int_t it = 0; it < fTimeMax; it++){
2584 sum[it] += phvalue[ir][ic][it];
2590 Double_t sumcharge = 0.0;
2591 for(Int_t it = 0; it < fTimeMax; it++){
2592 sumcharge += sum[it];
2593 if(sum[it] > 20.0) nbcl++;
2597 /////////////////////////////////////////////////////////
2599 ////////////////////////////////////////////////////////
2600 if(fDebugLevel > 0){
2601 if ( !fDebugStreamer ) {
2603 TDirectory *backup = gDirectory;
2604 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2605 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2608 Double_t amph0 = sum[0];
2609 Double_t amphlast = sum[fTimeMax-1];
2610 Double_t rms = TMath::RMS(fTimeMax,sum);
2611 Int_t goodtracklet = (Int_t) fGoodTracklet;
2612 for(Int_t it = 0; it < fTimeMax; it++){
2613 Double_t clustera = sum[it];
2615 (* fDebugStreamer) << "FillDAQa"<<
2616 "ampTotal="<<sumcharge<<
2619 "detector="<<idect<<
2621 "amphlast="<<amphlast<<
2622 "goodtracklet="<<goodtracklet<<
2623 "clustera="<<clustera<<
2630 ////////////////////////////////////////////////////////
2632 ///////////////////////////////////////////////////////
2633 if(sum[0] > 100.0) used = 1;
2634 if(nbcl < fNumberClusters) used = 1;
2635 if(nbcl > fNumberClustersf) used = 1;
2637 //if(fDetectorPreviousTrack == 15){
2638 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2640 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2642 for(Int_t it = 0; it < fTimeMax; it++){
2643 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2645 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2650 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2652 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2659 //____________Online trackling in AliTRDtrigger________________________________
2660 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2664 // Fill a simple average pulse height
2668 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2674 //____________Write_____________________________________________________
2675 //_____________________________________________________________________
2676 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2679 // Write infos to file
2683 if ( fDebugStreamer ) {
2684 delete fDebugStreamer;
2685 fDebugStreamer = 0x0;
2688 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2693 ,fNumberUsedPh[1]));
2695 TDirectory *backup = gDirectory;
2701 option = "recreate";
2703 TFile f(filename,option.Data());
2705 TStopwatch stopwatch;
2708 f.WriteTObject(fCalibraVector);
2713 f.WriteTObject(fCH2d);
2718 f.WriteTObject(fPH2d);
2723 f.WriteTObject(fPRF2d);
2726 if(fLinearFitterOn){
2727 AnalyseLinearFitter();
2728 f.WriteTObject(fLinearVdriftFit);
2733 if ( backup ) backup->cd();
2735 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2736 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2738 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2740 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2741 //___________________________________________probe the histos__________________________________________________
2742 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2745 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2746 // debug mode with 2 for TH2I and 3 for TProfile2D
2747 // It gives a pointer to a Double_t[7] with the info following...
2748 // [0] : number of calibration groups with entries
2749 // [1] : minimal number of entries found
2750 // [2] : calibration group number of the min
2751 // [3] : maximal number of entries found
2752 // [4] : calibration group number of the max
2753 // [5] : mean number of entries found
2754 // [6] : mean relative error
2757 Double_t *info = new Double_t[7];
2759 // Number of Xbins (detectors or groups of pads)
2760 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2761 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2764 Double_t nbwe = 0; //number of calibration groups with entries
2765 Double_t minentries = 0; //minimal number of entries found
2766 Double_t maxentries = 0; //maximal number of entries found
2767 Double_t placemin = 0; //calibration group number of the min
2768 Double_t placemax = -1; //calibration group number of the max
2769 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2770 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2772 Double_t counter = 0;
2775 TH1F *nbEntries = 0x0;//distribution of the number of entries
2776 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2777 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2779 // Beginning of the loop over the calibration groups
2780 for (Int_t idect = 0; idect < nbins; idect++) {
2782 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2783 projch->SetDirectory(0);
2785 // Number of entries for this calibration group
2786 Double_t nentries = 0.0;
2788 for (Int_t k = 0; k < nxbins; k++) {
2789 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2793 for (Int_t k = 0; k < nxbins; k++) {
2794 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2795 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2796 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2804 if((!((Bool_t)nbEntries)) && (nentries > 0)){
2805 nbEntries = new TH1F("Number of entries","Number of entries"
2806 ,100,(Int_t)nentries/2,nentries*2);
2807 nbEntries->SetDirectory(0);
2808 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
2810 nbEntriesPerGroup->SetDirectory(0);
2811 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
2812 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2813 nbEntriesPerSp->SetDirectory(0);
2816 if(nentries > 0) nbEntries->Fill(nentries);
2817 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2818 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2823 if(nentries > maxentries){
2824 maxentries = nentries;
2828 minentries = nentries;
2830 if(nentries < minentries){
2831 minentries = nentries;
2837 meanstats += nentries;
2839 }//calibration groups loop
2841 if(nbwe > 0) meanstats /= nbwe;
2842 if(counter > 0) meanrelativerror /= counter;
2844 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2845 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2846 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2847 AliInfo(Form("The mean number of entries is %f",meanstats));
2848 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2851 info[1] = minentries;
2853 info[3] = maxentries;
2855 info[5] = meanstats;
2856 info[6] = meanrelativerror;
2859 gStyle->SetPalette(1);
2860 gStyle->SetOptStat(1111);
2861 gStyle->SetPadBorderMode(0);
2862 gStyle->SetCanvasColor(10);
2863 gStyle->SetPadLeftMargin(0.13);
2864 gStyle->SetPadRightMargin(0.01);
2865 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2868 nbEntries->Draw("");
2870 nbEntriesPerSp->SetStats(0);
2871 nbEntriesPerSp->Draw("");
2872 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2874 nbEntriesPerGroup->SetStats(0);
2875 nbEntriesPerGroup->Draw("");
2881 //____________________________________________________________________________
2882 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2885 // Return a Int_t[4] with:
2886 // 0 Mean number of entries
2887 // 1 median of number of entries
2888 // 2 rms of number of entries
2889 // 3 number of group with entries
2892 Double_t *stat = new Double_t[4];
2895 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2896 Double_t *weight = new Double_t[nbofgroups];
2897 Int_t *nonul = new Int_t[nbofgroups];
2899 for(Int_t k = 0; k < nbofgroups; k++){
2900 if(fEntriesCH[k] > 0) {
2902 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2905 else weight[k] = 0.0;
2907 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2908 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2909 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2914 //____________________________________________________________________________
2915 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2918 // Return a Int_t[4] with:
2919 // 0 Mean number of entries
2920 // 1 median of number of entries
2921 // 2 rms of number of entries
2922 // 3 number of group with entries
2925 Double_t *stat = new Double_t[4];
2928 Int_t nbofgroups = 540;
2929 Double_t *weight = new Double_t[nbofgroups];
2930 Int_t *nonul = new Int_t[nbofgroups];
2932 for(Int_t k = 0; k < nbofgroups; k++){
2933 if(fEntriesLinearFitter[k] > 0) {
2935 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2938 else weight[k] = 0.0;
2940 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2941 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2942 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2947 //////////////////////////////////////////////////////////////////////////////////////
2949 //////////////////////////////////////////////////////////////////////////////////////
2950 //_____________________________________________________________________________
2951 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2954 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2955 // If fNgroupprf is zero then no binning in tnp
2959 name += fCalibraMode->GetNz(2);
2961 name += fCalibraMode->GetNrphi(2);
2965 if(fNgroupprf != 0){
2967 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2968 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2969 fPRF2d->SetYTitle("Det/pad groups");
2970 fPRF2d->SetXTitle("Position x/W [pad width units]");
2971 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2972 fPRF2d->SetStats(0);
2975 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2976 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2977 fPRF2d->SetYTitle("Det/pad groups");
2978 fPRF2d->SetXTitle("Position x/W [pad width units]");
2979 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2980 fPRF2d->SetStats(0);
2985 //_____________________________________________________________________________
2986 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2989 // Create the 2D histos
2993 name += fCalibraMode->GetNz(1);
2995 name += fCalibraMode->GetNrphi(1);
2997 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2998 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3000 fPH2d->SetYTitle("Det/pad groups");
3001 fPH2d->SetXTitle("time [#mus]");
3002 fPH2d->SetZTitle("<PH> [a.u.]");
3006 //_____________________________________________________________________________
3007 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3010 // Create the 2D histos
3014 name += fCalibraMode->GetNz(0);
3016 name += fCalibraMode->GetNrphi(0);
3018 fCH2d = new TH2I("CH2d",(const Char_t *) name
3019 ,fNumberBinCharge,0,300,nn,0,nn);
3020 fCH2d->SetYTitle("Det/pad groups");
3021 fCH2d->SetXTitle("charge deposit [a.u]");
3022 fCH2d->SetZTitle("counts");
3027 //////////////////////////////////////////////////////////////////////////////////
3028 // Set relative scale
3029 /////////////////////////////////////////////////////////////////////////////////
3030 //_____________________________________________________________________________
3031 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3034 // Set the factor that will divide the deposited charge
3035 // to fit in the histo range [0,300]
3038 if (RelativeScale > 0.0) {
3039 fRelativeScale = RelativeScale;
3042 AliInfo("RelativeScale must be strict positif!");
3046 //////////////////////////////////////////////////////////////////////////////////
3047 // Quick way to fill a histo
3048 //////////////////////////////////////////////////////////////////////////////////
3049 //_____________________________________________________________________
3050 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3053 // FillCH2d: Marian style
3056 //skip simply the value out of range
3057 if((y>=300.0) || (y<0.0)) return;
3059 //Calcul the y place
3060 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3061 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3064 fCH2d->GetArray()[place]++;
3068 //////////////////////////////////////////////////////////////////////////////////
3069 // Geometrical functions
3070 ///////////////////////////////////////////////////////////////////////////////////
3071 //_____________________________________________________________________________
3072 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3075 // Reconstruct the layer number from the detector number
3078 return ((Int_t) (d % 6));
3082 //_____________________________________________________________________________
3083 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3086 // Reconstruct the stack number from the detector number
3088 const Int_t kNlayer = 6;
3090 return ((Int_t) (d % 30) / kNlayer);
3094 //_____________________________________________________________________________
3095 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3098 // Reconstruct the sector number from the detector number
3102 return ((Int_t) (d / fg));
3105 ///////////////////////////////////////////////////////////////////////////////////
3106 // Getter functions for DAQ of the CH2d and the PH2d
3107 //////////////////////////////////////////////////////////////////////////////////
3108 //_____________________________________________________________________
3109 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3112 // return pointer to fPH2d TProfile2D
3113 // create a new TProfile2D if it doesn't exist allready
3119 fTimeMax = nbtimebin;
3120 fSf = samplefrequency;
3126 //_____________________________________________________________________
3127 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3130 // return pointer to fCH2d TH2I
3131 // create a new TH2I if it doesn't exist allready
3140 ////////////////////////////////////////////////////////////////////////////////////////////
3141 // Drift velocity calibration
3142 ///////////////////////////////////////////////////////////////////////////////////////////
3143 //_____________________________________________________________________
3144 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3147 // return pointer to TLinearFitter Calibration
3148 // if force is true create a new TLinearFitter if it doesn't exist allready
3151 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3152 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3155 // if we are forced and TLinearFitter doesn't yet exist create it
3157 // new TLinearFitter
3158 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3159 fLinearFitterArray.AddAt(linearfitter,detector);
3160 return linearfitter;
3163 //____________________________________________________________________________
3164 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3167 // Analyse array of linear fitter because can not be written
3168 // Store two arrays: one with the param the other one with the error param + number of entries
3171 for(Int_t k = 0; k < 540; k++){
3172 TLinearFitter *linearfitter = GetLinearFitter(k);
3173 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3174 TVectorD *par = new TVectorD(2);
3175 TVectorD pare = TVectorD(2);
3176 TVectorD *parE = new TVectorD(3);
3177 linearfitter->Eval();
3178 linearfitter->GetParameters(*par);
3179 linearfitter->GetErrors(pare);
3180 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3181 (*parE)[0] = pare[0]*ppointError;
3182 (*parE)[1] = pare[1]*ppointError;
3183 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3184 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3185 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);