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 "AliTRDrawStreamBase.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 ,fLimitChargeIntegration(kFALSE)
134 ,fCalibraMode(new AliTRDCalibraMode())
137 ,fDetectorPreviousTrack(-1)
141 ,fNumberClustersf(30)
147 ,fNumberBinCharge(100)
153 ,fGoodTracklet(kTRUE)
155 ,fEntriesLinearFitter(0x0)
160 ,fLinearFitterArray(540)
161 ,fLinearVdriftFit(0x0)
166 // Default constructor
170 // Init some default values
173 fNumberUsedCh[0] = 0;
174 fNumberUsedCh[1] = 0;
175 fNumberUsedPh[0] = 0;
176 fNumberUsedPh[1] = 0;
178 fGeo = new AliTRDgeometry();
182 //______________________________________________________________________________________
183 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
186 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
189 ,fPRF2dOn(c.fPRF2dOn)
190 ,fHisto2d(c.fHisto2d)
191 ,fVector2d(c.fVector2d)
192 ,fLinearFitterOn(c.fLinearFitterOn)
193 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
194 ,fRelativeScale(c.fRelativeScale)
195 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
196 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
199 ,fDebugLevel(c.fDebugLevel)
200 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
201 ,fMCMPrevious(c.fMCMPrevious)
202 ,fROBPrevious(c.fROBPrevious)
203 ,fNumberClusters(c.fNumberClusters)
204 ,fNumberClustersf(c.fNumberClustersf)
205 ,fProcent(c.fProcent)
206 ,fDifference(c.fDifference)
207 ,fNumberTrack(c.fNumberTrack)
208 ,fTimeMax(c.fTimeMax)
210 ,fNumberBinCharge(c.fNumberBinCharge)
211 ,fNumberBinPRF(c.fNumberBinPRF)
212 ,fNgroupprf(c.fNgroupprf)
216 ,fGoodTracklet(c.fGoodTracklet)
218 ,fEntriesLinearFitter(0x0)
223 ,fLinearFitterArray(540)
224 ,fLinearVdriftFit(0x0)
231 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
232 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
234 fPH2d = new TProfile2D(*c.fPH2d);
235 fPH2d->SetDirectory(0);
238 fPRF2d = new TProfile2D(*c.fPRF2d);
239 fPRF2d->SetDirectory(0);
242 fCH2d = new TH2I(*c.fCH2d);
243 fCH2d->SetDirectory(0);
245 if(c.fLinearVdriftFit){
246 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
249 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
250 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
255 fGeo = new AliTRDgeometry();
258 //____________________________________________________________________________________
259 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
262 // AliTRDCalibraFillHisto destructor
266 if ( fDebugStreamer ) delete fDebugStreamer;
268 if ( fCalDetGain ) delete fCalDetGain;
269 if ( fCalROCGain ) delete fCalROCGain;
276 //_____________________________________________________________________________
277 void AliTRDCalibraFillHisto::Destroy()
288 //_____________________________________________________________________________
289 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
292 // Delete DebugStreamer
295 if ( fDebugStreamer ) delete fDebugStreamer;
298 //_____________________________________________________________________________
299 void AliTRDCalibraFillHisto::ClearHistos()
319 //////////////////////////////////////////////////////////////////////////////////
320 // calibration with AliTRDtrackV1: Init, Update
321 //////////////////////////////////////////////////////////////////////////////////
322 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
323 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
326 // Init the histograms and stuff to be filled
331 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
333 AliInfo("Could not get calibDB");
336 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
338 AliInfo("Could not get CommonParam");
343 fTimeMax = cal->GetNumberOfTimeBins();
344 fSf = parCom->GetSamplingFrequency();
347 //calib object from database used for reconstruction
348 if(fCalDetGain) delete fCalDetGain;
349 fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
351 // Calcul Xbins Chambd0, Chamb2
352 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
353 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
354 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
356 // If vector method On initialised all the stuff
358 fCalibraVector = new AliTRDCalibraVector();
359 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
360 fCalibraVector->SetTimeMax(fTimeMax);
361 if(fNgroupprf != 0) {
362 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
363 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
366 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
367 fCalibraVector->SetPRFRange(1.5);
369 for(Int_t k = 0; k < 3; k++){
370 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
371 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
373 TString namech("Nz");
374 namech += fCalibraMode->GetNz(0);
376 namech += fCalibraMode->GetNrphi(0);
377 fCalibraVector->SetNameCH((const char* ) namech);
378 TString nameph("Nz");
379 nameph += fCalibraMode->GetNz(1);
381 nameph += fCalibraMode->GetNrphi(1);
382 fCalibraVector->SetNamePH((const char* ) nameph);
383 TString nameprf("Nz");
384 nameprf += fCalibraMode->GetNz(2);
386 nameprf += fCalibraMode->GetNrphi(2);
388 nameprf += fNgroupprf;
389 fCalibraVector->SetNamePRF((const char* ) nameprf);
392 // Create the 2D histos corresponding to the pad groupCalibration mode
395 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
396 ,fCalibraMode->GetNz(0)
397 ,fCalibraMode->GetNrphi(0)));
399 // Create the 2D histo
404 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
405 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
409 fEntriesCH = new Int_t[ntotal0];
410 for(Int_t k = 0; k < ntotal0; k++){
417 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
418 ,fCalibraMode->GetNz(1)
419 ,fCalibraMode->GetNrphi(1)));
421 // Create the 2D histo
426 fPHPlace = new Short_t[fTimeMax];
427 for (Int_t k = 0; k < fTimeMax; k++) {
430 fPHValue = new Float_t[fTimeMax];
431 for (Int_t k = 0; k < fTimeMax; k++) {
435 if (fLinearFitterOn) {
436 //fLinearFitterArray.Expand(540);
437 fLinearFitterArray.SetName("ArrayLinearFitters");
438 fEntriesLinearFitter = new Int_t[540];
439 for(Int_t k = 0; k < 540; k++){
440 fEntriesLinearFitter[k] = 0;
442 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
447 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
448 ,fCalibraMode->GetNz(2)
449 ,fCalibraMode->GetNrphi(2)));
450 // Create the 2D histo
452 CreatePRF2d(ntotal2);
459 //____________Offline tracking in the AliTRDtracker____________________________
460 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
463 // Use AliTRDtrack for the calibration
467 AliTRDcluster *cl = 0x0; // pointeur to cluster
468 Int_t index0 = 0; // index of the first cluster in the current chamber
469 Int_t index1 = 0; // index of the current cluster in the current chamber
471 //////////////////////////////
472 // loop over the clusters
473 ///////////////////////////////
474 while((cl = t->GetCluster(index1))){
476 /////////////////////////////////////////////////////////////////////////
477 // Fill the infos for the previous clusters if not the same detector
478 ////////////////////////////////////////////////////////////////////////
479 Int_t detector = cl->GetDetector();
480 if ((detector != fDetectorPreviousTrack) &&
481 (index0 != index1)) {
485 //If the same track, then look if the previous detector is in
486 //the same plane, if yes: not a good track
487 if ((GetPlane(detector) == GetPlane(fDetectorPreviousTrack))) {
491 // Fill only if the track doesn't touch a masked pad or doesn't
494 // drift velocity unables to cut bad tracklets
495 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
499 FillTheInfoOfTheTrackCH(index1-index0);
504 FillTheInfoOfTheTrackPH();
507 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
510 } // if a good tracklet
513 ResetfVariablestracklet();
516 } // Fill at the end the charge
519 /////////////////////////////////////////////////////////////
520 // Localise and take the calib gain object for the detector
521 ////////////////////////////////////////////////////////////
522 if (detector != fDetectorPreviousTrack) {
524 //Localise the detector
525 LocalisationDetectorXbins(detector);
528 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
530 AliInfo("Could not get calibDB");
535 if( fCalROCGain ) delete fCalROCGain;
536 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
540 // Reset the detectbjobsor
541 fDetectorPreviousTrack = detector;
543 ///////////////////////////////////////
544 // Store the info of the cluster
545 ///////////////////////////////////////
546 Int_t row = cl->GetPadRow();
547 Int_t col = cl->GetPadCol();
548 CheckGoodTracklet(detector,row,col);
549 Int_t group[2] = {0,0};
550 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
551 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
552 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
556 } // while on clusters
558 ///////////////////////////
559 // Fill the last plane
560 //////////////////////////
561 if( index0 != index1 ){
567 // drift velocity unables to cut bad tracklets
568 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
572 FillTheInfoOfTheTrackCH(index1-index0);
577 FillTheInfoOfTheTrackPH();
580 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
582 } // if a good tracklet
587 ResetfVariablestracklet();
592 //____________Offline tracking in the AliTRDtracker____________________________
593 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
596 // Use AliTRDtrackV1 for the calibration
600 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
601 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
604 ///////////////////////////
605 // loop over the tracklet
606 ///////////////////////////
607 for(Int_t itr = 0; itr < 6; itr++){
609 if(!(tracklet = t->GetTracklet(itr))) continue;
610 if(!tracklet->IsOK()) continue;
612 ResetfVariablestracklet();
614 //////////////////////////////////////////
615 // localisation of the tracklet and dqdl
616 //////////////////////////////////////////
617 Int_t plane = tracklet->GetPlane();
619 while(!(cl = tracklet->GetClusters(ic++))) continue;
620 Int_t detector = cl->GetDetector();
621 if (detector != fDetectorPreviousTrack) {
622 // don't use the rest of this track if in the same plane
623 if ((plane == GetPlane(fDetectorPreviousTrack))) {
626 //Localise the detector bin
627 LocalisationDetectorXbins(detector);
629 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
631 AliInfo("Could not get calibDB");
635 if( fCalROCGain ) delete fCalROCGain;
636 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
638 fDetectorPreviousTrack = detector;
642 ////////////////////////////
643 // loop over the clusters
644 ////////////////////////////
645 Int_t nbclusters = 0;
646 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
647 if(!(cl = tracklet->GetClusters(ic))) continue;
650 // Store the info bis of the tracklet
651 Int_t row = cl->GetPadRow();
652 Int_t col = cl->GetPadCol();
653 CheckGoodTracklet(detector,row,col);
654 Int_t group[2] = {0,0};
655 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
656 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
657 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(ic),group,row,col);
660 ////////////////////////////////////////
661 // Fill the stuffs if a good tracklet
662 ////////////////////////////////////////
665 // drift velocity unables to cut bad tracklets
666 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
670 FillTheInfoOfTheTrackCH(nbclusters);
675 FillTheInfoOfTheTrackPH();
678 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
680 } // if a good tracklet
688 ///////////////////////////////////////////////////////////////////////////////////
689 // Routine inside the update with AliTRDtrack
690 ///////////////////////////////////////////////////////////////////////////////////
691 //____________Offine tracking in the AliTRDtracker_____________________________
692 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
695 // Drift velocity calibration:
696 // Fit the clusters with a straight line
697 // From the slope find the drift velocity
700 //Number of points: if less than 3 return kFALSE
701 Int_t npoints = index1-index0;
702 if(npoints <= 2) return kFALSE;
707 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
708 // parameters of the track
709 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in the plan yx track
710 Double_t tgl = t->GetTglPlane(GetPlane(detector)); // dz/dl and not dz/dx!
711 Double_t tnp = 0.0; // tan angle in the plan xy track
712 if( TMath::Abs(snp) < 1.){
713 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
715 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
716 // tilting pad and cross row
717 Int_t crossrow = 0; // if it crosses a pad row
718 Int_t rowp = -1; // if it crosses a pad row
719 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
720 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
721 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
723 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
724 Double_t dydt = 0.0; // dydt tracklet after straight line fit
725 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
726 Double_t pointError = 0.0; // error after straight line fit
727 Int_t nbli = 0; // number linear fitter points
728 linearFitterTracklet.StoreData(kFALSE);
729 linearFitterTracklet.ClearPoints();
731 //////////////////////////////
732 // loop over clusters
733 ////////////////////////////
734 for(Int_t k = 0; k < npoints; k++){
736 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
737 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
738 Double_t ycluster = cl->GetY();
739 Int_t time = cl->GetPadTime();
740 Double_t timeis = time/fSf;
741 //Double_t q = cl->GetQ();
742 //See if cross two pad rows
743 Int_t row = cl->GetPadRow();
745 if(row != rowp) crossrow = 1;
747 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
752 //////////////////////////////
754 //////////////////////////////
755 if(nbli <= 2) return kFALSE;
757 linearFitterTracklet.Eval();
758 linearFitterTracklet.GetParameters(pars);
759 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
760 errorpar = linearFitterTracklet.GetParError(1)*pointError;
763 /////////////////////////////
765 ////////////////////////////
767 if ( !fDebugStreamer ) {
769 TDirectory *backup = gDirectory;
770 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
771 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
775 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
776 //"snpright="<<snpright<<
777 "npoints="<<npoints<<
779 "detector="<<detector<<
786 "crossrow="<<crossrow<<
787 "errorpar="<<errorpar<<
788 "pointError="<<pointError<<
792 Int_t nbclusters = index1-index0;
793 Int_t plane = GetPlane(fDetectorPreviousTrack);
795 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
796 //"snpright="<<snpright<<
797 "nbclusters="<<nbclusters<<
798 "detector="<<fDetectorPreviousTrack<<
804 ///////////////////////////
806 ///////////////////////////
807 if(npoints < fNumberClusters) return kFALSE;
808 if(npoints > fNumberClustersf) return kFALSE;
809 if(pointError >= 0.1) return kFALSE;
810 if(crossrow == 1) return kFALSE;
812 ////////////////////////////
814 ////////////////////////////
816 //Add to the linear fitter of the detector
817 if( TMath::Abs(snp) < 1.){
818 Double_t x = tnp-dzdx*tnt;
819 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
820 if(fLinearFitterDebugOn) {
821 fLinearVdriftFit->Update(detector,x,pars[1]);
823 fEntriesLinearFitter[detector]++;
829 //____________Offine tracking in the AliTRDtracker_____________________________
830 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
833 // Drift velocity calibration:
834 // Fit the clusters with a straight line
835 // From the slope find the drift velocity
838 ////////////////////////////////////////////////
839 //Number of points: if less than 3 return kFALSE
840 /////////////////////////////////////////////////
841 if(nbclusters <= 2) return kFALSE;
846 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
847 // results of the linear fit
848 Double_t dydt = 0.0; // dydt tracklet after straight line fit
849 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
850 Double_t pointError = 0.0; // error after straight line fit
851 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
852 Int_t crossrow = 0; // if it crosses a pad row
853 Int_t rowp = -1; // if it crosses a pad row
854 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
856 linearFitterTracklet.StoreData(kFALSE);
857 linearFitterTracklet.ClearPoints();
859 ///////////////////////////////////////////
860 // Take the parameters of the track
861 //////////////////////////////////////////
862 // take now the snp, tnp and tgl from the track
863 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
864 Double_t tnp = 0.0; // dy/dx at the end of the chamber
865 if( TMath::Abs(snp) < 1.){
866 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
868 Double_t tgl = tracklet->GetTgl(); // dz/dl
869 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
871 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
872 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
873 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
874 // at the end with correction due to linear fit
875 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
876 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
879 ////////////////////////////
880 // loop over the clusters
881 ////////////////////////////
883 AliTRDcluster *cl = 0x0;
884 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
885 if(!(cl = tracklet->GetClusters(ic))) continue;
886 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
888 Double_t ycluster = cl->GetY();
889 Int_t time = cl->GetPadTime();
890 Double_t timeis = time/fSf;
891 //See if cross two pad rows
892 Int_t row = cl->GetPadRow();
893 if(rowp==-1) rowp = row;
894 if(row != rowp) crossrow = 1;
896 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
901 ////////////////////////////////////
902 // Do the straight line fit now
903 ///////////////////////////////////
905 linearFitterTracklet.Eval();
906 linearFitterTracklet.GetParameters(pars);
907 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
908 errorpar = linearFitterTracklet.GetParError(1)*pointError;
911 ////////////////////////////////
913 ///////////////////////////////
917 if ( !fDebugStreamer ) {
919 TDirectory *backup = gDirectory;
920 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
921 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
925 Int_t plane = GetPlane(fDetectorPreviousTrack);
927 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
928 //"snpright="<<snpright<<
929 "nbclusters="<<nbclusters<<
930 "detector="<<fDetectorPreviousTrack<<
938 "crossrow="<<crossrow<<
939 "errorpar="<<errorpar<<
940 "pointError="<<pointError<<
945 /////////////////////////
947 ////////////////////////
949 if(nbclusters < fNumberClusters) return kFALSE;
950 if(nbclusters > fNumberClustersf) return kFALSE;
951 if(pointError >= 0.1) return kFALSE;
952 if(crossrow == 1) return kFALSE;
954 ///////////////////////
956 //////////////////////
959 //Add to the linear fitter of the detector
960 if( TMath::Abs(snp) < 1.){
961 Double_t x = tnp-dzdx*tnt;
962 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
963 if(fLinearFitterDebugOn) {
964 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
966 fEntriesLinearFitter[fDetectorPreviousTrack]++;
972 //____________Offine tracking in the AliTRDtracker_____________________________
973 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
976 // PRF width calibration
977 // Assume a Gaussian shape: determinate the position of the three pad clusters
978 // Fit with a straight line
979 // Take the fitted values for all the clusters (3 or 2 pad clusters)
980 // Fill the PRF as function of angle of the track
985 //////////////////////////
987 /////////////////////////
988 Int_t npoints = index1-index0; // number of total points
989 Int_t nb3pc = 0; // number of three pads clusters used for fit
990 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
991 // To see the difference due to the fit
992 Double_t *padPositions;
993 padPositions = new Double_t[npoints];
994 for(Int_t k = 0; k < npoints; k++){
995 padPositions[k] = 0.0;
997 // Take the tgl and snp with the track t now
998 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
999 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
1000 Float_t dzdx = 0.0; // dzdx
1002 if(TMath::Abs(snp) < 1.0){
1003 tnp = snp / (TMath::Sqrt(1-snp*snp));
1004 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1007 TLinearFitter fitter(2,"pol1");
1008 fitter.StoreData(kFALSE);
1009 fitter.ClearPoints();
1012 ///////////////////////////
1013 // calcul the tnp group
1014 ///////////////////////////
1015 Bool_t echec = kFALSE;
1016 Double_t shift = 0.0;
1017 //Calculate the shift in x coresponding to this tnp
1018 if(fNgroupprf != 0.0){
1019 shift = -3.0*(fNgroupprf-1)-1.5;
1020 Double_t limithigh = -0.2*(fNgroupprf-1);
1021 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1023 while(tnp > limithigh){
1029 if(echec) return kFALSE;
1032 //////////////////////
1034 /////////////////////
1035 for(Int_t k = 0; k < npoints; k++){
1037 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1038 Short_t *signals = cl->GetSignals();
1039 Double_t time = cl->GetPadTime();
1040 //Calculate x if possible
1041 Float_t xcenter = 0.0;
1042 Bool_t echec = kTRUE;
1043 if((time<=7) || (time>=21)) continue;
1044 // Center 3 balanced: position with the center of the pad
1045 if ((((Float_t) signals[3]) > 0.0) &&
1046 (((Float_t) signals[2]) > 0.0) &&
1047 (((Float_t) signals[4]) > 0.0)) {
1049 // Security if the denomiateur is 0
1050 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1051 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1052 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1053 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1054 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1060 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1062 //if no echec: calculate with the position of the pad
1063 // Position of the cluster
1064 Double_t padPosition = xcenter + cl->GetPadCol();
1065 padPositions[k] = padPosition;
1067 fitter.AddPoint(&time, padPosition,1);
1071 /////////////////////////////
1073 ////////////////////////////
1074 if(nb3pc < 3) return kFALSE;
1077 fitter.GetParameters(line);
1078 Float_t pointError = -1.0;
1079 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1081 /////////////////////////////////////////////////////
1082 // Now fill the PRF: second loop over clusters
1083 /////////////////////////////////////////////////////
1084 for(Int_t k = 0; k < npoints; k++){
1086 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1087 Short_t *signals = cl->GetSignals(); // signal
1088 Double_t time = cl->GetPadTime(); // time bin
1089 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1090 Float_t padPos = cl->GetPadCol(); // middle pad
1091 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1092 Float_t ycenter = 0.0; // relative center charge
1093 Float_t ymin = 0.0; // relative left charge
1094 Float_t ymax = 0.0; // relative right charge
1096 //Requiere simply two pads clusters at least
1097 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1098 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1099 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1100 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1101 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1102 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1106 Int_t rowcl = cl->GetPadRow(); // row of cluster
1107 Int_t colcl = cl->GetPadCol(); // col of cluster
1108 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1109 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1110 Float_t xcl = cl->GetY(); // y cluster
1111 Float_t qcl = cl->GetQ(); // charge cluster
1112 Int_t plane = GetPlane(detector); // plane
1113 Int_t chamber = GetChamber(detector); // chamber
1114 Double_t xdiff = dpad; // reconstructed position constant
1115 Double_t x = dpad; // reconstructed position moved
1116 Float_t ep = pointError; // error of fit
1117 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1118 Float_t signal3 = (Float_t)signals[3]; // signal
1119 Float_t signal2 = (Float_t)signals[2]; // signal
1120 Float_t signal4 = (Float_t)signals[4]; // signal
1121 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1123 //////////////////////////////
1125 /////////////////////////////
1126 if(fDebugLevel > 0){
1127 if ( !fDebugStreamer ) {
1129 TDirectory *backup = gDirectory;
1130 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1131 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1137 Float_t y = ycenter;
1138 (* fDebugStreamer) << "HandlePRFtracklet"<<
1139 "caligroup="<<caligroup<<
1140 "detector="<<detector<<
1142 "chamber="<<chamber<<
1143 "npoints="<<npoints<<
1152 "padPosition="<<padPositions[k]<<
1153 "padPosTracklet="<<padPosTracklet<<
1158 "signal1="<<signal1<<
1159 "signal2="<<signal2<<
1160 "signal3="<<signal3<<
1161 "signal4="<<signal4<<
1162 "signal5="<<signal5<<
1168 (* fDebugStreamer) << "HandlePRFtracklet"<<
1169 "caligroup="<<caligroup<<
1170 "detector="<<detector<<
1172 "chamber="<<chamber<<
1173 "npoints="<<npoints<<
1182 "padPosition="<<padPositions[k]<<
1183 "padPosTracklet="<<padPosTracklet<<
1188 "signal1="<<signal1<<
1189 "signal2="<<signal2<<
1190 "signal3="<<signal3<<
1191 "signal4="<<signal4<<
1192 "signal5="<<signal5<<
1198 (* fDebugStreamer) << "HandlePRFtracklet"<<
1199 "caligroup="<<caligroup<<
1200 "detector="<<detector<<
1202 "chamber="<<chamber<<
1203 "npoints="<<npoints<<
1212 "padPosition="<<padPositions[k]<<
1213 "padPosTracklet="<<padPosTracklet<<
1218 "signal1="<<signal1<<
1219 "signal2="<<signal2<<
1220 "signal3="<<signal3<<
1221 "signal4="<<signal4<<
1222 "signal5="<<signal5<<
1228 ////////////////////////////
1230 ///////////////////////////
1231 if(npoints < fNumberClusters) continue;
1232 if(npoints > fNumberClustersf) continue;
1233 if(nb3pc <= 5) continue;
1234 if((time >= 21) || (time < 7)) continue;
1235 if(TMath::Abs(snp) >= 1.0) continue;
1236 if(TMath::Abs(qcl) < 80) continue;
1238 ////////////////////////////
1240 ///////////////////////////
1242 if(TMath::Abs(dpad) < 1.5) {
1243 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1244 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1246 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1247 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1248 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1250 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1251 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1252 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1256 if(TMath::Abs(dpad) < 1.5) {
1257 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1258 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1260 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1261 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1262 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1264 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1265 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1266 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1273 //____________Offine tracking in the AliTRDtracker_____________________________
1274 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1277 // PRF width calibration
1278 // Assume a Gaussian shape: determinate the position of the three pad clusters
1279 // Fit with a straight line
1280 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1281 // Fill the PRF as function of angle of the track
1285 ///////////////////////////////////////////
1286 // Take the parameters of the track
1287 //////////////////////////////////////////
1288 // take now the snp, tnp and tgl from the track
1289 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1290 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1291 if( TMath::Abs(snp) < 1.){
1292 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1294 Double_t tgl = tracklet->GetTgl(); // dz/dl
1295 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1297 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1298 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1299 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1300 // at the end with correction due to linear fit
1301 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1302 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1304 ///////////////////////////////
1305 // Calculate tnp group shift
1306 ///////////////////////////////
1307 Bool_t echec = kFALSE;
1308 Double_t shift = 0.0;
1309 //Calculate the shift in x coresponding to this tnp
1310 if(fNgroupprf != 0.0){
1311 shift = -3.0*(fNgroupprf-1)-1.5;
1312 Double_t limithigh = -0.2*(fNgroupprf-1);
1313 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1315 while(tnp > limithigh){
1321 // do nothing if out of tnp range
1322 if(echec) return kFALSE;
1324 ///////////////////////
1326 //////////////////////
1328 Int_t nb3pc = 0; // number of three pads clusters used for fit
1329 Double_t *padPositions; // to see the difference between the fit and the 3 pad clusters position
1330 padPositions = new Double_t[AliTRDseed::knTimebins];
1331 for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1332 padPositions[k] = 0.0;
1334 TLinearFitter fitter(2,"pol1"); // TLinearFitter for the linear fit in the drift region
1335 fitter.StoreData(kFALSE);
1336 fitter.ClearPoints();
1338 ////////////////////////////
1339 // loop over the clusters
1340 ////////////////////////////
1341 AliTRDcluster *cl = 0x0;
1342 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1343 if(!(cl = tracklet->GetClusters(ic))) continue;
1345 Double_t time = cl->GetPadTime();
1346 if((time<=7) || (time>=21)) continue;
1347 Short_t *signals = cl->GetSignals();
1348 Float_t xcenter = 0.0;
1349 Bool_t echec = kTRUE;
1351 /////////////////////////////////////////////////////////////
1352 // Center 3 balanced: position with the center of the pad
1353 /////////////////////////////////////////////////////////////
1354 if ((((Float_t) signals[3]) > 0.0) &&
1355 (((Float_t) signals[2]) > 0.0) &&
1356 (((Float_t) signals[4]) > 0.0)) {
1358 // Security if the denomiateur is 0
1359 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1360 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1361 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1362 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1363 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1369 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1372 ////////////////////////////////////////////////////////
1373 //if no echec: calculate with the position of the pad
1374 // Position of the cluster
1375 // fill the linear fitter
1376 ///////////////////////////////////////////////////////
1377 Double_t padPosition = xcenter + cl->GetPadCol();
1378 padPositions[ic] = padPosition;
1380 fitter.AddPoint(&time, padPosition,1);
1386 //////////////////////////////
1387 // fit with a straight line
1388 /////////////////////////////
1389 if(nb3pc < 3) return kFALSE;
1392 fitter.GetParameters(line);
1393 Float_t pointError = -1.0;
1394 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1396 ////////////////////////////////////////////////
1397 // Fill the PRF: Second loop over clusters
1398 //////////////////////////////////////////////
1399 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1400 if(!(cl = tracklet->GetClusters(ic))) continue;
1402 Short_t *signals = cl->GetSignals(); // signal
1403 Double_t time = cl->GetPadTime(); // time bin
1404 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1405 Float_t padPos = cl->GetPadCol(); // middle pad
1406 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1407 Float_t ycenter = 0.0; // relative center charge
1408 Float_t ymin = 0.0; // relative left charge
1409 Float_t ymax = 0.0; // relative right charge
1411 ////////////////////////////////////////////////////////////////
1412 // Calculate ycenter, ymin and ymax even for two pad clusters
1413 ////////////////////////////////////////////////////////////////
1414 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1415 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1416 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1417 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1418 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1419 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1422 /////////////////////////
1423 // Calibration group
1424 ////////////////////////
1425 Int_t rowcl = cl->GetPadRow(); // row of cluster
1426 Int_t colcl = cl->GetPadCol(); // col of cluster
1427 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1428 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1429 Float_t xcl = cl->GetY(); // y cluster
1430 Float_t qcl = cl->GetQ(); // charge cluster
1431 Int_t plane = GetPlane(fDetectorPreviousTrack); // plane
1432 Int_t chamber = GetChamber(fDetectorPreviousTrack); // chamber
1433 Double_t xdiff = dpad; // reconstructed position constant
1434 Double_t x = dpad; // reconstructed position moved
1435 Float_t ep = pointError; // error of fit
1436 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1437 Float_t signal3 = (Float_t)signals[3]; // signal
1438 Float_t signal2 = (Float_t)signals[2]; // signal
1439 Float_t signal4 = (Float_t)signals[4]; // signal
1440 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1442 /////////////////////
1444 ////////////////////
1446 if(fDebugLevel > 0){
1447 if ( !fDebugStreamer ) {
1449 TDirectory *backup = gDirectory;
1450 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1451 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1456 Float_t y = ycenter;
1457 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1458 "caligroup="<<caligroup<<
1459 "detector="<<fDetectorPreviousTrack<<
1461 "chamber="<<chamber<<
1462 "npoints="<<nbclusters<<
1471 "padPosition="<<padPositions[ic]<<
1472 "padPosTracklet="<<padPosTracklet<<
1477 "signal1="<<signal1<<
1478 "signal2="<<signal2<<
1479 "signal3="<<signal3<<
1480 "signal4="<<signal4<<
1481 "signal5="<<signal5<<
1487 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1488 "caligroup="<<caligroup<<
1489 "detector="<<fDetectorPreviousTrack<<
1491 "chamber="<<chamber<<
1492 "npoints="<<nbclusters<<
1501 "padPosition="<<padPositions[ic]<<
1502 "padPosTracklet="<<padPosTracklet<<
1507 "signal1="<<signal1<<
1508 "signal2="<<signal2<<
1509 "signal3="<<signal3<<
1510 "signal4="<<signal4<<
1511 "signal5="<<signal5<<
1517 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1518 "caligroup="<<caligroup<<
1519 "detector="<<fDetectorPreviousTrack<<
1521 "chamber="<<chamber<<
1522 "npoints="<<nbclusters<<
1531 "padPosition="<<padPositions[ic]<<
1532 "padPosTracklet="<<padPosTracklet<<
1537 "signal1="<<signal1<<
1538 "signal2="<<signal2<<
1539 "signal3="<<signal3<<
1540 "signal4="<<signal4<<
1541 "signal5="<<signal5<<
1547 /////////////////////
1549 /////////////////////
1550 if(nbclusters < fNumberClusters) continue;
1551 if(nbclusters > fNumberClustersf) continue;
1552 if(nb3pc <= 5) continue;
1553 if((time >= 21) || (time < 7)) continue;
1554 if(TMath::Abs(qcl) < 80) continue;
1555 if( TMath::Abs(snp) > 1.) continue;
1558 ////////////////////////
1560 ///////////////////////
1562 if(TMath::Abs(dpad) < 1.5) {
1563 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1564 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1565 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1567 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1568 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1569 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1571 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1572 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1573 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1578 if(TMath::Abs(dpad) < 1.5) {
1579 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1580 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1582 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1583 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1584 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1586 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1587 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1588 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1591 } // second loop over clusters
1597 ///////////////////////////////////////////////////////////////////////////////////////
1598 // Pad row col stuff: see if masked or not
1599 ///////////////////////////////////////////////////////////////////////////////////////
1600 //_____________________________________________________________________________
1601 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t row, Int_t col)
1604 // See if we are not near a masked pad
1607 if (!IsPadOn(detector, col, row)) {
1608 fGoodTracklet = kFALSE;
1612 if (!IsPadOn(detector, col-1, row)) {
1613 fGoodTracklet = kFALSE;
1618 if (!IsPadOn(detector, col+1, row)) {
1619 fGoodTracklet = kFALSE;
1624 //_____________________________________________________________________________
1625 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1628 // Look in the choosen database if the pad is On.
1629 // If no the track will be "not good"
1632 // Get the parameter object
1633 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1635 AliInfo("Could not get calibDB");
1639 if (!cal->IsChamberInstalled(detector) ||
1640 cal->IsChamberMasked(detector) ||
1641 cal->IsPadMasked(detector,col,row)) {
1649 ///////////////////////////////////////////////////////////////////////////////////////
1650 // Calibration groups: calculate the number of groups, localise...
1651 ////////////////////////////////////////////////////////////////////////////////////////
1652 //_____________________________________________________________________________
1653 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1656 // Calculate the calibration group number for i
1659 // Row of the cluster and position in the pad groups
1661 if (fCalibraMode->GetNnZ(i) != 0) {
1662 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1666 // Col of the cluster and position in the pad groups
1668 if (fCalibraMode->GetNnRphi(i) != 0) {
1669 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1672 return posc*fCalibraMode->GetNfragZ(i)+posr;
1675 //____________________________________________________________________________________
1676 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1679 // Calculate the total number of calibration groups
1683 fCalibraMode->ModePadCalibration(2,i);
1684 fCalibraMode->ModePadFragmentation(0,2,0,i);
1685 fCalibraMode->SetDetChamb2(i);
1686 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1687 fCalibraMode->ModePadCalibration(0,i);
1688 fCalibraMode->ModePadFragmentation(0,0,0,i);
1689 fCalibraMode->SetDetChamb0(i);
1690 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1691 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1695 //_____________________________________________________________________________
1696 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1699 // Set the mode of calibration group in the z direction for the parameter i
1704 fCalibraMode->SetNz(i, Nz);
1707 AliInfo("You have to choose between 0 and 4");
1711 //_____________________________________________________________________________
1712 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1715 // Set the mode of calibration group in the rphi direction for the parameter i
1720 fCalibraMode->SetNrphi(i ,Nrphi);
1723 AliInfo("You have to choose between 0 and 6");
1727 //____________Set the pad calibration variables for the detector_______________
1728 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1731 // For the detector calcul the first Xbins and set the number of row
1732 // and col pads per calibration groups, the number of calibration
1733 // groups in the detector.
1736 // first Xbins of the detector
1738 fCalibraMode->CalculXBins(detector,0);
1741 fCalibraMode->CalculXBins(detector,1);
1744 fCalibraMode->CalculXBins(detector,2);
1747 // fragmentation of idect
1748 for (Int_t i = 0; i < 3; i++) {
1749 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1750 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1751 , (Int_t) GetChamber(detector)
1752 , (Int_t) GetSector(detector),i);
1758 //_____________________________________________________________________________
1759 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1762 // Should be between 0 and 6
1765 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1766 AliInfo("The number of groups must be between 0 and 6!");
1769 fNgroupprf = numberGroupsPRF;
1773 ///////////////////////////////////////////////////////////////////////////////////////////
1774 // Per tracklet: store or reset the info, fill the histos with the info
1775 //////////////////////////////////////////////////////////////////////////////////////////
1776 //_____________________________________________________________________________
1777 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1780 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1781 // Correct from the gain correction before
1784 // time bin of the cluster not corrected
1785 Int_t time = cl->GetPadTime();
1787 //Correct for the gain coefficient used in the database for reconstruction
1788 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1789 Float_t correction = 1.0;
1790 Float_t normalisation = 6.67;
1791 // we divide with gain in AliTRDclusterizer::Transform...
1792 if( correctthegain > 0 ) normalisation /= correctthegain;
1795 // take dd/dl corrected from the angle of the track
1796 correction = dqdl / (normalisation);
1799 // Fill the fAmpTotal with the charge
1801 if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1804 // Fill the fPHPlace and value
1806 fPHPlace[time] = group[1];
1807 fPHValue[time] = correction;
1811 //____________Offine tracking in the AliTRDtracker_____________________________
1812 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1815 // Reset values per tracklet
1818 //Reset good tracklet
1819 fGoodTracklet = kTRUE;
1821 // Reset the fPHValue
1823 //Reset the fPHValue and fPHPlace
1824 for (Int_t k = 0; k < fTimeMax; k++) {
1830 // Reset the fAmpTotal where we put value
1832 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1837 //____________Offine tracking in the AliTRDtracker_____________________________
1838 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1841 // For the offline tracking or mcm tracklets
1842 // This function will be called in the functions UpdateHistogram...
1843 // to fill the info of a track for the relativ gain calibration
1846 Int_t nb = 0; // Nombre de zones traversees
1847 Int_t fd = -1; // Premiere zone non nulle
1848 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1850 if(nbclusters < fNumberClusters) return;
1851 if(nbclusters > fNumberClustersf) return;
1854 // See if the track goes through different zones
1855 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1856 if (fAmpTotal[k] > 0.0) {
1857 totalcharge += fAmpTotal[k];
1870 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1872 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1873 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1876 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1880 if ((fAmpTotal[fd] > 0.0) &&
1881 (fAmpTotal[fd+1] > 0.0)) {
1882 // One of the two very big
1883 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1885 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1886 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1889 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1892 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1894 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1896 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1897 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1900 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1903 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1906 if (fCalibraMode->GetNfragZ(0) > 1) {
1907 if (fAmpTotal[fd] > 0.0) {
1908 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1909 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1910 // One of the two very big
1911 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1913 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1914 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1917 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1920 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1922 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1924 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1925 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1928 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1930 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1941 //____________Offine tracking in the AliTRDtracker_____________________________
1942 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1945 // For the offline tracking or mcm tracklets
1946 // This function will be called in the functions UpdateHistogram...
1947 // to fill the info of a track for the drift velocity calibration
1950 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1951 Int_t fd1 = -1; // Premiere zone non nulle
1952 Int_t fd2 = -1; // Deuxieme zone non nulle
1953 Int_t k1 = -1; // Debut de la premiere zone
1954 Int_t k2 = -1; // Debut de la seconde zone
1955 Int_t nbclusters = 0; // number of clusters
1959 // See if the track goes through different zones
1960 for (Int_t k = 0; k < fTimeMax; k++) {
1961 if (fPHValue[k] > 0.0) {
1967 if (fPHPlace[k] != fd1) {
1973 if (fPHPlace[k] != fd2) {
1980 if(nbclusters < fNumberClusters) return;
1981 if(nbclusters > fNumberClustersf) return;
1987 for (Int_t i = 0; i < fTimeMax; i++) {
1989 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1990 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1993 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1998 if ((fd1 == fd2+1) ||
2000 // One of the two fast all the think
2001 if (k2 > (k1+fDifference)) {
2002 //we choose to fill the fd1 with all the values
2004 for (Int_t i = 0; i < fTimeMax; i++) {
2006 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2009 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2013 if ((k2+fDifference) < fTimeMax) {
2014 //we choose to fill the fd2 with all the values
2016 for (Int_t i = 0; i < fTimeMax; i++) {
2018 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2021 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2026 // Two zones voisines sinon rien!
2027 if (fCalibraMode->GetNfragZ(1) > 1) {
2029 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2030 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2031 // One of the two fast all the think
2032 if (k2 > (k1+fDifference)) {
2033 //we choose to fill the fd1 with all the values
2035 for (Int_t i = 0; i < fTimeMax; i++) {
2037 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2040 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2044 if ((k2+fDifference) < fTimeMax) {
2045 //we choose to fill the fd2 with all the values
2047 for (Int_t i = 0; i < fTimeMax; i++) {
2049 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2052 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2058 // Two zones voisines sinon rien!
2060 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2061 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2062 // One of the two fast all the think
2063 if (k2 > (k1 + fDifference)) {
2064 //we choose to fill the fd1 with all the values
2066 for (Int_t i = 0; i < fTimeMax; i++) {
2068 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2071 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2075 if ((k2+fDifference) < fTimeMax) {
2076 //we choose to fill the fd2 with all the values
2078 for (Int_t i = 0; i < fTimeMax; i++) {
2080 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2083 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2094 //////////////////////////////////////////////////////////////////////////////////////////
2095 // DAQ process functions
2096 /////////////////////////////////////////////////////////////////////////////////////////
2097 //_____________________________________________________________________
2098 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2101 // Event Processing loop - AliTRDrawStreamBase
2102 // TestBeam 2007 version
2103 // 0 timebin problem
2107 // Same algorithm as TestBeam but different reader
2110 Int_t withInput = 1;
2112 Double_t phvalue[16][144][36];
2113 for(Int_t k = 0; k < 36; k++){
2114 for(Int_t j = 0; j < 16; j++){
2115 for(Int_t c = 0; c < 144; c++){
2116 phvalue[j][c][k] = 0.0;
2121 fDetectorPreviousTrack = -1;
2124 Int_t nbtimebin = 0;
2132 while (rawStream->Next()) {
2134 Int_t idetector = rawStream->GetDet(); // current detector
2135 Int_t imcm = rawStream->GetMCM(); // current MCM
2136 Int_t irob = rawStream->GetROB(); // current ROB
2138 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2141 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2145 for(Int_t k = 0; k < 36; k++){
2146 for(Int_t j = 0; j < 16; j++){
2147 for(Int_t c = 0; c < 144; c++){
2148 phvalue[j][c][k] = 0.0;
2154 fDetectorPreviousTrack = idetector;
2155 fMCMPrevious = imcm;
2156 fROBPrevious = irob;
2158 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2159 if(nbtimebin == 0) return 0;
2160 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2161 fTimeMax = nbtimebin;
2163 baseline = rawStream->GetCommonAdditive(); // common additive baseline
2165 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2166 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2167 Int_t col = rawStream->GetCol();
2168 Int_t row = rawStream->GetRow();
2171 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2174 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2176 for(Int_t itime = iTimeBin; itime < fin; itime++){
2177 phvalue[row][col][itime] = signal[n]-baseline;
2182 // fill the last one
2183 if(fDetectorPreviousTrack != -1){
2186 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2189 for(Int_t k = 0; k < 36; k++){
2190 for(Int_t j = 0; j < 16; j++){
2191 for(Int_t c = 0; c < 144; c++){
2192 phvalue[j][c][k] = 0.0;
2201 while (rawStream->Next()) {
2203 Int_t idetector = rawStream->GetDet(); // current detector
2204 Int_t imcm = rawStream->GetMCM(); // current MCM
2205 Int_t irob = rawStream->GetROB(); // current ROB
2207 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2210 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2213 for(Int_t k = 0; k < 36; k++){
2214 for(Int_t j = 0; j < 16; j++){
2215 for(Int_t c = 0; c < 144; c++){
2216 phvalue[j][c][k] = 0.0;
2222 fDetectorPreviousTrack = idetector;
2223 fMCMPrevious = imcm;
2224 fROBPrevious = irob;
2226 baseline = rawStream->GetCommonAdditive(); // common baseline
2228 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2229 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2230 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2231 Int_t col = rawStream->GetCol();
2232 Int_t row = rawStream->GetRow();
2234 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2237 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2239 for(Int_t itime = iTimeBin; itime < fin; itime++){
2240 phvalue[row][col][itime] = signal[n]-baseline;
2245 // fill the last one
2246 if(fDetectorPreviousTrack != -1){
2249 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2252 for(Int_t k = 0; k < 36; k++){
2253 for(Int_t j = 0; j < 16; j++){
2254 for(Int_t c = 0; c < 144; c++){
2255 phvalue[j][c][k] = 0.0;
2265 //_____________________________________________________________________
2266 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2269 // Event Processing loop - AliTRDrawStreamBase
2270 // Use old AliTRDmcmtracklet code
2271 // 0 timebin problem
2275 // Algorithm with mcm tracklet
2278 Int_t withInput = 1;
2280 AliTRDmcm mcm = AliTRDmcm(0);
2281 AliTRDtrigParam *param = AliTRDtrigParam::Instance();
2282 rawStream->SetSharedPadReadout(kTRUE);
2284 fDetectorPreviousTrack = -1;
2288 Int_t nbtimebin = 0;
2296 while (rawStream->Next()) {
2298 Int_t idetector = rawStream->GetDet(); // current detector
2299 Int_t imcm = rawStream->GetMCM(); // current MCM
2300 Int_t irob = rawStream->GetROB(); // current ROB
2301 row = rawStream->GetRow();
2304 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2307 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2313 mcm.SetChaId(idetector);
2315 mcm.SetColRange(0,21);
2318 if(fDetectorPreviousTrack == -1){
2321 mcm.SetChaId(idetector);
2323 mcm.SetColRange(0,21);
2327 fDetectorPreviousTrack = idetector;
2328 fMCMPrevious = imcm;
2329 fROBPrevious = irob;
2331 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2332 if(nbtimebin == 0) return 0;
2333 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2334 fTimeMax = nbtimebin;
2335 param->SetTimeRange(0,fTimeMax);
2337 baseline = rawStream->GetCommonAdditive(); // common additive baseline
2339 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2340 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2341 Int_t adc = rawStream->GetADC();
2344 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2347 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2349 for(Int_t itime = iTimeBin; itime < fin; itime++){
2350 mcm.SetADC(adc,itime,(signal[n]-baseline));
2355 // fill the last one
2356 if(fDetectorPreviousTrack != -1){
2359 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2364 mcm.SetRobId(fROBPrevious);
2365 mcm.SetChaId(fDetectorPreviousTrack);
2367 mcm.SetColRange(0,21);
2374 while (rawStream->Next()) {
2376 Int_t idetector = rawStream->GetDet(); // current detector
2377 Int_t imcm = rawStream->GetMCM(); // current MCM
2378 Int_t irob = rawStream->GetROB(); // current ROB
2379 row = rawStream->GetRow();
2381 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2384 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2390 mcm.SetChaId(idetector);
2392 mcm.SetColRange(0,21);
2396 fDetectorPreviousTrack = idetector;
2397 fMCMPrevious = imcm;
2398 fROBPrevious = irob;
2400 baseline = rawStream->GetCommonAdditive(); // common baseline
2402 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2403 param->SetTimeRange(0,fTimeMax);
2404 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
2405 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2406 Int_t adc = rawStream->GetADC();
2409 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
2412 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2414 for(Int_t itime = iTimeBin; itime < fin; itime++){
2415 mcm.SetADC(adc,itime,(signal[n]-baseline));
2420 // fill the last one
2421 if(fDetectorPreviousTrack != -1){
2424 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2428 mcm.SetRobId(fROBPrevious);
2429 mcm.SetChaId(fDetectorPreviousTrack);
2431 mcm.SetColRange(0,21);
2439 //_____________________________________________________________________
2440 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2443 // Event processing loop - AliRawReader
2444 // Testbeam 2007 version
2447 AliTRDrawStreamBase rawStream(rawReader);
2449 rawReader->Select("TRD");
2451 return ProcessEventDAQ(&rawStream, nocheck);
2454 //_________________________________________________________________________
2455 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2457 eventHeaderStruct *event,
2460 eventHeaderStruct* /*event*/,
2467 // process date event
2468 // Testbeam 2007 version
2471 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2472 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2476 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2481 //_____________________________________________________________________
2482 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliRawReader *rawReader, Bool_t nocheck)
2485 // Event processing loop - AliRawReader
2486 // use the old mcm traklet code
2489 AliTRDrawStreamBase rawStream(rawReader);
2491 rawReader->Select("TRD");
2493 return ProcessEventDAQV1(&rawStream, nocheck);
2495 //_________________________________________________________________________
2496 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(
2498 eventHeaderStruct *event,
2501 eventHeaderStruct* /*event*/,
2508 // process date event
2509 // use the old mcm tracklet code
2512 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2513 Int_t result=ProcessEventDAQV1(rawReader, nocheck);
2517 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2522 //////////////////////////////////////////////////////////////////////////////
2523 // Routine inside the DAQ process
2524 /////////////////////////////////////////////////////////////////////////////
2525 //_______________________________________________________________________
2526 Int_t AliTRDCalibraFillHisto::FillDAQ(AliTRDmcm *mcm){
2529 // Return 2 if some tracklets are found and used, 1 if nothing
2536 for (Int_t iSeed = 0; iSeed < 4; iSeed++) {
2538 if (mcm->GetSeedCol()[iSeed] < 0) {
2542 nbev += TestTracklet(mcm->GetChaId(),mcm->GetRow(),iSeed,mcm);
2547 if(nbev > 0) nbev = 2;
2553 //__________________________________________________________________________
2554 Int_t AliTRDCalibraFillHisto::TestTracklet( Int_t idet, Int_t row, Int_t iSeed, AliTRDmcm *mcm){
2557 // Build the tracklet and return if the tracklet if finally used or not (1/0)
2562 AliTRDmcmTracklet mcmtracklet = AliTRDmcmTracklet();
2563 //mcmtracklet.Reset();
2564 mcmtracklet.SetDetector(idet);
2565 mcmtracklet.SetRow(row);
2566 mcmtracklet.SetN(0);
2568 Int_t iCol, iCol1, iCol2, track[3];
2569 iCol = mcm->GetSeedCol()[iSeed]; // 0....20 (MCM)
2570 mcm->GetColRange(iCol1,iCol2); // range in the pad plane
2573 for (Int_t iTime = 0; iTime < fTimeMax; iTime++) {
2575 amp[0] = mcm->GetADC(iCol-1,iTime);
2576 amp[1] = mcm->GetADC(iCol ,iTime);
2577 amp[2] = mcm->GetADC(iCol+1,iTime);
2579 if(mcm->IsCluster(iCol,iTime)) {
2581 mcmtracklet.AddCluster(iCol+iCol1,iTime,amp,track);
2584 else if ((iCol+1+1) < 21) {
2586 amp[0] = mcm->GetADC(iCol-1+1,iTime);
2587 amp[1] = mcm->GetADC(iCol +1,iTime);
2588 amp[2] = mcm->GetADC(iCol+1+1,iTime);
2590 if(mcm->IsCluster(iCol+1,iTime)) {
2592 mcmtracklet.AddCluster(iCol+1+iCol1,iTime,amp,track);
2601 nbev = UpdateHistogramcm(&mcmtracklet);
2606 //____________Online trackling in AliTRDtrigger________________________________
2607 Int_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2610 // Return if the tracklet is finally used or not (1/0) for calibration
2615 fGoodTracklet = kTRUE;
2617 // Localisation of the Xbins involved
2618 Int_t idect = trk->GetDetector();
2619 Int_t idectrue = trk->GetDetector();
2622 Int_t nbclusters = trk->GetNclusters();
2624 // Eventuelle correction due to track angle in z direction
2625 Float_t correction = 1.0;
2626 if (fMcmCorrectAngle) {
2627 Float_t z = trk->GetRowz();
2628 Float_t r = trk->GetTime0();
2629 correction = r / TMath::Sqrt((r*r+z*z));
2633 Int_t row = trk->GetRow();
2636 // Boucle sur les clusters
2637 // Condition on number of cluster: don't come from the middle of the detector
2640 for(Int_t k =0; k < 36; k++) amph[k]=0.0;
2641 Double_t ampTotal = 0.0;
2643 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2645 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2646 Int_t time = trk->GetClusterTime(icl);
2647 Int_t col = trk->GetClusterCol(icl);
2649 CheckGoodTracklet(idect,row,col);
2651 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2652 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2653 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2655 ampTotal += (Float_t) (amp[0]+amp[1]+amp[2]);
2656 amph[time]=amp[0]+amp[1]+amp[2];
2658 if(fDebugLevel > 0){
2659 if ( !fDebugStreamer ) {
2661 TDirectory *backup = gDirectory;
2662 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2663 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2666 Double_t amp0 = amp[0];
2667 Double_t amp1 = amp[1];
2668 Double_t amp2 = amp[2];
2670 (* fDebugStreamer) << "UpdateHistogramcm0"<<
2671 "nbclusters="<<nbclusters<<
2678 "detector="<<idectrue<<
2682 } // Boucle clusters
2684 if((amph[0] > 100.0) || (!fGoodTracklet) || (trk->GetNclusters() < fNumberClusters) || (trk->GetNclusters() > fNumberClustersf)) used = 0;
2687 for(Int_t k = 0; k < fTimeMax; k++) UpdateDAQ(idect,0,0,k,amph[k],fTimeMax);
2688 //((TH2I *)GetCH2d()->Fill(ampTotal/30.0,idect));
2692 if(fDebugLevel > 0){
2693 if ( !fDebugStreamer ) {
2695 TDirectory *backup = gDirectory;
2696 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2697 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2700 Double_t amph0 = amph[0];
2701 Double_t amphlast = amph[fTimeMax-1];
2702 Double_t rms = TMath::RMS(fTimeMax,amph);
2703 Int_t goodtracklet = (Int_t) fGoodTracklet;
2705 (* fDebugStreamer) << "UpdateHistogramcm1"<<
2706 "nbclusters="<<nbclusters<<
2707 "ampTotal="<<ampTotal<<
2709 "detector="<<idectrue<<
2711 "amphlast="<<amphlast<<
2712 "goodtracklet="<<goodtracklet<<
2720 //_______________________________________________________________________
2721 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2724 // Look for the maximum by collapsing over the time
2725 // Sum over four pad col and two pad row
2731 Int_t idect = fDetectorPreviousTrack;
2733 for(Int_t tb = 0; tb < 36; tb++){
2737 fGoodTracklet = kTRUE;
2738 fDetectorPreviousTrack = 0;
2741 ///////////////////////////
2743 /////////////////////////
2747 Double_t integralMax = -1;
2749 for (Int_t ir = 1; ir <= 15; ir++)
2751 for (Int_t ic = 2; ic <= 142; ic++)
2753 Double_t integral = 0;
2754 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2756 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2758 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2759 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2762 for(Int_t tb = 0; tb< fTimeMax; tb++){
2763 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2769 if (integralMax < integral)
2773 integralMax = integral;
2774 } // check max integral
2778 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2780 if((imaxRow == 0) || (imaxCol == 0)) used=1;
2781 CheckGoodTracklet(fDetectorPreviousTrack,imaxRow,imaxCol);
2782 if(!fGoodTracklet) used = 1;;
2784 // /////////////////////////////////////////////////////
2785 // sum ober 2 row and 4 pad cols for each time bins
2786 // ////////////////////////////////////////////////////
2789 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2791 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2793 for(Int_t it = 0; it < fTimeMax; it++){
2794 sum[it] += phvalue[ir][ic][it];
2800 Double_t sumcharge = 0.0;
2801 for(Int_t it = 0; it < fTimeMax; it++){
2802 sumcharge += sum[it];
2803 if(sum[it] > 20.0) nbcl++;
2807 /////////////////////////////////////////////////////////
2809 ////////////////////////////////////////////////////////
2810 if(fDebugLevel > 0){
2811 if ( !fDebugStreamer ) {
2813 TDirectory *backup = gDirectory;
2814 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2815 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2818 Double_t amph0 = sum[0];
2819 Double_t amphlast = sum[fTimeMax-1];
2820 Double_t rms = TMath::RMS(fTimeMax,sum);
2821 Int_t goodtracklet = (Int_t) fGoodTracklet;
2822 for(Int_t it = 0; it < fTimeMax; it++){
2823 Double_t clustera = sum[it];
2825 (* fDebugStreamer) << "FillDAQa"<<
2826 "ampTotal="<<sumcharge<<
2829 "detector="<<idect<<
2831 "amphlast="<<amphlast<<
2832 "goodtracklet="<<goodtracklet<<
2833 "clustera="<<clustera<<
2840 ////////////////////////////////////////////////////////
2842 ///////////////////////////////////////////////////////
2843 if(sum[0] > 100.0) used = 1;
2844 if(nbcl < fNumberClusters) used = 1;
2845 if(nbcl > fNumberClustersf) used = 1;
2847 //if(fDetectorPreviousTrack == 15){
2848 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2850 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2852 for(Int_t it = 0; it < fTimeMax; it++){
2853 UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2857 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2864 //____________Online trackling in AliTRDtrigger________________________________
2865 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2869 // Fill a simple average pulse height
2873 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2879 //____________Write_____________________________________________________
2880 //_____________________________________________________________________
2881 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2884 // Write infos to file
2888 if ( fDebugStreamer ) {
2889 delete fDebugStreamer;
2890 fDebugStreamer = 0x0;
2893 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2898 ,fNumberUsedPh[1]));
2900 TDirectory *backup = gDirectory;
2906 option = "recreate";
2908 TFile f(filename,option.Data());
2910 TStopwatch stopwatch;
2913 f.WriteTObject(fCalibraVector);
2918 f.WriteTObject(fCH2d);
2923 f.WriteTObject(fPH2d);
2928 f.WriteTObject(fPRF2d);
2931 if(fLinearFitterOn){
2932 AnalyseLinearFitter();
2933 f.WriteTObject(fLinearVdriftFit);
2938 if ( backup ) backup->cd();
2940 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2941 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2943 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2945 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2946 //___________________________________________probe the histos__________________________________________________
2947 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2950 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2951 // debug mode with 2 for TH2I and 3 for TProfile2D
2952 // It gives a pointer to a Double_t[7] with the info following...
2953 // [0] : number of calibration groups with entries
2954 // [1] : minimal number of entries found
2955 // [2] : calibration group number of the min
2956 // [3] : maximal number of entries found
2957 // [4] : calibration group number of the max
2958 // [5] : mean number of entries found
2959 // [6] : mean relative error
2962 Double_t *info = new Double_t[7];
2964 // Number of Xbins (detectors or groups of pads)
2965 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2966 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2969 Double_t nbwe = 0; //number of calibration groups with entries
2970 Double_t minentries = 0; //minimal number of entries found
2971 Double_t maxentries = 0; //maximal number of entries found
2972 Double_t placemin = 0; //calibration group number of the min
2973 Double_t placemax = -1; //calibration group number of the max
2974 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2975 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2977 Double_t counter = 0;
2980 TH1F *nbEntries = 0x0;//distribution of the number of entries
2981 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2982 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2984 // Beginning of the loop over the calibration groups
2985 for (Int_t idect = 0; idect < nbins; idect++) {
2987 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2988 projch->SetDirectory(0);
2990 // Number of entries for this calibration group
2991 Double_t nentries = 0.0;
2993 for (Int_t k = 0; k < nxbins; k++) {
2994 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2998 for (Int_t k = 0; k < nxbins; k++) {
2999 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3000 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3001 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3009 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3010 nbEntries = new TH1F("Number of entries","Number of entries"
3011 ,100,(Int_t)nentries/2,nentries*2);
3012 nbEntries->SetDirectory(0);
3013 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3015 nbEntriesPerGroup->SetDirectory(0);
3016 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3017 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3018 nbEntriesPerSp->SetDirectory(0);
3021 if(nentries > 0) nbEntries->Fill(nentries);
3022 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3023 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3028 if(nentries > maxentries){
3029 maxentries = nentries;
3033 minentries = nentries;
3035 if(nentries < minentries){
3036 minentries = nentries;
3042 meanstats += nentries;
3044 }//calibration groups loop
3046 if(nbwe > 0) meanstats /= nbwe;
3047 if(counter > 0) meanrelativerror /= counter;
3049 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3050 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3051 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3052 AliInfo(Form("The mean number of entries is %f",meanstats));
3053 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3056 info[1] = minentries;
3058 info[3] = maxentries;
3060 info[5] = meanstats;
3061 info[6] = meanrelativerror;
3064 gStyle->SetPalette(1);
3065 gStyle->SetOptStat(1111);
3066 gStyle->SetPadBorderMode(0);
3067 gStyle->SetCanvasColor(10);
3068 gStyle->SetPadLeftMargin(0.13);
3069 gStyle->SetPadRightMargin(0.01);
3070 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3073 nbEntries->Draw("");
3075 nbEntriesPerSp->SetStats(0);
3076 nbEntriesPerSp->Draw("");
3077 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3079 nbEntriesPerGroup->SetStats(0);
3080 nbEntriesPerGroup->Draw("");
3086 //____________________________________________________________________________
3087 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3090 // Return a Int_t[4] with:
3091 // 0 Mean number of entries
3092 // 1 median of number of entries
3093 // 2 rms of number of entries
3094 // 3 number of group with entries
3097 Double_t *stat = new Double_t[4];
3100 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3101 Double_t *weight = new Double_t[nbofgroups];
3102 Int_t *nonul = new Int_t[nbofgroups];
3104 for(Int_t k = 0; k < nbofgroups; k++){
3105 if(fEntriesCH[k] > 0) {
3107 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3110 else weight[k] = 0.0;
3112 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3113 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3114 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3119 //____________________________________________________________________________
3120 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3123 // Return a Int_t[4] with:
3124 // 0 Mean number of entries
3125 // 1 median of number of entries
3126 // 2 rms of number of entries
3127 // 3 number of group with entries
3130 Double_t *stat = new Double_t[4];
3133 Int_t nbofgroups = 540;
3134 Double_t *weight = new Double_t[nbofgroups];
3135 Int_t *nonul = new Int_t[nbofgroups];
3137 for(Int_t k = 0; k < nbofgroups; k++){
3138 if(fEntriesLinearFitter[k] > 0) {
3140 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3143 else weight[k] = 0.0;
3145 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3146 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3147 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3152 //////////////////////////////////////////////////////////////////////////////////////
3154 //////////////////////////////////////////////////////////////////////////////////////
3155 //_____________________________________________________________________________
3156 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3159 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3160 // If fNgroupprf is zero then no binning in tnp
3164 name += fCalibraMode->GetNz(2);
3166 name += fCalibraMode->GetNrphi(2);
3170 if(fNgroupprf != 0){
3172 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3173 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3174 fPRF2d->SetYTitle("Det/pad groups");
3175 fPRF2d->SetXTitle("Position x/W [pad width units]");
3176 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3177 fPRF2d->SetStats(0);
3180 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3181 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3182 fPRF2d->SetYTitle("Det/pad groups");
3183 fPRF2d->SetXTitle("Position x/W [pad width units]");
3184 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3185 fPRF2d->SetStats(0);
3190 //_____________________________________________________________________________
3191 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3194 // Create the 2D histos
3198 name += fCalibraMode->GetNz(1);
3200 name += fCalibraMode->GetNrphi(1);
3202 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3203 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3205 fPH2d->SetYTitle("Det/pad groups");
3206 fPH2d->SetXTitle("time [#mus]");
3207 fPH2d->SetZTitle("<PH> [a.u.]");
3211 //_____________________________________________________________________________
3212 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3215 // Create the 2D histos
3219 name += fCalibraMode->GetNz(0);
3221 name += fCalibraMode->GetNrphi(0);
3223 fCH2d = new TH2I("CH2d",(const Char_t *) name
3224 ,fNumberBinCharge,0,300,nn,0,nn);
3225 fCH2d->SetYTitle("Det/pad groups");
3226 fCH2d->SetXTitle("charge deposit [a.u]");
3227 fCH2d->SetZTitle("counts");
3232 //////////////////////////////////////////////////////////////////////////////////
3233 // Set relative scale
3234 /////////////////////////////////////////////////////////////////////////////////
3235 //_____________________________________________________________________________
3236 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3239 // Set the factor that will divide the deposited charge
3240 // to fit in the histo range [0,300]
3243 if (RelativeScale > 0.0) {
3244 fRelativeScale = RelativeScale;
3247 AliInfo("RelativeScale must be strict positif!");
3251 //////////////////////////////////////////////////////////////////////////////////
3252 // Quick way to fill a histo
3253 //////////////////////////////////////////////////////////////////////////////////
3254 //_____________________________________________________________________
3255 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3258 // FillCH2d: Marian style
3261 //skip simply the value out of range
3262 if((y>=300.0) || (y<0.0)) return;
3264 //Calcul the y place
3265 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3266 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3269 fCH2d->GetArray()[place]++;
3273 //////////////////////////////////////////////////////////////////////////////////
3274 // Geometrical functions
3275 ///////////////////////////////////////////////////////////////////////////////////
3276 //_____________________________________________________________________________
3277 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3280 // Reconstruct the plane number from the detector number
3283 return ((Int_t) (d % 6));
3287 //_____________________________________________________________________________
3288 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3291 // Reconstruct the chamber number from the detector number
3295 return ((Int_t) (d % 30) / fgkNplan);
3299 //_____________________________________________________________________________
3300 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3303 // Reconstruct the sector number from the detector number
3307 return ((Int_t) (d / fg));
3310 ///////////////////////////////////////////////////////////////////////////////////
3311 // Getter functions for DAQ of the CH2d and the PH2d
3312 //////////////////////////////////////////////////////////////////////////////////
3313 //_____________________________________________________________________
3314 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3317 // return pointer to fPH2d TProfile2D
3318 // create a new TProfile2D if it doesn't exist allready
3324 fTimeMax = nbtimebin;
3325 fSf = samplefrequency;
3331 //_____________________________________________________________________
3332 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3335 // return pointer to fCH2d TH2I
3336 // create a new TH2I if it doesn't exist allready
3345 ////////////////////////////////////////////////////////////////////////////////////////////
3346 // Drift velocity calibration
3347 ///////////////////////////////////////////////////////////////////////////////////////////
3348 //_____________________________________________________________________
3349 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3352 // return pointer to TLinearFitter Calibration
3353 // if force is true create a new TLinearFitter if it doesn't exist allready
3356 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3357 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3360 // if we are forced and TLinearFitter doesn't yet exist create it
3362 // new TLinearFitter
3363 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3364 fLinearFitterArray.AddAt(linearfitter,detector);
3365 return linearfitter;
3368 //____________________________________________________________________________
3369 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3372 // Analyse array of linear fitter because can not be written
3373 // Store two arrays: one with the param the other one with the error param + number of entries
3376 for(Int_t k = 0; k < 540; k++){
3377 TLinearFitter *linearfitter = GetLinearFitter(k);
3378 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3379 TVectorD *par = new TVectorD(2);
3380 TVectorD pare = TVectorD(2);
3381 TVectorD *parE = new TVectorD(3);
3382 linearfitter->Eval();
3383 linearfitter->GetParameters(*par);
3384 linearfitter->GetErrors(pare);
3385 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3386 (*parE)[0] = pare[0]*ppointError;
3387 (*parE)[1] = pare[1]*ppointError;
3388 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3389 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3390 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);