1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
20 // AliTRDCalibraFillHisto
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fills histos or vectors.
24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
25 // The user has to choose with the functions SetNz and SetNrphi the precision of the calibration (see AliTRDCalibraMode).
26 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
27 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
28 // in the function "FollowBackProlongation" (AliTRDtracker)
29 // Per default the functions to fill are off.
32 // R. Bailhache (R.Bailhache@gsi.de)
34 //////////////////////////////////////////////////////////////////////////////////////
36 #include <TProfile2D.h>
41 #include <TObjArray.h>
46 #include <TStopwatch.h>
48 #include <TDirectory.h>
49 #include <TTreeStream.h>
54 #include "AliTRDCalibraFillHisto.h"
55 #include "AliTRDCalibraMode.h"
56 #include "AliTRDCalibraVector.h"
57 #include "AliTRDCalibraVdriftLinearFit.h"
58 #include "AliTRDcalibDB.h"
59 #include "AliTRDCommonParam.h"
60 #include "AliTRDmcmTracklet.h"
61 #include "AliTRDmcm.h"
62 #include "AliTRDtrigParam.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcluster.h"
65 #include "AliTRDtrackV1.h"
66 #include "AliTRDrawStreamTB.h"
67 #include "AliRawReader.h"
68 #include "AliRawReaderDate.h"
69 #include "AliTRDgeometry.h"
70 #include "./Cal/AliTRDCalROC.h"
71 #include "./Cal/AliTRDCalDet.h"
78 ClassImp(AliTRDCalibraFillHisto)
80 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
81 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
83 //_____________singleton implementation_________________________________________________
84 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
87 // Singleton implementation
90 if (fgTerminated != kFALSE) {
94 if (fgInstance == 0) {
95 fgInstance = new AliTRDCalibraFillHisto();
102 //______________________________________________________________________________________
103 void AliTRDCalibraFillHisto::Terminate()
106 // Singleton implementation
107 // Deletes the instance of this class
110 fgTerminated = kTRUE;
112 if (fgInstance != 0) {
119 //______________________________________________________________________________________
120 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
123 ,fMcmCorrectAngle(kFALSE)
129 ,fLinearFitterOn(kFALSE)
130 ,fLinearFitterDebugOn(kFALSE)
132 ,fThresholdClusterPRF2(15.0)
133 ,fCalibraMode(new AliTRDCalibraMode())
136 ,fDetectorPreviousTrack(-1)
140 ,fNumberClustersf(30)
146 ,fNumberBinCharge(100)
152 ,fGoodTracklet(kTRUE)
154 ,fEntriesLinearFitter(0x0)
159 ,fLinearFitterArray(540)
160 ,fLinearVdriftFit(0x0)
165 // Default constructor
169 // Init some default values
172 fNumberUsedCh[0] = 0;
173 fNumberUsedCh[1] = 0;
174 fNumberUsedPh[0] = 0;
175 fNumberUsedPh[1] = 0;
177 fGeo = new AliTRDgeometry();
181 //______________________________________________________________________________________
182 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
185 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
188 ,fPRF2dOn(c.fPRF2dOn)
189 ,fHisto2d(c.fHisto2d)
190 ,fVector2d(c.fVector2d)
191 ,fLinearFitterOn(c.fLinearFitterOn)
192 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
193 ,fRelativeScale(c.fRelativeScale)
194 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
197 ,fDebugLevel(c.fDebugLevel)
198 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
199 ,fMCMPrevious(c.fMCMPrevious)
200 ,fROBPrevious(c.fROBPrevious)
201 ,fNumberClusters(c.fNumberClusters)
202 ,fNumberClustersf(c.fNumberClustersf)
203 ,fProcent(c.fProcent)
204 ,fDifference(c.fDifference)
205 ,fNumberTrack(c.fNumberTrack)
206 ,fTimeMax(c.fTimeMax)
208 ,fNumberBinCharge(c.fNumberBinCharge)
209 ,fNumberBinPRF(c.fNumberBinPRF)
210 ,fNgroupprf(c.fNgroupprf)
214 ,fGoodTracklet(c.fGoodTracklet)
216 ,fEntriesLinearFitter(0x0)
221 ,fLinearFitterArray(540)
222 ,fLinearVdriftFit(0x0)
229 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
230 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
232 fPH2d = new TProfile2D(*c.fPH2d);
233 fPH2d->SetDirectory(0);
236 fPRF2d = new TProfile2D(*c.fPRF2d);
237 fPRF2d->SetDirectory(0);
240 fCH2d = new TH2I(*c.fCH2d);
241 fCH2d->SetDirectory(0);
243 if(c.fLinearVdriftFit){
244 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
247 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
248 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
253 fGeo = new AliTRDgeometry();
256 //____________________________________________________________________________________
257 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
260 // AliTRDCalibraFillHisto destructor
264 if ( fDebugStreamer ) delete fDebugStreamer;
266 if ( fCalDetGain ) delete fCalDetGain;
267 if ( fCalROCGain ) delete fCalROCGain;
274 //_____________________________________________________________________________
275 void AliTRDCalibraFillHisto::Destroy()
286 //_____________________________________________________________________________
287 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
290 // Delete DebugStreamer
293 if ( fDebugStreamer ) delete fDebugStreamer;
296 //_____________________________________________________________________________
297 void AliTRDCalibraFillHisto::ClearHistos()
317 //////////////////////////////////////////////////////////////////////////////////
318 // calibration with AliTRDtrackV1: Init, Update
319 //////////////////////////////////////////////////////////////////////////////////
320 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
321 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
324 // Init the histograms and stuff to be filled
329 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
331 AliInfo("Could not get calibDB");
334 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
336 AliInfo("Could not get CommonParam");
341 fTimeMax = cal->GetNumberOfTimeBins();
342 fSf = parCom->GetSamplingFrequency();
345 //calib object from database used for reconstruction
346 if(fCalDetGain) delete fCalDetGain;
347 fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
349 // Calcul Xbins Chambd0, Chamb2
350 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
351 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
352 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
354 // If vector method On initialised all the stuff
356 fCalibraVector = new AliTRDCalibraVector();
357 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
358 fCalibraVector->SetTimeMax(fTimeMax);
359 if(fNgroupprf != 0) {
360 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
361 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
364 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
365 fCalibraVector->SetPRFRange(1.5);
367 for(Int_t k = 0; k < 3; k++){
368 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
369 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
371 TString namech("Nz");
372 namech += fCalibraMode->GetNz(0);
374 namech += fCalibraMode->GetNrphi(0);
375 fCalibraVector->SetNameCH((const char* ) namech);
376 TString nameph("Nz");
377 nameph += fCalibraMode->GetNz(1);
379 nameph += fCalibraMode->GetNrphi(1);
380 fCalibraVector->SetNamePH((const char* ) nameph);
381 TString nameprf("Nz");
382 nameprf += fCalibraMode->GetNz(2);
384 nameprf += fCalibraMode->GetNrphi(2);
386 nameprf += fNgroupprf;
387 fCalibraVector->SetNamePRF((const char* ) nameprf);
390 // Create the 2D histos corresponding to the pad groupCalibration mode
393 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
394 ,fCalibraMode->GetNz(0)
395 ,fCalibraMode->GetNrphi(0)));
397 // Create the 2D histo
402 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
403 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
407 fEntriesCH = new Int_t[ntotal0];
408 for(Int_t k = 0; k < ntotal0; k++){
415 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
416 ,fCalibraMode->GetNz(1)
417 ,fCalibraMode->GetNrphi(1)));
419 // Create the 2D histo
424 fPHPlace = new Short_t[fTimeMax];
425 for (Int_t k = 0; k < fTimeMax; k++) {
428 fPHValue = new Float_t[fTimeMax];
429 for (Int_t k = 0; k < fTimeMax; k++) {
433 if (fLinearFitterOn) {
434 //fLinearFitterArray.Expand(540);
435 fLinearFitterArray.SetName("ArrayLinearFitters");
436 fEntriesLinearFitter = new Int_t[540];
437 for(Int_t k = 0; k < 540; k++){
438 fEntriesLinearFitter[k] = 0;
440 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
445 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
446 ,fCalibraMode->GetNz(2)
447 ,fCalibraMode->GetNrphi(2)));
448 // Create the 2D histo
450 CreatePRF2d(ntotal2);
457 //____________Offline tracking in the AliTRDtracker____________________________
458 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
461 // Use AliTRDtrack for the calibration
465 AliTRDcluster *cl = 0x0; // pointeur to cluster
466 Int_t index0 = 0; // index of the first cluster in the current chamber
467 Int_t index1 = 0; // index of the current cluster in the current chamber
469 //////////////////////////////
470 // loop over the clusters
471 ///////////////////////////////
472 while((cl = t->GetCluster(index1))){
474 /////////////////////////////////////////////////////////////////////////
475 // Fill the infos for the previous clusters if not the same detector
476 ////////////////////////////////////////////////////////////////////////
477 Int_t detector = cl->GetDetector();
478 if ((detector != fDetectorPreviousTrack) &&
479 (index0 != index1)) {
483 //If the same track, then look if the previous detector is in
484 //the same plane, if yes: not a good track
485 if ((GetPlane(detector) == GetPlane(fDetectorPreviousTrack))) {
489 // Fill only if the track doesn't touch a masked pad or doesn't
492 // drift velocity unables to cut bad tracklets
493 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
497 FillTheInfoOfTheTrackCH(index1-index0);
502 FillTheInfoOfTheTrackPH();
505 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
508 } // if a good tracklet
511 ResetfVariablestracklet();
514 } // Fill at the end the charge
517 /////////////////////////////////////////////////////////////
518 // Localise and take the calib gain object for the detector
519 ////////////////////////////////////////////////////////////
520 if (detector != fDetectorPreviousTrack) {
522 //Localise the detector
523 LocalisationDetectorXbins(detector);
526 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
528 AliInfo("Could not get calibDB");
533 if( fCalROCGain ) delete fCalROCGain;
534 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
538 // Reset the detectbjobsor
539 fDetectorPreviousTrack = detector;
541 ///////////////////////////////////////
542 // Store the info of the cluster
543 ///////////////////////////////////////
544 Int_t row = cl->GetPadRow();
545 Int_t col = cl->GetPadCol();
546 CheckGoodTracklet(detector,row,col);
547 Int_t group[2] = {0,0};
548 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
549 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
550 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
554 } // while on clusters
556 ///////////////////////////
557 // Fill the last plane
558 //////////////////////////
559 if( index0 != index1 ){
565 // drift velocity unables to cut bad tracklets
566 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
570 FillTheInfoOfTheTrackCH(index1-index0);
575 FillTheInfoOfTheTrackPH();
578 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
580 } // if a good tracklet
585 ResetfVariablestracklet();
590 //____________Offline tracking in the AliTRDtracker____________________________
591 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
594 // Use AliTRDtrackV1 for the calibration
598 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
599 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
602 ///////////////////////////
603 // loop over the tracklet
604 ///////////////////////////
605 for(Int_t itr = 0; itr < 6; itr++){
607 if(!(tracklet = t->GetTracklet(itr))) continue;
608 if(!tracklet->IsOK()) continue;
610 ResetfVariablestracklet();
612 //////////////////////////////////////////
613 // localisation of the tracklet and dqdl
614 //////////////////////////////////////////
615 Int_t plane = tracklet->GetPlane();
617 while(!(cl = tracklet->GetClusters(ic++))) continue;
618 Int_t detector = cl->GetDetector();
619 if (detector != fDetectorPreviousTrack) {
620 // don't use the rest of this track if in the same plane
621 if ((plane == GetPlane(fDetectorPreviousTrack))) {
624 //Localise the detector bin
625 LocalisationDetectorXbins(detector);
627 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
629 AliInfo("Could not get calibDB");
633 if( fCalROCGain ) delete fCalROCGain;
634 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
636 fDetectorPreviousTrack = detector;
640 ////////////////////////////
641 // loop over the clusters
642 ////////////////////////////
643 Int_t nbclusters = 0;
644 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
645 if(!(cl = tracklet->GetClusters(ic))) continue;
648 // Store the info bis of the tracklet
649 Int_t row = cl->GetPadRow();
650 Int_t col = cl->GetPadCol();
651 CheckGoodTracklet(detector,row,col);
652 Int_t group[2] = {0,0};
653 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
654 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
655 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(ic),group,row,col);
658 ////////////////////////////////////////
659 // Fill the stuffs if a good tracklet
660 ////////////////////////////////////////
663 // drift velocity unables to cut bad tracklets
664 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
668 FillTheInfoOfTheTrackCH(nbclusters);
673 FillTheInfoOfTheTrackPH();
676 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
678 } // if a good tracklet
686 ///////////////////////////////////////////////////////////////////////////////////
687 // Routine inside the update with AliTRDtrack
688 ///////////////////////////////////////////////////////////////////////////////////
689 //____________Offine tracking in the AliTRDtracker_____________________________
690 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
693 // Drift velocity calibration:
694 // Fit the clusters with a straight line
695 // From the slope find the drift velocity
698 //Number of points: if less than 3 return kFALSE
699 Int_t npoints = index1-index0;
700 if(npoints <= 2) return kFALSE;
705 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
706 // parameters of the track
707 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in the plan yx track
708 Double_t tgl = t->GetTglPlane(GetPlane(detector)); // dz/dl and not dz/dx!
709 Double_t tnp = 0.0; // tan angle in the plan xy track
710 if( TMath::Abs(snp) < 1.){
711 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
713 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
714 // tilting pad and cross row
715 Int_t crossrow = 0; // if it crosses a pad row
716 Int_t rowp = -1; // if it crosses a pad row
717 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
718 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
719 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
721 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
722 Double_t dydt = 0.0; // dydt tracklet after straight line fit
723 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
724 Double_t pointError = 0.0; // error after straight line fit
725 Int_t nbli = 0; // number linear fitter points
726 linearFitterTracklet.StoreData(kFALSE);
727 linearFitterTracklet.ClearPoints();
729 //////////////////////////////
730 // loop over clusters
731 ////////////////////////////
732 for(Int_t k = 0; k < npoints; k++){
734 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
735 if(!cl->IsInChamber()) continue;
736 Double_t ycluster = cl->GetY();
737 Int_t time = cl->GetPadTime();
738 Double_t timeis = time/fSf;
739 //Double_t q = cl->GetQ();
740 //See if cross two pad rows
741 Int_t row = cl->GetPadRow();
743 if(row != rowp) crossrow = 1;
745 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
750 //////////////////////////////
752 //////////////////////////////
753 if(nbli <= 2) return kFALSE;
755 linearFitterTracklet.Eval();
756 linearFitterTracklet.GetParameters(pars);
757 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
758 errorpar = linearFitterTracklet.GetParError(1)*pointError;
761 /////////////////////////////
763 ////////////////////////////
765 if ( !fDebugStreamer ) {
767 TDirectory *backup = gDirectory;
768 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
769 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
773 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
774 //"snpright="<<snpright<<
775 "npoints="<<npoints<<
777 "detector="<<detector<<
784 "crossrow="<<crossrow<<
785 "errorpar="<<errorpar<<
786 "pointError="<<pointError<<
790 Int_t nbclusters = index1-index0;
791 Int_t plane = GetPlane(fDetectorPreviousTrack);
793 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
794 //"snpright="<<snpright<<
795 "nbclusters="<<nbclusters<<
796 "detector="<<fDetectorPreviousTrack<<
802 ///////////////////////////
804 ///////////////////////////
805 if(npoints < fNumberClusters) return kFALSE;
806 if(npoints > fNumberClustersf) return kFALSE;
807 if(pointError >= 0.1) return kFALSE;
808 if(crossrow == 1) return kFALSE;
810 ////////////////////////////
812 ////////////////////////////
814 //Add to the linear fitter of the detector
815 if( TMath::Abs(snp) < 1.){
816 Double_t x = tnp-dzdx*tnt;
817 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
818 if(fLinearFitterDebugOn) {
819 fLinearVdriftFit->Update(detector,x,pars[1]);
821 fEntriesLinearFitter[detector]++;
827 //____________Offine tracking in the AliTRDtracker_____________________________
828 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
831 // Drift velocity calibration:
832 // Fit the clusters with a straight line
833 // From the slope find the drift velocity
836 ////////////////////////////////////////////////
837 //Number of points: if less than 3 return kFALSE
838 /////////////////////////////////////////////////
839 if(nbclusters <= 2) return kFALSE;
844 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
845 // results of the linear fit
846 Double_t dydt = 0.0; // dydt tracklet after straight line fit
847 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
848 Double_t pointError = 0.0; // error after straight line fit
849 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
850 Int_t crossrow = 0; // if it crosses a pad row
851 Int_t rowp = -1; // if it crosses a pad row
852 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
854 linearFitterTracklet.StoreData(kFALSE);
855 linearFitterTracklet.ClearPoints();
857 ///////////////////////////////////////////
858 // Take the parameters of the track
859 //////////////////////////////////////////
860 // take now the snp, tnp and tgl from the track
861 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
862 Double_t tnp = 0.0; // dy/dx at the end of the chamber
863 if( TMath::Abs(snp) < 1.){
864 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
866 Double_t tgl = tracklet->GetTgl(); // dz/dl
867 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
869 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
870 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
871 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
872 // at the end with correction due to linear fit
873 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
874 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
877 ////////////////////////////
878 // loop over the clusters
879 ////////////////////////////
881 AliTRDcluster *cl = 0x0;
882 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
883 if(!(cl = tracklet->GetClusters(ic))) continue;
884 if(!cl->IsInChamber()) continue;
886 Double_t ycluster = cl->GetY();
887 Int_t time = cl->GetPadTime();
888 Double_t timeis = time/fSf;
889 //See if cross two pad rows
890 Int_t row = cl->GetPadRow();
891 if(rowp==-1) rowp = row;
892 if(row != rowp) crossrow = 1;
894 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
899 ////////////////////////////////////
900 // Do the straight line fit now
901 ///////////////////////////////////
903 linearFitterTracklet.Eval();
904 linearFitterTracklet.GetParameters(pars);
905 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
906 errorpar = linearFitterTracklet.GetParError(1)*pointError;
909 ////////////////////////////////
911 ///////////////////////////////
915 if ( !fDebugStreamer ) {
917 TDirectory *backup = gDirectory;
918 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
919 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
923 Int_t plane = GetPlane(fDetectorPreviousTrack);
925 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
926 //"snpright="<<snpright<<
927 "nbclusters="<<nbclusters<<
928 "detector="<<fDetectorPreviousTrack<<
936 "crossrow="<<crossrow<<
937 "errorpar="<<errorpar<<
938 "pointError="<<pointError<<
943 /////////////////////////
945 ////////////////////////
947 if(nbclusters < fNumberClusters) return kFALSE;
948 if(nbclusters > fNumberClustersf) return kFALSE;
949 if(pointError >= 0.1) return kFALSE;
950 if(crossrow == 1) return kFALSE;
952 ///////////////////////
954 //////////////////////
957 //Add to the linear fitter of the detector
958 if( TMath::Abs(snp) < 1.){
959 Double_t x = tnp-dzdx*tnt;
960 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
961 if(fLinearFitterDebugOn) {
962 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
964 fEntriesLinearFitter[fDetectorPreviousTrack]++;
970 //____________Offine tracking in the AliTRDtracker_____________________________
971 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
974 // PRF width calibration
975 // Assume a Gaussian shape: determinate the position of the three pad clusters
976 // Fit with a straight line
977 // Take the fitted values for all the clusters (3 or 2 pad clusters)
978 // Fill the PRF as function of angle of the track
983 //////////////////////////
985 /////////////////////////
986 Int_t npoints = index1-index0; // number of total points
987 Int_t nb3pc = 0; // number of three pads clusters used for fit
988 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
989 // To see the difference due to the fit
990 Double_t *padPositions;
991 padPositions = new Double_t[npoints];
992 for(Int_t k = 0; k < npoints; k++){
993 padPositions[k] = 0.0;
995 // Take the tgl and snp with the track t now
996 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
997 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
998 Float_t dzdx = 0.0; // dzdx
1000 if(TMath::Abs(snp) < 1.0){
1001 tnp = snp / (TMath::Sqrt(1-snp*snp));
1002 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1005 TLinearFitter fitter(2,"pol1");
1006 fitter.StoreData(kFALSE);
1007 fitter.ClearPoints();
1010 ///////////////////////////
1011 // calcul the tnp group
1012 ///////////////////////////
1013 Bool_t echec = kFALSE;
1014 Double_t shift = 0.0;
1015 //Calculate the shift in x coresponding to this tnp
1016 if(fNgroupprf != 0.0){
1017 shift = -3.0*(fNgroupprf-1)-1.5;
1018 Double_t limithigh = -0.2*(fNgroupprf-1);
1019 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1021 while(tnp > limithigh){
1027 if(echec) return kFALSE;
1030 //////////////////////
1032 /////////////////////
1033 for(Int_t k = 0; k < npoints; k++){
1035 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1036 Short_t *signals = cl->GetSignals();
1037 Double_t time = cl->GetPadTime();
1038 //Calculate x if possible
1039 Float_t xcenter = 0.0;
1040 Bool_t echec = kTRUE;
1041 if((time<=7) || (time>=21)) continue;
1042 // Center 3 balanced: position with the center of the pad
1043 if ((((Float_t) signals[3]) > 0.0) &&
1044 (((Float_t) signals[2]) > 0.0) &&
1045 (((Float_t) signals[4]) > 0.0)) {
1047 // Security if the denomiateur is 0
1048 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1049 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1050 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1051 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1052 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1058 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1060 //if no echec: calculate with the position of the pad
1061 // Position of the cluster
1062 Double_t padPosition = xcenter + cl->GetPadCol();
1063 padPositions[k] = padPosition;
1065 fitter.AddPoint(&time, padPosition,1);
1069 /////////////////////////////
1071 ////////////////////////////
1072 if(nb3pc < 3) return kFALSE;
1075 fitter.GetParameters(line);
1076 Float_t pointError = -1.0;
1077 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1079 /////////////////////////////////////////////////////
1080 // Now fill the PRF: second loop over clusters
1081 /////////////////////////////////////////////////////
1082 for(Int_t k = 0; k < npoints; k++){
1084 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1085 Short_t *signals = cl->GetSignals(); // signal
1086 Double_t time = cl->GetPadTime(); // time bin
1087 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1088 Float_t padPos = cl->GetPadCol(); // middle pad
1089 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1090 Float_t ycenter = 0.0; // relative center charge
1091 Float_t ymin = 0.0; // relative left charge
1092 Float_t ymax = 0.0; // relative right charge
1094 //Requiere simply two pads clusters at least
1095 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1096 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1097 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1098 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1099 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1100 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1104 Int_t rowcl = cl->GetPadRow(); // row of cluster
1105 Int_t colcl = cl->GetPadCol(); // col of cluster
1106 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1107 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1108 Float_t xcl = cl->GetY(); // y cluster
1109 Float_t qcl = cl->GetQ(); // charge cluster
1110 Int_t plane = GetPlane(detector); // plane
1111 Int_t chamber = GetChamber(detector); // chamber
1112 Double_t xdiff = dpad; // reconstructed position constant
1113 Double_t x = dpad; // reconstructed position moved
1114 Float_t ep = pointError; // error of fit
1115 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1116 Float_t signal3 = (Float_t)signals[3]; // signal
1117 Float_t signal2 = (Float_t)signals[2]; // signal
1118 Float_t signal4 = (Float_t)signals[4]; // signal
1119 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1121 //////////////////////////////
1123 /////////////////////////////
1124 if(fDebugLevel > 0){
1125 if ( !fDebugStreamer ) {
1127 TDirectory *backup = gDirectory;
1128 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1129 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1135 Float_t y = ycenter;
1136 (* fDebugStreamer) << "HandlePRFtracklet"<<
1137 "caligroup="<<caligroup<<
1138 "detector="<<detector<<
1140 "chamber="<<chamber<<
1141 "npoints="<<npoints<<
1150 "padPosition="<<padPositions[k]<<
1151 "padPosTracklet="<<padPosTracklet<<
1156 "signal1="<<signal1<<
1157 "signal2="<<signal2<<
1158 "signal3="<<signal3<<
1159 "signal4="<<signal4<<
1160 "signal5="<<signal5<<
1166 (* fDebugStreamer) << "HandlePRFtracklet"<<
1167 "caligroup="<<caligroup<<
1168 "detector="<<detector<<
1170 "chamber="<<chamber<<
1171 "npoints="<<npoints<<
1180 "padPosition="<<padPositions[k]<<
1181 "padPosTracklet="<<padPosTracklet<<
1186 "signal1="<<signal1<<
1187 "signal2="<<signal2<<
1188 "signal3="<<signal3<<
1189 "signal4="<<signal4<<
1190 "signal5="<<signal5<<
1196 (* fDebugStreamer) << "HandlePRFtracklet"<<
1197 "caligroup="<<caligroup<<
1198 "detector="<<detector<<
1200 "chamber="<<chamber<<
1201 "npoints="<<npoints<<
1210 "padPosition="<<padPositions[k]<<
1211 "padPosTracklet="<<padPosTracklet<<
1216 "signal1="<<signal1<<
1217 "signal2="<<signal2<<
1218 "signal3="<<signal3<<
1219 "signal4="<<signal4<<
1220 "signal5="<<signal5<<
1226 ////////////////////////////
1228 ///////////////////////////
1229 if(npoints < fNumberClusters) continue;
1230 if(npoints > fNumberClustersf) continue;
1231 if(nb3pc <= 5) continue;
1232 if((time >= 21) || (time < 7)) continue;
1233 if(TMath::Abs(snp) >= 1.0) continue;
1234 if(TMath::Abs(qcl) < 80) continue;
1236 ////////////////////////////
1238 ///////////////////////////
1240 if(TMath::Abs(dpad) < 1.5) {
1241 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1242 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1244 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1245 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1246 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1248 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1249 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1250 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1254 if(TMath::Abs(dpad) < 1.5) {
1255 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1256 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1258 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1259 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1260 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1262 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1263 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1264 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1271 //____________Offine tracking in the AliTRDtracker_____________________________
1272 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1275 // PRF width calibration
1276 // Assume a Gaussian shape: determinate the position of the three pad clusters
1277 // Fit with a straight line
1278 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1279 // Fill the PRF as function of angle of the track
1283 ///////////////////////////////////////////
1284 // Take the parameters of the track
1285 //////////////////////////////////////////
1286 // take now the snp, tnp and tgl from the track
1287 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1288 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1289 if( TMath::Abs(snp) < 1.){
1290 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1292 Double_t tgl = tracklet->GetTgl(); // dz/dl
1293 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1295 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1296 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1297 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1298 // at the end with correction due to linear fit
1299 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1300 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1302 ///////////////////////////////
1303 // Calculate tnp group shift
1304 ///////////////////////////////
1305 Bool_t echec = kFALSE;
1306 Double_t shift = 0.0;
1307 //Calculate the shift in x coresponding to this tnp
1308 if(fNgroupprf != 0.0){
1309 shift = -3.0*(fNgroupprf-1)-1.5;
1310 Double_t limithigh = -0.2*(fNgroupprf-1);
1311 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1313 while(tnp > limithigh){
1319 // do nothing if out of tnp range
1320 if(echec) return kFALSE;
1322 ///////////////////////
1324 //////////////////////
1326 Int_t nb3pc = 0; // number of three pads clusters used for fit
1327 Double_t *padPositions; // to see the difference between the fit and the 3 pad clusters position
1328 padPositions = new Double_t[AliTRDseed::knTimebins];
1329 for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1330 padPositions[k] = 0.0;
1332 TLinearFitter fitter(2,"pol1"); // TLinearFitter for the linear fit in the drift region
1333 fitter.StoreData(kFALSE);
1334 fitter.ClearPoints();
1336 ////////////////////////////
1337 // loop over the clusters
1338 ////////////////////////////
1339 AliTRDcluster *cl = 0x0;
1340 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1341 if(!(cl = tracklet->GetClusters(ic))) continue;
1343 Double_t time = cl->GetPadTime();
1344 if((time<=7) || (time>=21)) continue;
1345 Short_t *signals = cl->GetSignals();
1346 Float_t xcenter = 0.0;
1347 Bool_t echec = kTRUE;
1349 /////////////////////////////////////////////////////////////
1350 // Center 3 balanced: position with the center of the pad
1351 /////////////////////////////////////////////////////////////
1352 if ((((Float_t) signals[3]) > 0.0) &&
1353 (((Float_t) signals[2]) > 0.0) &&
1354 (((Float_t) signals[4]) > 0.0)) {
1356 // Security if the denomiateur is 0
1357 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1358 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1359 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1360 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1361 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1367 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1370 ////////////////////////////////////////////////////////
1371 //if no echec: calculate with the position of the pad
1372 // Position of the cluster
1373 // fill the linear fitter
1374 ///////////////////////////////////////////////////////
1375 Double_t padPosition = xcenter + cl->GetPadCol();
1376 padPositions[ic] = padPosition;
1378 fitter.AddPoint(&time, padPosition,1);
1384 //////////////////////////////
1385 // fit with a straight line
1386 /////////////////////////////
1387 if(nb3pc < 3) return kFALSE;
1390 fitter.GetParameters(line);
1391 Float_t pointError = -1.0;
1392 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1394 ////////////////////////////////////////////////
1395 // Fill the PRF: Second loop over clusters
1396 //////////////////////////////////////////////
1397 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1398 if(!(cl = tracklet->GetClusters(ic))) continue;
1400 Short_t *signals = cl->GetSignals(); // signal
1401 Double_t time = cl->GetPadTime(); // time bin
1402 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1403 Float_t padPos = cl->GetPadCol(); // middle pad
1404 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1405 Float_t ycenter = 0.0; // relative center charge
1406 Float_t ymin = 0.0; // relative left charge
1407 Float_t ymax = 0.0; // relative right charge
1409 ////////////////////////////////////////////////////////////////
1410 // Calculate ycenter, ymin and ymax even for two pad clusters
1411 ////////////////////////////////////////////////////////////////
1412 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1413 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1414 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1415 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1416 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1417 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1420 /////////////////////////
1421 // Calibration group
1422 ////////////////////////
1423 Int_t rowcl = cl->GetPadRow(); // row of cluster
1424 Int_t colcl = cl->GetPadCol(); // col of cluster
1425 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1426 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1427 Float_t xcl = cl->GetY(); // y cluster
1428 Float_t qcl = cl->GetQ(); // charge cluster
1429 Int_t plane = GetPlane(fDetectorPreviousTrack); // plane
1430 Int_t chamber = GetChamber(fDetectorPreviousTrack); // chamber
1431 Double_t xdiff = dpad; // reconstructed position constant
1432 Double_t x = dpad; // reconstructed position moved
1433 Float_t ep = pointError; // error of fit
1434 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1435 Float_t signal3 = (Float_t)signals[3]; // signal
1436 Float_t signal2 = (Float_t)signals[2]; // signal
1437 Float_t signal4 = (Float_t)signals[4]; // signal
1438 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1440 /////////////////////
1442 ////////////////////
1444 if(fDebugLevel > 0){
1445 if ( !fDebugStreamer ) {
1447 TDirectory *backup = gDirectory;
1448 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1449 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1454 Float_t y = ycenter;
1455 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1456 "caligroup="<<caligroup<<
1457 "detector="<<fDetectorPreviousTrack<<
1459 "chamber="<<chamber<<
1460 "npoints="<<nbclusters<<
1469 "padPosition="<<padPositions[ic]<<
1470 "padPosTracklet="<<padPosTracklet<<
1475 "signal1="<<signal1<<
1476 "signal2="<<signal2<<
1477 "signal3="<<signal3<<
1478 "signal4="<<signal4<<
1479 "signal5="<<signal5<<
1485 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1486 "caligroup="<<caligroup<<
1487 "detector="<<fDetectorPreviousTrack<<
1489 "chamber="<<chamber<<
1490 "npoints="<<nbclusters<<
1499 "padPosition="<<padPositions[ic]<<
1500 "padPosTracklet="<<padPosTracklet<<
1505 "signal1="<<signal1<<
1506 "signal2="<<signal2<<
1507 "signal3="<<signal3<<
1508 "signal4="<<signal4<<
1509 "signal5="<<signal5<<
1515 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1516 "caligroup="<<caligroup<<
1517 "detector="<<fDetectorPreviousTrack<<
1519 "chamber="<<chamber<<
1520 "npoints="<<nbclusters<<
1529 "padPosition="<<padPositions[ic]<<
1530 "padPosTracklet="<<padPosTracklet<<
1535 "signal1="<<signal1<<
1536 "signal2="<<signal2<<
1537 "signal3="<<signal3<<
1538 "signal4="<<signal4<<
1539 "signal5="<<signal5<<
1545 /////////////////////
1547 /////////////////////
1548 if(nbclusters < fNumberClusters) continue;
1549 if(nbclusters > fNumberClustersf) continue;
1550 if(nb3pc <= 5) continue;
1551 if((time >= 21) || (time < 7)) continue;
1552 if(TMath::Abs(qcl) < 80) continue;
1553 if( TMath::Abs(snp) > 1.) continue;
1556 ////////////////////////
1558 ///////////////////////
1560 if(TMath::Abs(dpad) < 1.5) {
1561 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1562 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1563 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1565 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1566 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1567 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1569 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1570 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1571 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1576 if(TMath::Abs(dpad) < 1.5) {
1577 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1578 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1580 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1581 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1582 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1584 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1585 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1586 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1589 } // second loop over clusters
1595 ///////////////////////////////////////////////////////////////////////////////////////
1596 // Pad row col stuff: see if masked or not
1597 ///////////////////////////////////////////////////////////////////////////////////////
1598 //_____________________________________________________________________________
1599 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t row, Int_t col)
1602 // See if we are not near a masked pad
1605 if (!IsPadOn(detector, col, row)) {
1606 fGoodTracklet = kFALSE;
1610 if (!IsPadOn(detector, col-1, row)) {
1611 fGoodTracklet = kFALSE;
1616 if (!IsPadOn(detector, col+1, row)) {
1617 fGoodTracklet = kFALSE;
1622 //_____________________________________________________________________________
1623 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1626 // Look in the choosen database if the pad is On.
1627 // If no the track will be "not good"
1630 // Get the parameter object
1631 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1633 AliInfo("Could not get calibDB");
1637 if (!cal->IsChamberInstalled(detector) ||
1638 cal->IsChamberMasked(detector) ||
1639 cal->IsPadMasked(detector,col,row)) {
1647 ///////////////////////////////////////////////////////////////////////////////////////
1648 // Calibration groups: calculate the number of groups, localise...
1649 ////////////////////////////////////////////////////////////////////////////////////////
1650 //_____________________________________________________________________________
1651 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1654 // Calculate the calibration group number for i
1657 // Row of the cluster and position in the pad groups
1659 if (fCalibraMode->GetNnZ(i) != 0) {
1660 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1664 // Col of the cluster and position in the pad groups
1666 if (fCalibraMode->GetNnRphi(i) != 0) {
1667 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1670 return posc*fCalibraMode->GetNfragZ(i)+posr;
1673 //____________________________________________________________________________________
1674 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1677 // Calculate the total number of calibration groups
1681 fCalibraMode->ModePadCalibration(2,i);
1682 fCalibraMode->ModePadFragmentation(0,2,0,i);
1683 fCalibraMode->SetDetChamb2(i);
1684 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1685 fCalibraMode->ModePadCalibration(0,i);
1686 fCalibraMode->ModePadFragmentation(0,0,0,i);
1687 fCalibraMode->SetDetChamb0(i);
1688 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1689 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1693 //_____________________________________________________________________________
1694 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1697 // Set the mode of calibration group in the z direction for the parameter i
1702 fCalibraMode->SetNz(i, Nz);
1705 AliInfo("You have to choose between 0 and 4");
1709 //_____________________________________________________________________________
1710 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1713 // Set the mode of calibration group in the rphi direction for the parameter i
1718 fCalibraMode->SetNrphi(i ,Nrphi);
1721 AliInfo("You have to choose between 0 and 6");
1725 //____________Set the pad calibration variables for the detector_______________
1726 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1729 // For the detector calcul the first Xbins and set the number of row
1730 // and col pads per calibration groups, the number of calibration
1731 // groups in the detector.
1734 // first Xbins of the detector
1736 fCalibraMode->CalculXBins(detector,0);
1739 fCalibraMode->CalculXBins(detector,1);
1742 fCalibraMode->CalculXBins(detector,2);
1745 // fragmentation of idect
1746 for (Int_t i = 0; i < 3; i++) {
1747 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1748 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1749 , (Int_t) GetChamber(detector)
1750 , (Int_t) GetSector(detector),i);
1756 //_____________________________________________________________________________
1757 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1760 // Should be between 0 and 6
1763 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1764 AliInfo("The number of groups must be between 0 and 6!");
1767 fNgroupprf = numberGroupsPRF;
1771 ///////////////////////////////////////////////////////////////////////////////////////////
1772 // Per tracklet: store or reset the info, fill the histos with the info
1773 //////////////////////////////////////////////////////////////////////////////////////////
1774 //_____________________________________________________________________________
1775 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1778 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1779 // Correct from the gain correction before
1782 // time bin of the cluster not corrected
1783 Int_t time = cl->GetPadTime();
1785 //Correct for the gain coefficient used in the database for reconstruction
1786 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1787 Float_t correction = 1.0;
1788 Float_t normalisation = 6.67;
1789 // we divide with gain in AliTRDclusterizer::Transform...
1790 if( correctthegain > 0 ) normalisation /= correctthegain;
1793 // take dd/dl corrected from the angle of the track
1794 correction = dqdl / (normalisation);
1797 // Fill the fAmpTotal with the charge
1798 if (fCH2dOn && cl->IsInChamber()) {
1799 fAmpTotal[(Int_t) group[0]] += correction;
1802 // Fill the fPHPlace and value
1804 fPHPlace[time] = group[1];
1805 fPHValue[time] = correction;
1809 //____________Offine tracking in the AliTRDtracker_____________________________
1810 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1813 // Reset values per tracklet
1816 //Reset good tracklet
1817 fGoodTracklet = kTRUE;
1819 // Reset the fPHValue
1821 //Reset the fPHValue and fPHPlace
1822 for (Int_t k = 0; k < fTimeMax; k++) {
1828 // Reset the fAmpTotal where we put value
1830 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1835 //____________Offine tracking in the AliTRDtracker_____________________________
1836 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1839 // For the offline tracking or mcm tracklets
1840 // This function will be called in the functions UpdateHistogram...
1841 // to fill the info of a track for the relativ gain calibration
1844 Int_t nb = 0; // Nombre de zones traversees
1845 Int_t fd = -1; // Premiere zone non nulle
1846 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1848 if(nbclusters < fNumberClusters) return;
1849 if(nbclusters > fNumberClustersf) return;
1852 // See if the track goes through different zones
1853 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1854 if (fAmpTotal[k] > 0.0) {
1855 totalcharge += fAmpTotal[k];
1868 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1870 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1871 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1874 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1878 if ((fAmpTotal[fd] > 0.0) &&
1879 (fAmpTotal[fd+1] > 0.0)) {
1880 // One of the two very big
1881 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1883 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1884 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1887 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1890 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1892 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1894 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1895 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1898 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1901 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1904 if (fCalibraMode->GetNfragZ(0) > 1) {
1905 if (fAmpTotal[fd] > 0.0) {
1906 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1907 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1908 // One of the two very big
1909 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1911 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1912 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1915 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1918 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1920 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1922 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1923 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1926 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1928 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1939 //____________Offine tracking in the AliTRDtracker_____________________________
1940 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1943 // For the offline tracking or mcm tracklets
1944 // This function will be called in the functions UpdateHistogram...
1945 // to fill the info of a track for the drift velocity calibration
1948 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1949 Int_t fd1 = -1; // Premiere zone non nulle
1950 Int_t fd2 = -1; // Deuxieme zone non nulle
1951 Int_t k1 = -1; // Debut de la premiere zone
1952 Int_t k2 = -1; // Debut de la seconde zone
1953 Int_t nbclusters = 0; // number of clusters
1957 // See if the track goes through different zones
1958 for (Int_t k = 0; k < fTimeMax; k++) {
1959 if (fPHValue[k] > 0.0) {
1965 if (fPHPlace[k] != fd1) {
1971 if (fPHPlace[k] != fd2) {
1978 if(nbclusters < fNumberClusters) return;
1979 if(nbclusters > fNumberClustersf) return;
1985 for (Int_t i = 0; i < fTimeMax; i++) {
1987 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1988 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1991 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1996 if ((fd1 == fd2+1) ||
1998 // One of the two fast all the think
1999 if (k2 > (k1+fDifference)) {
2000 //we choose to fill the fd1 with all the values
2002 for (Int_t i = 0; i < fTimeMax; i++) {
2004 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2007 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2011 if ((k2+fDifference) < fTimeMax) {
2012 //we choose to fill the fd2 with all the values
2014 for (Int_t i = 0; i < fTimeMax; i++) {
2016 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2019 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2024 // Two zones voisines sinon rien!
2025 if (fCalibraMode->GetNfragZ(1) > 1) {
2027 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2028 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2029 // One of the two fast all the think
2030 if (k2 > (k1+fDifference)) {
2031 //we choose to fill the fd1 with all the values
2033 for (Int_t i = 0; i < fTimeMax; i++) {
2035 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2038 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2042 if ((k2+fDifference) < fTimeMax) {
2043 //we choose to fill the fd2 with all the values
2045 for (Int_t i = 0; i < fTimeMax; i++) {
2047 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2050 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2056 // Two zones voisines sinon rien!
2058 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2059 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2060 // One of the two fast all the think
2061 if (k2 > (k1 + fDifference)) {
2062 //we choose to fill the fd1 with all the values
2064 for (Int_t i = 0; i < fTimeMax; i++) {
2066 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2069 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2073 if ((k2+fDifference) < fTimeMax) {
2074 //we choose to fill the fd2 with all the values
2076 for (Int_t i = 0; i < fTimeMax; i++) {
2078 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2081 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2092 //////////////////////////////////////////////////////////////////////////////////////////
2093 // DAQ process functions
2094 /////////////////////////////////////////////////////////////////////////////////////////
2095 //_____________________________________________________________________
2096 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamTB *rawStream, Bool_t nocheck)
2099 // Event Processing loop - AliTRDrawStreamTB
2100 // TestBeam 2007 version
2101 // 0 timebin problem
2105 // Same algorithm as TestBeam but different reader
2108 Int_t withInput = 1;
2110 Double_t phvalue[16][144][36];
2111 for(Int_t k = 0; k < 36; k++){
2112 for(Int_t j = 0; j < 16; j++){
2113 for(Int_t c = 0; c < 144; c++){
2114 phvalue[j][c][k] = 0.0;
2119 fDetectorPreviousTrack = -1;
2122 Int_t nbtimebin = 0;
2130 while (rawStream->Next()) {
2132 Int_t idetector = rawStream->GetDet(); // current detector
2133 Int_t imcm = rawStream->GetMCM(); // current MCM
2134 Int_t irob = rawStream->GetROB(); // current ROB
2136 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2139 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2143 for(Int_t k = 0; k < 36; k++){
2144 for(Int_t j = 0; j < 16; j++){
2145 for(Int_t c = 0; c < 144; c++){
2146 phvalue[j][c][k] = 0.0;
2152 fDetectorPreviousTrack = idetector;
2153 fMCMPrevious = imcm;
2154 fROBPrevious = irob;
2156 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2157 if(nbtimebin == 0) return 0;
2158 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2159 fTimeMax = nbtimebin;
2161 baseline = rawStream->GetCommonAdditive(); // common additive baseline
2163 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2164 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2165 Int_t col = rawStream->GetCol();
2166 Int_t row = rawStream->GetRow();
2169 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2172 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2174 for(Int_t itime = iTimeBin; itime < fin; itime++){
2175 phvalue[row][col][itime] = signal[n]-baseline;
2180 // fill the last one
2181 if(fDetectorPreviousTrack != -1){
2184 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2187 for(Int_t k = 0; k < 36; k++){
2188 for(Int_t j = 0; j < 16; j++){
2189 for(Int_t c = 0; c < 144; c++){
2190 phvalue[j][c][k] = 0.0;
2199 while (rawStream->Next()) {
2201 Int_t idetector = rawStream->GetDet(); // current detector
2202 Int_t imcm = rawStream->GetMCM(); // current MCM
2203 Int_t irob = rawStream->GetROB(); // current ROB
2205 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2208 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2211 for(Int_t k = 0; k < 36; k++){
2212 for(Int_t j = 0; j < 16; j++){
2213 for(Int_t c = 0; c < 144; c++){
2214 phvalue[j][c][k] = 0.0;
2220 fDetectorPreviousTrack = idetector;
2221 fMCMPrevious = imcm;
2222 fROBPrevious = irob;
2224 baseline = rawStream->GetCommonAdditive(); // common baseline
2226 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2227 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2228 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2229 Int_t col = rawStream->GetCol();
2230 Int_t row = rawStream->GetRow();
2232 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2235 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2237 for(Int_t itime = iTimeBin; itime < fin; itime++){
2238 phvalue[row][col][itime] = signal[n]-baseline;
2243 // fill the last one
2244 if(fDetectorPreviousTrack != -1){
2247 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2250 for(Int_t k = 0; k < 36; k++){
2251 for(Int_t j = 0; j < 16; j++){
2252 for(Int_t c = 0; c < 144; c++){
2253 phvalue[j][c][k] = 0.0;
2263 //_____________________________________________________________________
2264 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliTRDrawStreamTB *rawStream, Bool_t nocheck)
2267 // Event Processing loop - AliTRDrawStreamTB
2268 // Use old AliTRDmcmtracklet code
2269 // 0 timebin problem
2273 // Algorithm with mcm tracklet
2276 Int_t withInput = 1;
2278 AliTRDmcm mcm = AliTRDmcm(0);
2279 AliTRDtrigParam *param = AliTRDtrigParam::Instance();
2280 rawStream->SetSharedPadReadout(kTRUE);
2282 fDetectorPreviousTrack = -1;
2286 Int_t nbtimebin = 0;
2294 while (rawStream->Next()) {
2296 Int_t idetector = rawStream->GetDet(); // current detector
2297 Int_t imcm = rawStream->GetMCM(); // current MCM
2298 Int_t irob = rawStream->GetROB(); // current ROB
2299 row = rawStream->GetRow();
2302 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2305 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2311 mcm.SetChaId(idetector);
2313 mcm.SetColRange(0,21);
2316 if(fDetectorPreviousTrack == -1){
2319 mcm.SetChaId(idetector);
2321 mcm.SetColRange(0,21);
2325 fDetectorPreviousTrack = idetector;
2326 fMCMPrevious = imcm;
2327 fROBPrevious = irob;
2329 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2330 if(nbtimebin == 0) return 0;
2331 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2332 fTimeMax = nbtimebin;
2333 param->SetTimeRange(0,fTimeMax);
2335 baseline = rawStream->GetCommonAdditive(); // common additive baseline
2337 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2338 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2339 Int_t adc = rawStream->GetADC();
2342 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2345 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2347 for(Int_t itime = iTimeBin; itime < fin; itime++){
2348 mcm.SetADC(adc,itime,(signal[n]-baseline));
2353 // fill the last one
2354 if(fDetectorPreviousTrack != -1){
2357 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2362 mcm.SetRobId(fROBPrevious);
2363 mcm.SetChaId(fDetectorPreviousTrack);
2365 mcm.SetColRange(0,21);
2372 while (rawStream->Next()) {
2374 Int_t idetector = rawStream->GetDet(); // current detector
2375 Int_t imcm = rawStream->GetMCM(); // current MCM
2376 Int_t irob = rawStream->GetROB(); // current ROB
2377 row = rawStream->GetRow();
2379 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2382 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2388 mcm.SetChaId(idetector);
2390 mcm.SetColRange(0,21);
2394 fDetectorPreviousTrack = idetector;
2395 fMCMPrevious = imcm;
2396 fROBPrevious = irob;
2398 baseline = rawStream->GetCommonAdditive(); // common baseline
2400 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2401 param->SetTimeRange(0,fTimeMax);
2402 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2403 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2404 Int_t adc = rawStream->GetADC();
2407 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2410 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2412 for(Int_t itime = iTimeBin; itime < fin; itime++){
2413 mcm.SetADC(adc,itime,(signal[n]-baseline));
2418 // fill the last one
2419 if(fDetectorPreviousTrack != -1){
2422 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2426 mcm.SetRobId(fROBPrevious);
2427 mcm.SetChaId(fDetectorPreviousTrack);
2429 mcm.SetColRange(0,21);
2437 //_____________________________________________________________________
2438 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2441 // Event processing loop - AliRawReader
2442 // Testbeam 2007 version
2445 AliTRDrawStreamTB rawStream(rawReader);
2447 rawReader->Select("TRD");
2449 return ProcessEventDAQ(&rawStream, nocheck);
2452 //_________________________________________________________________________
2453 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2455 eventHeaderStruct *event,
2458 eventHeaderStruct* /*event*/,
2465 // process date event
2466 // Testbeam 2007 version
2469 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2470 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2474 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2479 //_____________________________________________________________________
2480 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliRawReader *rawReader, Bool_t nocheck)
2483 // Event processing loop - AliRawReader
2484 // use the old mcm traklet code
2487 AliTRDrawStreamTB rawStream(rawReader);
2489 rawReader->Select("TRD");
2491 return ProcessEventDAQV1(&rawStream, nocheck);
2493 //_________________________________________________________________________
2494 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(
2496 eventHeaderStruct *event,
2499 eventHeaderStruct* /*event*/,
2506 // process date event
2507 // use the old mcm tracklet code
2510 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2511 Int_t result=ProcessEventDAQV1(rawReader, nocheck);
2515 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2520 //////////////////////////////////////////////////////////////////////////////
2521 // Routine inside the DAQ process
2522 /////////////////////////////////////////////////////////////////////////////
2523 //_______________________________________________________________________
2524 Int_t AliTRDCalibraFillHisto::FillDAQ(AliTRDmcm *mcm){
2527 // Return 2 if some tracklets are found and used, 1 if nothing
2534 for (Int_t iSeed = 0; iSeed < 4; iSeed++) {
2536 if (mcm->GetSeedCol()[iSeed] < 0) {
2540 nbev += TestTracklet(mcm->GetChaId(),mcm->GetRow(),iSeed,mcm);
2545 if(nbev > 0) nbev = 2;
2551 //__________________________________________________________________________
2552 Int_t AliTRDCalibraFillHisto::TestTracklet( Int_t idet, Int_t row, Int_t iSeed, AliTRDmcm *mcm){
2555 // Build the tracklet and return if the tracklet if finally used or not (1/0)
2560 AliTRDmcmTracklet mcmtracklet = AliTRDmcmTracklet();
2561 //mcmtracklet.Reset();
2562 mcmtracklet.SetDetector(idet);
2563 mcmtracklet.SetRow(row);
2564 mcmtracklet.SetN(0);
2566 Int_t iCol, iCol1, iCol2, track[3];
2567 iCol = mcm->GetSeedCol()[iSeed]; // 0....20 (MCM)
2568 mcm->GetColRange(iCol1,iCol2); // range in the pad plane
2571 for (Int_t iTime = 0; iTime < fTimeMax; iTime++) {
2573 amp[0] = mcm->GetADC(iCol-1,iTime);
2574 amp[1] = mcm->GetADC(iCol ,iTime);
2575 amp[2] = mcm->GetADC(iCol+1,iTime);
2577 if(mcm->IsCluster(iCol,iTime)) {
2579 mcmtracklet.AddCluster(iCol+iCol1,iTime,amp,track);
2582 else if ((iCol+1+1) < 21) {
2584 amp[0] = mcm->GetADC(iCol-1+1,iTime);
2585 amp[1] = mcm->GetADC(iCol +1,iTime);
2586 amp[2] = mcm->GetADC(iCol+1+1,iTime);
2588 if(mcm->IsCluster(iCol+1,iTime)) {
2590 mcmtracklet.AddCluster(iCol+1+iCol1,iTime,amp,track);
2599 nbev = UpdateHistogramcm(&mcmtracklet);
2604 //____________Online trackling in AliTRDtrigger________________________________
2605 Int_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2608 // Return if the tracklet is finally used or not (1/0) for calibration
2613 fGoodTracklet = kTRUE;
2615 // Localisation of the Xbins involved
2616 Int_t idect = trk->GetDetector();
2617 Int_t idectrue = trk->GetDetector();
2620 Int_t nbclusters = trk->GetNclusters();
2622 // Eventuelle correction due to track angle in z direction
2623 Float_t correction = 1.0;
2624 if (fMcmCorrectAngle) {
2625 Float_t z = trk->GetRowz();
2626 Float_t r = trk->GetTime0();
2627 correction = r / TMath::Sqrt((r*r+z*z));
2631 Int_t row = trk->GetRow();
2634 // Boucle sur les clusters
2635 // Condition on number of cluster: don't come from the middle of the detector
2638 for(Int_t k =0; k < 36; k++) amph[k]=0.0;
2639 Double_t ampTotal = 0.0;
2641 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2643 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2644 Int_t time = trk->GetClusterTime(icl);
2645 Int_t col = trk->GetClusterCol(icl);
2647 CheckGoodTracklet(idect,row,col);
2649 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2650 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2651 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2653 ampTotal += (Float_t) (amp[0]+amp[1]+amp[2]);
2654 amph[time]=amp[0]+amp[1]+amp[2];
2656 if(fDebugLevel > 0){
2657 if ( !fDebugStreamer ) {
2659 TDirectory *backup = gDirectory;
2660 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2661 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2664 Double_t amp0 = amp[0];
2665 Double_t amp1 = amp[1];
2666 Double_t amp2 = amp[2];
2668 (* fDebugStreamer) << "UpdateHistogramcm0"<<
2669 "nbclusters="<<nbclusters<<
2676 "detector="<<idectrue<<
2680 } // Boucle clusters
2682 if((amph[0] > 100.0) || (!fGoodTracklet) || (trk->GetNclusters() < fNumberClusters) || (trk->GetNclusters() > fNumberClustersf)) used = 0;
2685 for(Int_t k = 0; k < fTimeMax; k++) UpdateDAQ(idect,0,0,k,amph[k],fTimeMax);
2686 //((TH2I *)GetCH2d()->Fill(ampTotal/30.0,idect));
2690 if(fDebugLevel > 0){
2691 if ( !fDebugStreamer ) {
2693 TDirectory *backup = gDirectory;
2694 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2695 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2698 Double_t amph0 = amph[0];
2699 Double_t amphlast = amph[fTimeMax-1];
2700 Double_t rms = TMath::RMS(fTimeMax,amph);
2701 Int_t goodtracklet = (Int_t) fGoodTracklet;
2703 (* fDebugStreamer) << "UpdateHistogramcm1"<<
2704 "nbclusters="<<nbclusters<<
2705 "ampTotal="<<ampTotal<<
2707 "detector="<<idectrue<<
2709 "amphlast="<<amphlast<<
2710 "goodtracklet="<<goodtracklet<<
2718 //_______________________________________________________________________
2719 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2722 // Look for the maximum by collapsing over the time
2723 // Sum over four pad col and two pad row
2729 Int_t idect = fDetectorPreviousTrack;
2731 for(Int_t tb = 0; tb < 36; tb++){
2735 fGoodTracklet = kTRUE;
2736 fDetectorPreviousTrack = 0;
2739 ///////////////////////////
2741 /////////////////////////
2745 Double_t integralMax = -1;
2747 for (Int_t ir = 1; ir <= 15; ir++)
2749 for (Int_t ic = 2; ic <= 142; ic++)
2751 Double_t integral = 0;
2752 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2754 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2756 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2757 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2760 for(Int_t tb = 0; tb< fTimeMax; tb++){
2761 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2767 if (integralMax < integral)
2771 integralMax = integral;
2772 } // check max integral
2776 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2778 if((imaxRow == 0) || (imaxCol == 0)) used=1;
2779 CheckGoodTracklet(fDetectorPreviousTrack,imaxRow,imaxCol);
2780 if(!fGoodTracklet) used = 1;;
2782 // /////////////////////////////////////////////////////
2783 // sum ober 2 row and 4 pad cols for each time bins
2784 // ////////////////////////////////////////////////////
2787 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2789 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2791 for(Int_t it = 0; it < fTimeMax; it++){
2792 sum[it] += phvalue[ir][ic][it];
2798 Double_t sumcharge = 0.0;
2799 for(Int_t it = 0; it < fTimeMax; it++){
2800 sumcharge += sum[it];
2801 if(sum[it] > 20.0) nbcl++;
2805 /////////////////////////////////////////////////////////
2807 ////////////////////////////////////////////////////////
2808 if(fDebugLevel > 0){
2809 if ( !fDebugStreamer ) {
2811 TDirectory *backup = gDirectory;
2812 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2813 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2816 Double_t amph0 = sum[0];
2817 Double_t amphlast = sum[fTimeMax-1];
2818 Double_t rms = TMath::RMS(fTimeMax,sum);
2819 Int_t goodtracklet = (Int_t) fGoodTracklet;
2820 for(Int_t it = 0; it < fTimeMax; it++){
2821 Double_t clustera = sum[it];
2823 (* fDebugStreamer) << "FillDAQa"<<
2824 "ampTotal="<<sumcharge<<
2827 "detector="<<idect<<
2829 "amphlast="<<amphlast<<
2830 "goodtracklet="<<goodtracklet<<
2831 "clustera="<<clustera<<
2838 ////////////////////////////////////////////////////////
2840 ///////////////////////////////////////////////////////
2841 if(sum[0] > 100.0) used = 1;
2842 if(nbcl < fNumberClusters) used = 1;
2843 if(nbcl > fNumberClustersf) used = 1;
2845 //if(fDetectorPreviousTrack == 15){
2846 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2848 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2850 for(Int_t it = 0; it < fTimeMax; it++){
2851 UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2855 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2862 //____________Online trackling in AliTRDtrigger________________________________
2863 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2867 // Fill a simple average pulse height
2871 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2877 //____________Write_____________________________________________________
2878 //_____________________________________________________________________
2879 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2882 // Write infos to file
2886 if ( fDebugStreamer ) {
2887 delete fDebugStreamer;
2888 fDebugStreamer = 0x0;
2891 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2896 ,fNumberUsedPh[1]));
2898 TDirectory *backup = gDirectory;
2904 option = "recreate";
2906 TFile f(filename,option.Data());
2908 TStopwatch stopwatch;
2911 f.WriteTObject(fCalibraVector);
2916 f.WriteTObject(fCH2d);
2921 f.WriteTObject(fPH2d);
2926 f.WriteTObject(fPRF2d);
2929 if(fLinearFitterOn){
2930 AnalyseLinearFitter();
2931 f.WriteTObject(fLinearVdriftFit);
2936 if ( backup ) backup->cd();
2938 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2939 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2941 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2943 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2944 //___________________________________________probe the histos__________________________________________________
2945 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2948 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2949 // debug mode with 2 for TH2I and 3 for TProfile2D
2950 // It gives a pointer to a Double_t[7] with the info following...
2951 // [0] : number of calibration groups with entries
2952 // [1] : minimal number of entries found
2953 // [2] : calibration group number of the min
2954 // [3] : maximal number of entries found
2955 // [4] : calibration group number of the max
2956 // [5] : mean number of entries found
2957 // [6] : mean relative error
2960 Double_t *info = new Double_t[7];
2962 // Number of Xbins (detectors or groups of pads)
2963 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2964 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2967 Double_t nbwe = 0; //number of calibration groups with entries
2968 Double_t minentries = 0; //minimal number of entries found
2969 Double_t maxentries = 0; //maximal number of entries found
2970 Double_t placemin = 0; //calibration group number of the min
2971 Double_t placemax = -1; //calibration group number of the max
2972 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2973 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2975 Double_t counter = 0;
2978 TH1F *nbEntries = 0x0;//distribution of the number of entries
2979 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2980 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2982 // Beginning of the loop over the calibration groups
2983 for (Int_t idect = 0; idect < nbins; idect++) {
2985 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2986 projch->SetDirectory(0);
2988 // Number of entries for this calibration group
2989 Double_t nentries = 0.0;
2991 for (Int_t k = 0; k < nxbins; k++) {
2992 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2996 for (Int_t k = 0; k < nxbins; k++) {
2997 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2998 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2999 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3007 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3008 nbEntries = new TH1F("Number of entries","Number of entries"
3009 ,100,(Int_t)nentries/2,nentries*2);
3010 nbEntries->SetDirectory(0);
3011 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3013 nbEntriesPerGroup->SetDirectory(0);
3014 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3015 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3016 nbEntriesPerSp->SetDirectory(0);
3019 if(nentries > 0) nbEntries->Fill(nentries);
3020 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3021 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3026 if(nentries > maxentries){
3027 maxentries = nentries;
3031 minentries = nentries;
3033 if(nentries < minentries){
3034 minentries = nentries;
3040 meanstats += nentries;
3042 }//calibration groups loop
3044 if(nbwe > 0) meanstats /= nbwe;
3045 if(counter > 0) meanrelativerror /= counter;
3047 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3048 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3049 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3050 AliInfo(Form("The mean number of entries is %f",meanstats));
3051 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3054 info[1] = minentries;
3056 info[3] = maxentries;
3058 info[5] = meanstats;
3059 info[6] = meanrelativerror;
3062 gStyle->SetPalette(1);
3063 gStyle->SetOptStat(1111);
3064 gStyle->SetPadBorderMode(0);
3065 gStyle->SetCanvasColor(10);
3066 gStyle->SetPadLeftMargin(0.13);
3067 gStyle->SetPadRightMargin(0.01);
3068 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3071 nbEntries->Draw("");
3073 nbEntriesPerSp->SetStats(0);
3074 nbEntriesPerSp->Draw("");
3075 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3077 nbEntriesPerGroup->SetStats(0);
3078 nbEntriesPerGroup->Draw("");
3084 //____________________________________________________________________________
3085 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3088 // Return a Int_t[4] with:
3089 // 0 Mean number of entries
3090 // 1 median of number of entries
3091 // 2 rms of number of entries
3092 // 3 number of group with entries
3095 Double_t *stat = new Double_t[4];
3098 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3099 Double_t *weight = new Double_t[nbofgroups];
3100 Int_t *nonul = new Int_t[nbofgroups];
3102 for(Int_t k = 0; k < nbofgroups; k++){
3103 if(fEntriesCH[k] > 0) {
3105 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3108 else weight[k] = 0.0;
3110 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3111 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3112 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3117 //____________________________________________________________________________
3118 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3121 // Return a Int_t[4] with:
3122 // 0 Mean number of entries
3123 // 1 median of number of entries
3124 // 2 rms of number of entries
3125 // 3 number of group with entries
3128 Double_t *stat = new Double_t[4];
3131 Int_t nbofgroups = 540;
3132 Double_t *weight = new Double_t[nbofgroups];
3133 Int_t *nonul = new Int_t[nbofgroups];
3135 for(Int_t k = 0; k < nbofgroups; k++){
3136 if(fEntriesLinearFitter[k] > 0) {
3138 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3141 else weight[k] = 0.0;
3143 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3144 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3145 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3150 //////////////////////////////////////////////////////////////////////////////////////
3152 //////////////////////////////////////////////////////////////////////////////////////
3153 //_____________________________________________________________________________
3154 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3157 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3158 // If fNgroupprf is zero then no binning in tnp
3162 name += fCalibraMode->GetNz(2);
3164 name += fCalibraMode->GetNrphi(2);
3168 if(fNgroupprf != 0){
3170 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3171 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3172 fPRF2d->SetYTitle("Det/pad groups");
3173 fPRF2d->SetXTitle("Position x/W [pad width units]");
3174 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3175 fPRF2d->SetStats(0);
3178 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3179 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3180 fPRF2d->SetYTitle("Det/pad groups");
3181 fPRF2d->SetXTitle("Position x/W [pad width units]");
3182 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3183 fPRF2d->SetStats(0);
3188 //_____________________________________________________________________________
3189 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3192 // Create the 2D histos
3196 name += fCalibraMode->GetNz(1);
3198 name += fCalibraMode->GetNrphi(1);
3200 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3201 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3203 fPH2d->SetYTitle("Det/pad groups");
3204 fPH2d->SetXTitle("time [#mus]");
3205 fPH2d->SetZTitle("<PH> [a.u.]");
3209 //_____________________________________________________________________________
3210 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3213 // Create the 2D histos
3217 name += fCalibraMode->GetNz(0);
3219 name += fCalibraMode->GetNrphi(0);
3221 fCH2d = new TH2I("CH2d",(const Char_t *) name
3222 ,fNumberBinCharge,0,300,nn,0,nn);
3223 fCH2d->SetYTitle("Det/pad groups");
3224 fCH2d->SetXTitle("charge deposit [a.u]");
3225 fCH2d->SetZTitle("counts");
3230 //////////////////////////////////////////////////////////////////////////////////
3231 // Set relative scale
3232 /////////////////////////////////////////////////////////////////////////////////
3233 //_____________________________________________________________________________
3234 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3237 // Set the factor that will divide the deposited charge
3238 // to fit in the histo range [0,300]
3241 if (RelativeScale > 0.0) {
3242 fRelativeScale = RelativeScale;
3245 AliInfo("RelativeScale must be strict positif!");
3249 //////////////////////////////////////////////////////////////////////////////////
3250 // Quick way to fill a histo
3251 //////////////////////////////////////////////////////////////////////////////////
3252 //_____________________________________________________________________
3253 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3256 // FillCH2d: Marian style
3259 //skip simply the value out of range
3260 if((y>=300.0) || (y<0.0)) return;
3262 //Calcul the y place
3263 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3264 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3267 fCH2d->GetArray()[place]++;
3271 //////////////////////////////////////////////////////////////////////////////////
3272 // Geometrical functions
3273 ///////////////////////////////////////////////////////////////////////////////////
3274 //_____________________________________________________________________________
3275 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3278 // Reconstruct the plane number from the detector number
3281 return ((Int_t) (d % 6));
3285 //_____________________________________________________________________________
3286 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3289 // Reconstruct the chamber number from the detector number
3293 return ((Int_t) (d % 30) / fgkNplan);
3297 //_____________________________________________________________________________
3298 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3301 // Reconstruct the sector number from the detector number
3305 return ((Int_t) (d / fg));
3308 ///////////////////////////////////////////////////////////////////////////////////
3309 // Getter functions for DAQ of the CH2d and the PH2d
3310 //////////////////////////////////////////////////////////////////////////////////
3311 //_____________________________________________________________________
3312 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3315 // return pointer to fPH2d TProfile2D
3316 // create a new TProfile2D if it doesn't exist allready
3322 fTimeMax = nbtimebin;
3323 fSf = samplefrequency;
3329 //_____________________________________________________________________
3330 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3333 // return pointer to fCH2d TH2I
3334 // create a new TH2I if it doesn't exist allready
3343 ////////////////////////////////////////////////////////////////////////////////////////////
3344 // Drift velocity calibration
3345 ///////////////////////////////////////////////////////////////////////////////////////////
3346 //_____________________________________________________________________
3347 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3350 // return pointer to TLinearFitter Calibration
3351 // if force is true create a new TLinearFitter if it doesn't exist allready
3354 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3355 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3358 // if we are forced and TLinearFitter doesn't yet exist create it
3360 // new TLinearFitter
3361 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3362 fLinearFitterArray.AddAt(linearfitter,detector);
3363 return linearfitter;
3366 //____________________________________________________________________________
3367 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3370 // Analyse array of linear fitter because can not be written
3371 // Store two arrays: one with the param the other one with the error param + number of entries
3374 for(Int_t k = 0; k < 540; k++){
3375 TLinearFitter *linearfitter = GetLinearFitter(k);
3376 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3377 TVectorD *par = new TVectorD(2);
3378 TVectorD pare = TVectorD(2);
3379 TVectorD *parE = new TVectorD(3);
3380 linearfitter->Eval();
3381 linearfitter->GetParameters(*par);
3382 linearfitter->GetErrors(pare);
3383 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3384 (*parE)[0] = pare[0]*ppointError;
3385 (*parE)[1] = pare[1]*ppointError;
3386 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3387 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3388 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);