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();
346 fNumberClustersf = fTimeMax;
347 fNumberClusters = (Int_t)(0.6*fTimeMax);
349 //calib object from database used for reconstruction
350 if(fCalDetGain) delete fCalDetGain;
351 fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
353 // Calcul Xbins Chambd0, Chamb2
354 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
355 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
356 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
358 // If vector method On initialised all the stuff
360 fCalibraVector = new AliTRDCalibraVector();
361 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
362 fCalibraVector->SetTimeMax(fTimeMax);
363 if(fNgroupprf != 0) {
364 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
365 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
368 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
369 fCalibraVector->SetPRFRange(1.5);
371 for(Int_t k = 0; k < 3; k++){
372 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
373 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
375 TString namech("Nz");
376 namech += fCalibraMode->GetNz(0);
378 namech += fCalibraMode->GetNrphi(0);
379 fCalibraVector->SetNameCH((const char* ) namech);
380 TString nameph("Nz");
381 nameph += fCalibraMode->GetNz(1);
383 nameph += fCalibraMode->GetNrphi(1);
384 fCalibraVector->SetNamePH((const char* ) nameph);
385 TString nameprf("Nz");
386 nameprf += fCalibraMode->GetNz(2);
388 nameprf += fCalibraMode->GetNrphi(2);
390 nameprf += fNgroupprf;
391 fCalibraVector->SetNamePRF((const char* ) nameprf);
394 // Create the 2D histos corresponding to the pad groupCalibration mode
397 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
398 ,fCalibraMode->GetNz(0)
399 ,fCalibraMode->GetNrphi(0)));
401 // Create the 2D histo
406 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
407 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
411 fEntriesCH = new Int_t[ntotal0];
412 for(Int_t k = 0; k < ntotal0; k++){
419 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
420 ,fCalibraMode->GetNz(1)
421 ,fCalibraMode->GetNrphi(1)));
423 // Create the 2D histo
428 fPHPlace = new Short_t[fTimeMax];
429 for (Int_t k = 0; k < fTimeMax; k++) {
432 fPHValue = new Float_t[fTimeMax];
433 for (Int_t k = 0; k < fTimeMax; k++) {
437 if (fLinearFitterOn) {
438 //fLinearFitterArray.Expand(540);
439 fLinearFitterArray.SetName("ArrayLinearFitters");
440 fEntriesLinearFitter = new Int_t[540];
441 for(Int_t k = 0; k < 540; k++){
442 fEntriesLinearFitter[k] = 0;
444 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
449 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
450 ,fCalibraMode->GetNz(2)
451 ,fCalibraMode->GetNrphi(2)));
452 // Create the 2D histo
454 CreatePRF2d(ntotal2);
461 //____________Offline tracking in the AliTRDtracker____________________________
462 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
465 // Use AliTRDtrack for the calibration
469 AliTRDcluster *cl = 0x0; // pointeur to cluster
470 Int_t index0 = 0; // index of the first cluster in the current chamber
471 Int_t index1 = 0; // index of the current cluster in the current chamber
473 //////////////////////////////
474 // loop over the clusters
475 ///////////////////////////////
476 while((cl = t->GetCluster(index1))){
478 /////////////////////////////////////////////////////////////////////////
479 // Fill the infos for the previous clusters if not the same detector
480 ////////////////////////////////////////////////////////////////////////
481 Int_t detector = cl->GetDetector();
482 if ((detector != fDetectorPreviousTrack) &&
483 (index0 != index1)) {
487 //If the same track, then look if the previous detector is in
488 //the same plane, if yes: not a good track
489 if ((GetPlane(detector) == GetPlane(fDetectorPreviousTrack))) {
493 // Fill only if the track doesn't touch a masked pad or doesn't
496 // drift velocity unables to cut bad tracklets
497 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
501 FillTheInfoOfTheTrackCH(index1-index0);
506 FillTheInfoOfTheTrackPH();
509 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
512 } // if a good tracklet
515 ResetfVariablestracklet();
518 } // Fill at the end the charge
521 /////////////////////////////////////////////////////////////
522 // Localise and take the calib gain object for the detector
523 ////////////////////////////////////////////////////////////
524 if (detector != fDetectorPreviousTrack) {
526 //Localise the detector
527 LocalisationDetectorXbins(detector);
530 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
532 AliInfo("Could not get calibDB");
537 if( fCalROCGain ) delete fCalROCGain;
538 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
542 // Reset the detectbjobsor
543 fDetectorPreviousTrack = detector;
545 ///////////////////////////////////////
546 // Store the info of the cluster
547 ///////////////////////////////////////
548 Int_t row = cl->GetPadRow();
549 Int_t col = cl->GetPadCol();
550 CheckGoodTrackletV1(cl);
551 Int_t group[2] = {0,0};
552 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
553 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
554 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
558 } // while on clusters
560 ///////////////////////////
561 // Fill the last plane
562 //////////////////////////
563 if( index0 != index1 ){
569 // drift velocity unables to cut bad tracklets
570 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
574 FillTheInfoOfTheTrackCH(index1-index0);
579 FillTheInfoOfTheTrackPH();
582 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
584 } // if a good tracklet
589 ResetfVariablestracklet();
594 //____________Offline tracking in the AliTRDtracker____________________________
595 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
598 // Use AliTRDtrackV1 for the calibration
602 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
603 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
606 ///////////////////////////
607 // loop over the tracklet
608 ///////////////////////////
609 for(Int_t itr = 0; itr < 6; itr++){
611 if(!(tracklet = t->GetTracklet(itr))) continue;
612 if(!tracklet->IsOK()) continue;
614 ResetfVariablestracklet();
616 //////////////////////////////////////////
617 // localisation of the tracklet and dqdl
618 //////////////////////////////////////////
619 Int_t plane = tracklet->GetPlane();
621 while(!(cl = tracklet->GetClusters(ic++))) continue;
622 Int_t detector = cl->GetDetector();
623 if (detector != fDetectorPreviousTrack) {
624 // don't use the rest of this track if in the same plane
625 if ((plane == GetPlane(fDetectorPreviousTrack))) {
628 //Localise the detector bin
629 LocalisationDetectorXbins(detector);
631 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
633 AliInfo("Could not get calibDB");
637 if( fCalROCGain ) delete fCalROCGain;
638 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
640 fDetectorPreviousTrack = detector;
644 ////////////////////////////
645 // loop over the clusters
646 ////////////////////////////
647 Int_t nbclusters = 0;
648 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
649 if(!(cl = tracklet->GetClusters(ic))) continue;
652 // Store the info bis of the tracklet
653 Int_t row = cl->GetPadRow();
654 Int_t col = cl->GetPadCol();
655 CheckGoodTrackletV1(cl);
656 Int_t group[2] = {0,0};
657 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
658 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
659 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(ic),group,row,col);
662 ////////////////////////////////////////
663 // Fill the stuffs if a good tracklet
664 ////////////////////////////////////////
667 // drift velocity unables to cut bad tracklets
668 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
672 FillTheInfoOfTheTrackCH(nbclusters);
677 FillTheInfoOfTheTrackPH();
680 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
682 } // if a good tracklet
690 ///////////////////////////////////////////////////////////////////////////////////
691 // Routine inside the update with AliTRDtrack
692 ///////////////////////////////////////////////////////////////////////////////////
693 //____________Offine tracking in the AliTRDtracker_____________________________
694 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
697 // Drift velocity calibration:
698 // Fit the clusters with a straight line
699 // From the slope find the drift velocity
702 //Number of points: if less than 3 return kFALSE
703 Int_t npoints = index1-index0;
704 if(npoints <= 2) return kFALSE;
709 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
710 // parameters of the track
711 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in the plan yx track
712 Double_t tgl = t->GetTglPlane(GetPlane(detector)); // dz/dl and not dz/dx!
713 Double_t tnp = 0.0; // tan angle in the plan xy track
714 if( TMath::Abs(snp) < 1.){
715 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
717 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
718 // tilting pad and cross row
719 Int_t crossrow = 0; // if it crosses a pad row
720 Int_t rowp = -1; // if it crosses a pad row
721 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
722 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
723 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
725 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
726 Double_t dydt = 0.0; // dydt tracklet after straight line fit
727 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
728 Double_t pointError = 0.0; // error after straight line fit
729 Int_t nbli = 0; // number linear fitter points
730 linearFitterTracklet.StoreData(kFALSE);
731 linearFitterTracklet.ClearPoints();
733 //////////////////////////////
734 // loop over clusters
735 ////////////////////////////
736 for(Int_t k = 0; k < npoints; k++){
738 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
739 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
740 Double_t ycluster = cl->GetY();
741 Int_t time = cl->GetPadTime();
742 Double_t timeis = time/fSf;
743 //Double_t q = cl->GetQ();
744 //See if cross two pad rows
745 Int_t row = cl->GetPadRow();
747 if(row != rowp) crossrow = 1;
749 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
754 //////////////////////////////
756 //////////////////////////////
757 if(nbli <= 2) return kFALSE;
759 linearFitterTracklet.Eval();
760 linearFitterTracklet.GetParameters(pars);
761 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
762 errorpar = linearFitterTracklet.GetParError(1)*pointError;
765 /////////////////////////////
767 ////////////////////////////
769 if ( !fDebugStreamer ) {
771 TDirectory *backup = gDirectory;
772 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
773 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
777 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
778 //"snpright="<<snpright<<
779 "npoints="<<npoints<<
781 "detector="<<detector<<
788 "crossrow="<<crossrow<<
789 "errorpar="<<errorpar<<
790 "pointError="<<pointError<<
794 Int_t nbclusters = index1-index0;
795 Int_t plane = GetPlane(fDetectorPreviousTrack);
797 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
798 //"snpright="<<snpright<<
799 "nbclusters="<<nbclusters<<
800 "detector="<<fDetectorPreviousTrack<<
806 ///////////////////////////
808 ///////////////////////////
809 if(npoints < fNumberClusters) return kFALSE;
810 if(npoints > fNumberClustersf) return kFALSE;
811 if(pointError >= 0.1) return kFALSE;
812 if(crossrow == 1) return kFALSE;
814 ////////////////////////////
816 ////////////////////////////
818 //Add to the linear fitter of the detector
819 if( TMath::Abs(snp) < 1.){
820 Double_t x = tnp-dzdx*tnt;
821 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
822 if(fLinearFitterDebugOn) {
823 fLinearVdriftFit->Update(detector,x,pars[1]);
825 fEntriesLinearFitter[detector]++;
831 //____________Offine tracking in the AliTRDtracker_____________________________
832 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
835 // Drift velocity calibration:
836 // Fit the clusters with a straight line
837 // From the slope find the drift velocity
840 ////////////////////////////////////////////////
841 //Number of points: if less than 3 return kFALSE
842 /////////////////////////////////////////////////
843 if(nbclusters <= 2) return kFALSE;
848 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
849 // results of the linear fit
850 Double_t dydt = 0.0; // dydt tracklet after straight line fit
851 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
852 Double_t pointError = 0.0; // error after straight line fit
853 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
854 Int_t crossrow = 0; // if it crosses a pad row
855 Int_t rowp = -1; // if it crosses a pad row
856 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
858 linearFitterTracklet.StoreData(kFALSE);
859 linearFitterTracklet.ClearPoints();
861 ///////////////////////////////////////////
862 // Take the parameters of the track
863 //////////////////////////////////////////
864 // take now the snp, tnp and tgl from the track
865 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
866 Double_t tnp = 0.0; // dy/dx at the end of the chamber
867 if( TMath::Abs(snp) < 1.){
868 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
870 Double_t tgl = tracklet->GetTgl(); // dz/dl
871 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
873 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
874 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
875 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
876 // at the end with correction due to linear fit
877 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
878 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
881 ////////////////////////////
882 // loop over the clusters
883 ////////////////////////////
885 AliTRDcluster *cl = 0x0;
886 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
887 if(!(cl = tracklet->GetClusters(ic))) continue;
888 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
890 Double_t ycluster = cl->GetY();
891 Int_t time = cl->GetPadTime();
892 Double_t timeis = time/fSf;
893 //See if cross two pad rows
894 Int_t row = cl->GetPadRow();
895 if(rowp==-1) rowp = row;
896 if(row != rowp) crossrow = 1;
898 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
903 ////////////////////////////////////
904 // Do the straight line fit now
905 ///////////////////////////////////
907 linearFitterTracklet.Eval();
908 linearFitterTracklet.GetParameters(pars);
909 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
910 errorpar = linearFitterTracklet.GetParError(1)*pointError;
913 ////////////////////////////////
915 ///////////////////////////////
919 if ( !fDebugStreamer ) {
921 TDirectory *backup = gDirectory;
922 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
923 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
927 Int_t plane = GetPlane(fDetectorPreviousTrack);
929 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
930 //"snpright="<<snpright<<
931 "nbclusters="<<nbclusters<<
932 "detector="<<fDetectorPreviousTrack<<
940 "crossrow="<<crossrow<<
941 "errorpar="<<errorpar<<
942 "pointError="<<pointError<<
947 /////////////////////////
949 ////////////////////////
951 if(nbclusters < fNumberClusters) return kFALSE;
952 if(nbclusters > fNumberClustersf) return kFALSE;
953 if(pointError >= 0.1) return kFALSE;
954 if(crossrow == 1) return kFALSE;
956 ///////////////////////
958 //////////////////////
961 //Add to the linear fitter of the detector
962 if( TMath::Abs(snp) < 1.){
963 Double_t x = tnp-dzdx*tnt;
964 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
965 if(fLinearFitterDebugOn) {
966 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
968 fEntriesLinearFitter[fDetectorPreviousTrack]++;
974 //____________Offine tracking in the AliTRDtracker_____________________________
975 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
978 // PRF width calibration
979 // Assume a Gaussian shape: determinate the position of the three pad clusters
980 // Fit with a straight line
981 // Take the fitted values for all the clusters (3 or 2 pad clusters)
982 // Fill the PRF as function of angle of the track
987 //////////////////////////
989 /////////////////////////
990 Int_t npoints = index1-index0; // number of total points
991 Int_t nb3pc = 0; // number of three pads clusters used for fit
992 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
993 // To see the difference due to the fit
994 Double_t *padPositions;
995 padPositions = new Double_t[npoints];
996 for(Int_t k = 0; k < npoints; k++){
997 padPositions[k] = 0.0;
999 // Take the tgl and snp with the track t now
1000 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
1001 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
1002 Float_t dzdx = 0.0; // dzdx
1004 if(TMath::Abs(snp) < 1.0){
1005 tnp = snp / (TMath::Sqrt(1-snp*snp));
1006 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1009 TLinearFitter fitter(2,"pol1");
1010 fitter.StoreData(kFALSE);
1011 fitter.ClearPoints();
1014 ///////////////////////////
1015 // calcul the tnp group
1016 ///////////////////////////
1017 Bool_t echec = kFALSE;
1018 Double_t shift = 0.0;
1019 //Calculate the shift in x coresponding to this tnp
1020 if(fNgroupprf != 0.0){
1021 shift = -3.0*(fNgroupprf-1)-1.5;
1022 Double_t limithigh = -0.2*(fNgroupprf-1);
1023 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1025 while(tnp > limithigh){
1031 if(echec) return kFALSE;
1034 //////////////////////
1036 /////////////////////
1037 for(Int_t k = 0; k < npoints; k++){
1039 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1040 Short_t *signals = cl->GetSignals();
1041 Double_t time = cl->GetPadTime();
1042 //Calculate x if possible
1043 Float_t xcenter = 0.0;
1044 Bool_t echec = kTRUE;
1045 if((time<=7) || (time>=21)) continue;
1046 // Center 3 balanced: position with the center of the pad
1047 if ((((Float_t) signals[3]) > 0.0) &&
1048 (((Float_t) signals[2]) > 0.0) &&
1049 (((Float_t) signals[4]) > 0.0)) {
1051 // Security if the denomiateur is 0
1052 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1053 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1054 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1055 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1056 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1062 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1064 //if no echec: calculate with the position of the pad
1065 // Position of the cluster
1066 Double_t padPosition = xcenter + cl->GetPadCol();
1067 padPositions[k] = padPosition;
1069 fitter.AddPoint(&time, padPosition,1);
1073 /////////////////////////////
1075 ////////////////////////////
1076 if(nb3pc < 3) return kFALSE;
1079 fitter.GetParameters(line);
1080 Float_t pointError = -1.0;
1081 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1083 /////////////////////////////////////////////////////
1084 // Now fill the PRF: second loop over clusters
1085 /////////////////////////////////////////////////////
1086 for(Int_t k = 0; k < npoints; k++){
1088 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1089 Short_t *signals = cl->GetSignals(); // signal
1090 Double_t time = cl->GetPadTime(); // time bin
1091 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1092 Float_t padPos = cl->GetPadCol(); // middle pad
1093 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1094 Float_t ycenter = 0.0; // relative center charge
1095 Float_t ymin = 0.0; // relative left charge
1096 Float_t ymax = 0.0; // relative right charge
1098 //Requiere simply two pads clusters at least
1099 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1100 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1101 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1102 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1103 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1104 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1108 Int_t rowcl = cl->GetPadRow(); // row of cluster
1109 Int_t colcl = cl->GetPadCol(); // col of cluster
1110 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1111 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1112 Float_t xcl = cl->GetY(); // y cluster
1113 Float_t qcl = cl->GetQ(); // charge cluster
1114 Int_t plane = GetPlane(detector); // plane
1115 Int_t chamber = GetChamber(detector); // chamber
1116 Double_t xdiff = dpad; // reconstructed position constant
1117 Double_t x = dpad; // reconstructed position moved
1118 Float_t ep = pointError; // error of fit
1119 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1120 Float_t signal3 = (Float_t)signals[3]; // signal
1121 Float_t signal2 = (Float_t)signals[2]; // signal
1122 Float_t signal4 = (Float_t)signals[4]; // signal
1123 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1125 //////////////////////////////
1127 /////////////////////////////
1128 if(fDebugLevel > 0){
1129 if ( !fDebugStreamer ) {
1131 TDirectory *backup = gDirectory;
1132 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1133 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1139 Float_t y = ycenter;
1140 (* fDebugStreamer) << "HandlePRFtracklet"<<
1141 "caligroup="<<caligroup<<
1142 "detector="<<detector<<
1144 "chamber="<<chamber<<
1145 "npoints="<<npoints<<
1154 "padPosition="<<padPositions[k]<<
1155 "padPosTracklet="<<padPosTracklet<<
1160 "signal1="<<signal1<<
1161 "signal2="<<signal2<<
1162 "signal3="<<signal3<<
1163 "signal4="<<signal4<<
1164 "signal5="<<signal5<<
1170 (* fDebugStreamer) << "HandlePRFtracklet"<<
1171 "caligroup="<<caligroup<<
1172 "detector="<<detector<<
1174 "chamber="<<chamber<<
1175 "npoints="<<npoints<<
1184 "padPosition="<<padPositions[k]<<
1185 "padPosTracklet="<<padPosTracklet<<
1190 "signal1="<<signal1<<
1191 "signal2="<<signal2<<
1192 "signal3="<<signal3<<
1193 "signal4="<<signal4<<
1194 "signal5="<<signal5<<
1200 (* fDebugStreamer) << "HandlePRFtracklet"<<
1201 "caligroup="<<caligroup<<
1202 "detector="<<detector<<
1204 "chamber="<<chamber<<
1205 "npoints="<<npoints<<
1214 "padPosition="<<padPositions[k]<<
1215 "padPosTracklet="<<padPosTracklet<<
1220 "signal1="<<signal1<<
1221 "signal2="<<signal2<<
1222 "signal3="<<signal3<<
1223 "signal4="<<signal4<<
1224 "signal5="<<signal5<<
1230 ////////////////////////////
1232 ///////////////////////////
1233 if(npoints < fNumberClusters) continue;
1234 if(npoints > fNumberClustersf) continue;
1235 if(nb3pc <= 5) continue;
1236 if((time >= 21) || (time < 7)) continue;
1237 if(TMath::Abs(snp) >= 1.0) continue;
1238 if(TMath::Abs(qcl) < 80) continue;
1240 ////////////////////////////
1242 ///////////////////////////
1244 if(TMath::Abs(dpad) < 1.5) {
1245 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1246 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1248 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1249 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1250 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1252 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1253 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1254 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1258 if(TMath::Abs(dpad) < 1.5) {
1259 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1260 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1262 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1263 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1264 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1266 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1267 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1268 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1275 //____________Offine tracking in the AliTRDtracker_____________________________
1276 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1279 // PRF width calibration
1280 // Assume a Gaussian shape: determinate the position of the three pad clusters
1281 // Fit with a straight line
1282 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1283 // Fill the PRF as function of angle of the track
1287 ///////////////////////////////////////////
1288 // Take the parameters of the track
1289 //////////////////////////////////////////
1290 // take now the snp, tnp and tgl from the track
1291 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1292 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1293 if( TMath::Abs(snp) < 1.){
1294 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1296 Double_t tgl = tracklet->GetTgl(); // dz/dl
1297 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1299 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1300 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1301 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1302 // at the end with correction due to linear fit
1303 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1304 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1306 ///////////////////////////////
1307 // Calculate tnp group shift
1308 ///////////////////////////////
1309 Bool_t echec = kFALSE;
1310 Double_t shift = 0.0;
1311 //Calculate the shift in x coresponding to this tnp
1312 if(fNgroupprf != 0.0){
1313 shift = -3.0*(fNgroupprf-1)-1.5;
1314 Double_t limithigh = -0.2*(fNgroupprf-1);
1315 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1317 while(tnp > limithigh){
1323 // do nothing if out of tnp range
1324 if(echec) return kFALSE;
1326 ///////////////////////
1328 //////////////////////
1330 Int_t nb3pc = 0; // number of three pads clusters used for fit
1331 Double_t *padPositions; // to see the difference between the fit and the 3 pad clusters position
1332 padPositions = new Double_t[AliTRDseed::knTimebins];
1333 for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1334 padPositions[k] = 0.0;
1336 TLinearFitter fitter(2,"pol1"); // TLinearFitter for the linear fit in the drift region
1337 fitter.StoreData(kFALSE);
1338 fitter.ClearPoints();
1340 ////////////////////////////
1341 // loop over the clusters
1342 ////////////////////////////
1343 AliTRDcluster *cl = 0x0;
1344 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1345 if(!(cl = tracklet->GetClusters(ic))) continue;
1347 Double_t time = cl->GetPadTime();
1348 if((time<=7) || (time>=21)) continue;
1349 Short_t *signals = cl->GetSignals();
1350 Float_t xcenter = 0.0;
1351 Bool_t echec = kTRUE;
1353 /////////////////////////////////////////////////////////////
1354 // Center 3 balanced: position with the center of the pad
1355 /////////////////////////////////////////////////////////////
1356 if ((((Float_t) signals[3]) > 0.0) &&
1357 (((Float_t) signals[2]) > 0.0) &&
1358 (((Float_t) signals[4]) > 0.0)) {
1360 // Security if the denomiateur is 0
1361 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1362 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1363 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1364 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1365 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1371 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1374 ////////////////////////////////////////////////////////
1375 //if no echec: calculate with the position of the pad
1376 // Position of the cluster
1377 // fill the linear fitter
1378 ///////////////////////////////////////////////////////
1379 Double_t padPosition = xcenter + cl->GetPadCol();
1380 padPositions[ic] = padPosition;
1382 fitter.AddPoint(&time, padPosition,1);
1388 //////////////////////////////
1389 // fit with a straight line
1390 /////////////////////////////
1391 if(nb3pc < 3) return kFALSE;
1394 fitter.GetParameters(line);
1395 Float_t pointError = -1.0;
1396 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
1398 ////////////////////////////////////////////////
1399 // Fill the PRF: Second loop over clusters
1400 //////////////////////////////////////////////
1401 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1402 if(!(cl = tracklet->GetClusters(ic))) continue;
1404 Short_t *signals = cl->GetSignals(); // signal
1405 Double_t time = cl->GetPadTime(); // time bin
1406 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1407 Float_t padPos = cl->GetPadCol(); // middle pad
1408 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1409 Float_t ycenter = 0.0; // relative center charge
1410 Float_t ymin = 0.0; // relative left charge
1411 Float_t ymax = 0.0; // relative right charge
1413 ////////////////////////////////////////////////////////////////
1414 // Calculate ycenter, ymin and ymax even for two pad clusters
1415 ////////////////////////////////////////////////////////////////
1416 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1417 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1418 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1419 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1420 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1421 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1424 /////////////////////////
1425 // Calibration group
1426 ////////////////////////
1427 Int_t rowcl = cl->GetPadRow(); // row of cluster
1428 Int_t colcl = cl->GetPadCol(); // col of cluster
1429 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1430 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1431 Float_t xcl = cl->GetY(); // y cluster
1432 Float_t qcl = cl->GetQ(); // charge cluster
1433 Int_t plane = GetPlane(fDetectorPreviousTrack); // plane
1434 Int_t chamber = GetChamber(fDetectorPreviousTrack); // chamber
1435 Double_t xdiff = dpad; // reconstructed position constant
1436 Double_t x = dpad; // reconstructed position moved
1437 Float_t ep = pointError; // error of fit
1438 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1439 Float_t signal3 = (Float_t)signals[3]; // signal
1440 Float_t signal2 = (Float_t)signals[2]; // signal
1441 Float_t signal4 = (Float_t)signals[4]; // signal
1442 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1444 /////////////////////
1446 ////////////////////
1448 if(fDebugLevel > 0){
1449 if ( !fDebugStreamer ) {
1451 TDirectory *backup = gDirectory;
1452 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1453 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1458 Float_t y = ycenter;
1459 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1460 "caligroup="<<caligroup<<
1461 "detector="<<fDetectorPreviousTrack<<
1463 "chamber="<<chamber<<
1464 "npoints="<<nbclusters<<
1473 "padPosition="<<padPositions[ic]<<
1474 "padPosTracklet="<<padPosTracklet<<
1479 "signal1="<<signal1<<
1480 "signal2="<<signal2<<
1481 "signal3="<<signal3<<
1482 "signal4="<<signal4<<
1483 "signal5="<<signal5<<
1489 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1490 "caligroup="<<caligroup<<
1491 "detector="<<fDetectorPreviousTrack<<
1493 "chamber="<<chamber<<
1494 "npoints="<<nbclusters<<
1503 "padPosition="<<padPositions[ic]<<
1504 "padPosTracklet="<<padPosTracklet<<
1509 "signal1="<<signal1<<
1510 "signal2="<<signal2<<
1511 "signal3="<<signal3<<
1512 "signal4="<<signal4<<
1513 "signal5="<<signal5<<
1519 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1520 "caligroup="<<caligroup<<
1521 "detector="<<fDetectorPreviousTrack<<
1523 "chamber="<<chamber<<
1524 "npoints="<<nbclusters<<
1533 "padPosition="<<padPositions[ic]<<
1534 "padPosTracklet="<<padPosTracklet<<
1539 "signal1="<<signal1<<
1540 "signal2="<<signal2<<
1541 "signal3="<<signal3<<
1542 "signal4="<<signal4<<
1543 "signal5="<<signal5<<
1549 /////////////////////
1551 /////////////////////
1552 if(nbclusters < fNumberClusters) continue;
1553 if(nbclusters > fNumberClustersf) continue;
1554 if(nb3pc <= 5) continue;
1555 if((time >= 21) || (time < 7)) continue;
1556 if(TMath::Abs(qcl) < 80) continue;
1557 if( TMath::Abs(snp) > 1.) continue;
1560 ////////////////////////
1562 ///////////////////////
1564 if(TMath::Abs(dpad) < 1.5) {
1565 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1566 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1567 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1569 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1570 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1571 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1573 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1574 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1575 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1580 if(TMath::Abs(dpad) < 1.5) {
1581 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1582 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1584 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1585 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1586 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1588 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1589 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1590 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1593 } // second loop over clusters
1599 ///////////////////////////////////////////////////////////////////////////////////////
1600 // Pad row col stuff: see if masked or not
1601 ///////////////////////////////////////////////////////////////////////////////////////
1602 //_____________________________________________________________________________
1603 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1606 // See if we are not near a masked pad
1609 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1613 //_____________________________________________________________________________
1614 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1617 // See if we are not near a masked pad
1620 if (!IsPadOn(detector, col, row)) {
1621 fGoodTracklet = kFALSE;
1625 if (!IsPadOn(detector, col-1, row)) {
1626 fGoodTracklet = kFALSE;
1631 if (!IsPadOn(detector, col+1, row)) {
1632 fGoodTracklet = kFALSE;
1637 //_____________________________________________________________________________
1638 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1641 // Look in the choosen database if the pad is On.
1642 // If no the track will be "not good"
1645 // Get the parameter object
1646 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1648 AliInfo("Could not get calibDB");
1652 if (!cal->IsChamberInstalled(detector) ||
1653 cal->IsChamberMasked(detector) ||
1654 cal->IsPadMasked(detector,col,row)) {
1662 ///////////////////////////////////////////////////////////////////////////////////////
1663 // Calibration groups: calculate the number of groups, localise...
1664 ////////////////////////////////////////////////////////////////////////////////////////
1665 //_____________________________________________________________________________
1666 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1669 // Calculate the calibration group number for i
1672 // Row of the cluster and position in the pad groups
1674 if (fCalibraMode->GetNnZ(i) != 0) {
1675 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1679 // Col of the cluster and position in the pad groups
1681 if (fCalibraMode->GetNnRphi(i) != 0) {
1682 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1685 return posc*fCalibraMode->GetNfragZ(i)+posr;
1688 //____________________________________________________________________________________
1689 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1692 // Calculate the total number of calibration groups
1696 fCalibraMode->ModePadCalibration(2,i);
1697 fCalibraMode->ModePadFragmentation(0,2,0,i);
1698 fCalibraMode->SetDetChamb2(i);
1699 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1700 fCalibraMode->ModePadCalibration(0,i);
1701 fCalibraMode->ModePadFragmentation(0,0,0,i);
1702 fCalibraMode->SetDetChamb0(i);
1703 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1704 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1708 //_____________________________________________________________________________
1709 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1712 // Set the mode of calibration group in the z direction for the parameter i
1717 fCalibraMode->SetNz(i, Nz);
1720 AliInfo("You have to choose between 0 and 4");
1724 //_____________________________________________________________________________
1725 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1728 // Set the mode of calibration group in the rphi direction for the parameter i
1733 fCalibraMode->SetNrphi(i ,Nrphi);
1736 AliInfo("You have to choose between 0 and 6");
1740 //____________Set the pad calibration variables for the detector_______________
1741 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1744 // For the detector calcul the first Xbins and set the number of row
1745 // and col pads per calibration groups, the number of calibration
1746 // groups in the detector.
1749 // first Xbins of the detector
1751 fCalibraMode->CalculXBins(detector,0);
1754 fCalibraMode->CalculXBins(detector,1);
1757 fCalibraMode->CalculXBins(detector,2);
1760 // fragmentation of idect
1761 for (Int_t i = 0; i < 3; i++) {
1762 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1763 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1764 , (Int_t) GetChamber(detector)
1765 , (Int_t) GetSector(detector),i);
1771 //_____________________________________________________________________________
1772 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1775 // Should be between 0 and 6
1778 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1779 AliInfo("The number of groups must be between 0 and 6!");
1782 fNgroupprf = numberGroupsPRF;
1786 ///////////////////////////////////////////////////////////////////////////////////////////
1787 // Per tracklet: store or reset the info, fill the histos with the info
1788 //////////////////////////////////////////////////////////////////////////////////////////
1789 //_____________________________________________________________________________
1790 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1793 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1794 // Correct from the gain correction before
1797 // time bin of the cluster not corrected
1798 Int_t time = cl->GetPadTime();
1800 //Correct for the gain coefficient used in the database for reconstruction
1801 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1802 Float_t correction = 1.0;
1803 Float_t normalisation = 6.67;
1804 // we divide with gain in AliTRDclusterizer::Transform...
1805 if( correctthegain > 0 ) normalisation /= correctthegain;
1808 // take dd/dl corrected from the angle of the track
1809 correction = dqdl / (normalisation);
1812 // Fill the fAmpTotal with the charge
1814 if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1817 // Fill the fPHPlace and value
1819 fPHPlace[time] = group[1];
1820 fPHValue[time] = correction;
1824 //____________Offine tracking in the AliTRDtracker_____________________________
1825 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1828 // Reset values per tracklet
1831 //Reset good tracklet
1832 fGoodTracklet = kTRUE;
1834 // Reset the fPHValue
1836 //Reset the fPHValue and fPHPlace
1837 for (Int_t k = 0; k < fTimeMax; k++) {
1843 // Reset the fAmpTotal where we put value
1845 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1850 //____________Offine tracking in the AliTRDtracker_____________________________
1851 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1854 // For the offline tracking or mcm tracklets
1855 // This function will be called in the functions UpdateHistogram...
1856 // to fill the info of a track for the relativ gain calibration
1859 Int_t nb = 0; // Nombre de zones traversees
1860 Int_t fd = -1; // Premiere zone non nulle
1861 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1863 if(nbclusters < fNumberClusters) return;
1864 if(nbclusters > fNumberClustersf) return;
1867 // See if the track goes through different zones
1868 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1869 if (fAmpTotal[k] > 0.0) {
1870 totalcharge += fAmpTotal[k];
1883 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
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);
1893 if ((fAmpTotal[fd] > 0.0) &&
1894 (fAmpTotal[fd+1] > 0.0)) {
1895 // One of the two very big
1896 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1898 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1899 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1902 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1905 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1907 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1909 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1910 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1913 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1916 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1919 if (fCalibraMode->GetNfragZ(0) > 1) {
1920 if (fAmpTotal[fd] > 0.0) {
1921 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1922 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1923 // One of the two very big
1924 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1926 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1927 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1930 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1933 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1935 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1937 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1938 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1941 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1943 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1954 //____________Offine tracking in the AliTRDtracker_____________________________
1955 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1958 // For the offline tracking or mcm tracklets
1959 // This function will be called in the functions UpdateHistogram...
1960 // to fill the info of a track for the drift velocity calibration
1963 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1964 Int_t fd1 = -1; // Premiere zone non nulle
1965 Int_t fd2 = -1; // Deuxieme zone non nulle
1966 Int_t k1 = -1; // Debut de la premiere zone
1967 Int_t k2 = -1; // Debut de la seconde zone
1968 Int_t nbclusters = 0; // number of clusters
1972 // See if the track goes through different zones
1973 for (Int_t k = 0; k < fTimeMax; k++) {
1974 if (fPHValue[k] > 0.0) {
1980 if (fPHPlace[k] != fd1) {
1986 if (fPHPlace[k] != fd2) {
1993 if(nbclusters < fNumberClusters) return;
1994 if(nbclusters > fNumberClustersf) return;
2000 for (Int_t i = 0; i < fTimeMax; i++) {
2002 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2003 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2006 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2011 if ((fd1 == fd2+1) ||
2013 // One of the two fast all the think
2014 if (k2 > (k1+fDifference)) {
2015 //we choose to fill the fd1 with all the values
2017 for (Int_t i = 0; i < fTimeMax; i++) {
2019 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2022 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2026 if ((k2+fDifference) < fTimeMax) {
2027 //we choose to fill the fd2 with all the values
2029 for (Int_t i = 0; i < fTimeMax; i++) {
2031 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2034 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2039 // Two zones voisines sinon rien!
2040 if (fCalibraMode->GetNfragZ(1) > 1) {
2042 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2043 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2044 // One of the two fast all the think
2045 if (k2 > (k1+fDifference)) {
2046 //we choose to fill the fd1 with all the values
2048 for (Int_t i = 0; i < fTimeMax; i++) {
2050 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2053 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2057 if ((k2+fDifference) < fTimeMax) {
2058 //we choose to fill the fd2 with all the values
2060 for (Int_t i = 0; i < fTimeMax; i++) {
2062 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2065 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2071 // Two zones voisines sinon rien!
2073 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2074 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2075 // One of the two fast all the think
2076 if (k2 > (k1 + fDifference)) {
2077 //we choose to fill the fd1 with all the values
2079 for (Int_t i = 0; i < fTimeMax; i++) {
2081 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2084 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2088 if ((k2+fDifference) < fTimeMax) {
2089 //we choose to fill the fd2 with all the values
2091 for (Int_t i = 0; i < fTimeMax; i++) {
2093 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2096 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2107 //////////////////////////////////////////////////////////////////////////////////////////
2108 // DAQ process functions
2109 /////////////////////////////////////////////////////////////////////////////////////////
2110 //_____________________________________________________________________
2111 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2114 // Event Processing loop - AliTRDrawStreamBase
2115 // TestBeam 2007 version
2116 // 0 timebin problem
2120 // Same algorithm as TestBeam but different reader
2123 Int_t withInput = 1;
2125 Double_t phvalue[16][144][36];
2126 for(Int_t k = 0; k < 36; k++){
2127 for(Int_t j = 0; j < 16; j++){
2128 for(Int_t c = 0; c < 144; c++){
2129 phvalue[j][c][k] = 0.0;
2134 fDetectorPreviousTrack = -1;
2137 Int_t nbtimebin = 0;
2138 Int_t baseline = 10;
2145 while (rawStream->Next()) {
2147 Int_t idetector = rawStream->GetDet(); // current detector
2148 Int_t imcm = rawStream->GetMCM(); // current MCM
2149 Int_t irob = rawStream->GetROB(); // current ROB
2151 //printf("Detector %d\n",idetector);
2153 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2156 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2160 for(Int_t k = 0; k < 36; k++){
2161 for(Int_t j = 0; j < 16; j++){
2162 for(Int_t c = 0; c < 144; c++){
2163 phvalue[j][c][k] = 0.0;
2169 fDetectorPreviousTrack = idetector;
2170 fMCMPrevious = imcm;
2171 fROBPrevious = irob;
2173 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2174 if(nbtimebin == 0) return 0;
2175 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2176 fTimeMax = nbtimebin;
2178 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2179 fNumberClustersf = fTimeMax;
2180 fNumberClusters = (Int_t)(0.6*fTimeMax);
2183 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2184 Int_t col = rawStream->GetCol();
2185 Int_t row = rawStream->GetRow();
2188 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2191 for(Int_t itime = 0; itime < nbtimebin; itime++){
2192 phvalue[row][col][itime] = signal[itime]-baseline;
2196 // fill the last one
2197 if(fDetectorPreviousTrack != -1){
2200 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2203 for(Int_t k = 0; k < 36; k++){
2204 for(Int_t j = 0; j < 16; j++){
2205 for(Int_t c = 0; c < 144; c++){
2206 phvalue[j][c][k] = 0.0;
2215 while (rawStream->Next()) {
2217 Int_t idetector = rawStream->GetDet(); // current detector
2218 Int_t imcm = rawStream->GetMCM(); // current MCM
2219 Int_t irob = rawStream->GetROB(); // current ROB
2221 //printf("Detector %d\n",idetector);
2223 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2226 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2229 for(Int_t k = 0; k < 36; k++){
2230 for(Int_t j = 0; j < 16; j++){
2231 for(Int_t c = 0; c < 144; c++){
2232 phvalue[j][c][k] = 0.0;
2238 fDetectorPreviousTrack = idetector;
2239 fMCMPrevious = imcm;
2240 fROBPrevious = irob;
2242 //baseline = rawStream->GetCommonAdditive(); // common baseline
2244 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2245 fNumberClustersf = fTimeMax;
2246 fNumberClusters = (Int_t)(0.6*fTimeMax);
2247 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2248 Int_t col = rawStream->GetCol();
2249 Int_t row = rawStream->GetRow();
2252 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2254 for(Int_t itime = 0; itime < fTimeMax; itime++){
2255 phvalue[row][col][itime] = signal[itime]-baseline;
2259 // fill the last one
2260 if(fDetectorPreviousTrack != -1){
2263 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2266 for(Int_t k = 0; k < 36; k++){
2267 for(Int_t j = 0; j < 16; j++){
2268 for(Int_t c = 0; c < 144; c++){
2269 phvalue[j][c][k] = 0.0;
2279 //_____________________________________________________________________
2280 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2283 // Event Processing loop - AliTRDrawStreamBase
2284 // Use old AliTRDmcmtracklet code
2285 // 0 timebin problem
2289 // Algorithm with mcm tracklet
2292 Int_t withInput = 1;
2294 AliTRDmcm mcm = AliTRDmcm(0);
2295 AliTRDtrigParam *param = AliTRDtrigParam::Instance();
2296 rawStream->SetSharedPadReadout(kTRUE);
2298 fDetectorPreviousTrack = -1;
2302 Int_t nbtimebin = 0;
2303 Int_t baseline = 10;
2310 while (rawStream->Next()) {
2312 Int_t idetector = rawStream->GetDet(); // current detector
2313 Int_t imcm = rawStream->GetMCM(); // current MCM
2314 Int_t irob = rawStream->GetROB(); // current ROB
2315 row = rawStream->GetRow();
2318 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2321 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2327 mcm.SetChaId(idetector);
2329 mcm.SetColRange(0,21);
2332 if(fDetectorPreviousTrack == -1){
2335 mcm.SetChaId(idetector);
2337 mcm.SetColRange(0,21);
2341 fDetectorPreviousTrack = idetector;
2342 fMCMPrevious = imcm;
2343 fROBPrevious = irob;
2345 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2346 if(nbtimebin == 0) return 0;
2347 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2348 fTimeMax = nbtimebin;
2349 fNumberClustersf = fTimeMax;
2350 fNumberClusters = (Int_t)(0.6*fTimeMax);
2351 param->SetTimeRange(0,fTimeMax);
2353 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2355 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2356 Int_t adc = rawStream->GetADC();
2359 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2361 for(Int_t itime = 0; itime < nbtimebin; itime++){
2362 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2366 // fill the last one
2367 if(fDetectorPreviousTrack != -1){
2370 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2375 mcm.SetRobId(fROBPrevious);
2376 mcm.SetChaId(fDetectorPreviousTrack);
2378 mcm.SetColRange(0,21);
2385 while (rawStream->Next()) {
2387 Int_t idetector = rawStream->GetDet(); // current detector
2388 Int_t imcm = rawStream->GetMCM(); // current MCM
2389 Int_t irob = rawStream->GetROB(); // current ROB
2390 row = rawStream->GetRow();
2392 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2395 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2401 mcm.SetChaId(idetector);
2403 mcm.SetColRange(0,21);
2407 fDetectorPreviousTrack = idetector;
2408 fMCMPrevious = imcm;
2409 fROBPrevious = irob;
2411 //baseline = rawStream->GetCommonAdditive(); // common baseline
2413 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2414 fNumberClustersf = fTimeMax;
2415 fNumberClusters = (Int_t)(0.6*fTimeMax);
2416 param->SetTimeRange(0,fTimeMax);
2417 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2418 Int_t adc = rawStream->GetADC();
2421 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2423 for(Int_t itime = 0; itime < fTimeMax; itime++){
2424 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2428 // fill the last one
2429 if(fDetectorPreviousTrack != -1){
2432 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2436 mcm.SetRobId(fROBPrevious);
2437 mcm.SetChaId(fDetectorPreviousTrack);
2439 mcm.SetColRange(0,21);
2447 //_____________________________________________________________________
2448 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2451 // Event processing loop - AliRawReader
2452 // Testbeam 2007 version
2455 AliTRDrawStreamBase rawStream(rawReader);
2457 rawReader->Select("TRD");
2459 return ProcessEventDAQ(&rawStream, nocheck);
2462 //_________________________________________________________________________
2463 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2465 eventHeaderStruct *event,
2468 eventHeaderStruct* /*event*/,
2475 // process date event
2476 // Testbeam 2007 version
2479 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2480 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2484 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2489 //_____________________________________________________________________
2490 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliRawReader *rawReader, Bool_t nocheck)
2493 // Event processing loop - AliRawReader
2494 // use the old mcm traklet code
2497 AliTRDrawStreamBase rawStream(rawReader);
2499 rawReader->Select("TRD");
2501 return ProcessEventDAQV1(&rawStream, nocheck);
2503 //_________________________________________________________________________
2504 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(
2506 eventHeaderStruct *event,
2509 eventHeaderStruct* /*event*/,
2516 // process date event
2517 // use the old mcm tracklet code
2520 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2521 Int_t result=ProcessEventDAQV1(rawReader, nocheck);
2525 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2530 //////////////////////////////////////////////////////////////////////////////
2531 // Routine inside the DAQ process
2532 /////////////////////////////////////////////////////////////////////////////
2533 //_______________________________________________________________________
2534 Int_t AliTRDCalibraFillHisto::FillDAQ(AliTRDmcm *mcm){
2537 // Return 2 if some tracklets are found and used, 1 if nothing
2544 for (Int_t iSeed = 0; iSeed < 4; iSeed++) {
2546 if (mcm->GetSeedCol()[iSeed] < 0) {
2550 nbev += TestTracklet(mcm->GetChaId(),mcm->GetRow(),iSeed,mcm);
2555 if(nbev > 0) nbev = 2;
2561 //__________________________________________________________________________
2562 Int_t AliTRDCalibraFillHisto::TestTracklet( Int_t idet, Int_t row, Int_t iSeed, AliTRDmcm *mcm){
2565 // Build the tracklet and return if the tracklet if finally used or not (1/0)
2570 AliTRDmcmTracklet mcmtracklet = AliTRDmcmTracklet();
2571 //mcmtracklet.Reset();
2572 mcmtracklet.SetDetector(idet);
2573 mcmtracklet.SetRow(row);
2574 mcmtracklet.SetN(0);
2576 Int_t iCol, iCol1, iCol2, track[3];
2577 iCol = mcm->GetSeedCol()[iSeed]; // 0....20 (MCM)
2578 mcm->GetColRange(iCol1,iCol2); // range in the pad plane
2581 for (Int_t iTime = 0; iTime < fTimeMax; iTime++) {
2583 amp[0] = mcm->GetADC(iCol-1,iTime);
2584 amp[1] = mcm->GetADC(iCol ,iTime);
2585 amp[2] = mcm->GetADC(iCol+1,iTime);
2587 if(mcm->IsCluster(iCol,iTime)) {
2589 mcmtracklet.AddCluster(iCol+iCol1,iTime,amp,track);
2592 else if ((iCol+1+1) < 21) {
2594 amp[0] = mcm->GetADC(iCol-1+1,iTime);
2595 amp[1] = mcm->GetADC(iCol +1,iTime);
2596 amp[2] = mcm->GetADC(iCol+1+1,iTime);
2598 if(mcm->IsCluster(iCol+1,iTime)) {
2600 mcmtracklet.AddCluster(iCol+1+iCol1,iTime,amp,track);
2609 nbev = UpdateHistogramcm(&mcmtracklet);
2614 //____________Online trackling in AliTRDtrigger________________________________
2615 Int_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2618 // Return if the tracklet is finally used or not (1/0) for calibration
2623 //fGoodTracklet = kTRUE;
2625 // Localisation of the Xbins involved
2626 Int_t idect = trk->GetDetector();
2627 Int_t idectrue = trk->GetDetector();
2630 Int_t nbclusters = trk->GetNclusters();
2632 // Eventuelle correction due to track angle in z direction
2633 Float_t correction = 1.0;
2634 if (fMcmCorrectAngle) {
2635 Float_t z = trk->GetRowz();
2636 Float_t r = trk->GetTime0();
2637 correction = r / TMath::Sqrt((r*r+z*z));
2641 Int_t row = trk->GetRow();
2644 // Boucle sur les clusters
2645 // Condition on number of cluster: don't come from the middle of the detector
2648 for(Int_t k =0; k < 36; k++) amph[k]=0.0;
2649 Double_t ampTotal = 0.0;
2651 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2653 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2654 Int_t time = trk->GetClusterTime(icl);
2655 Int_t col = trk->GetClusterCol(icl);
2657 //CheckGoodTrackletV0(idect,row,col);
2659 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2660 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2661 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2663 ampTotal += (Float_t) (amp[0]+amp[1]+amp[2]);
2664 amph[time]=amp[0]+amp[1]+amp[2];
2666 if(fDebugLevel > 0){
2667 if ( !fDebugStreamer ) {
2669 TDirectory *backup = gDirectory;
2670 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2671 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2674 Double_t amp0 = amp[0];
2675 Double_t amp1 = amp[1];
2676 Double_t amp2 = amp[2];
2678 (* fDebugStreamer) << "UpdateHistogramcm0"<<
2679 "nbclusters="<<nbclusters<<
2686 "detector="<<idectrue<<
2690 } // Boucle clusters
2692 if((amph[0] > 100.0) || (!fGoodTracklet) || (trk->GetNclusters() < fNumberClusters) || (trk->GetNclusters() > fNumberClustersf)) used = 0;
2695 for(Int_t k = 0; k < fTimeMax; k++) UpdateDAQ(idect,0,0,k,amph[k],fTimeMax);
2696 //((TH2I *)GetCH2d()->Fill(ampTotal/30.0,idect));
2700 if(fDebugLevel > 0){
2701 if ( !fDebugStreamer ) {
2703 TDirectory *backup = gDirectory;
2704 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2705 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2708 Double_t amph0 = amph[0];
2709 Double_t amphlast = amph[fTimeMax-1];
2710 Double_t rms = TMath::RMS(fTimeMax,amph);
2711 Int_t goodtracklet = (Int_t) fGoodTracklet;
2713 (* fDebugStreamer) << "UpdateHistogramcm1"<<
2714 "nbclusters="<<nbclusters<<
2715 "ampTotal="<<ampTotal<<
2717 "detector="<<idectrue<<
2719 "amphlast="<<amphlast<<
2720 "goodtracklet="<<goodtracklet<<
2728 //_______________________________________________________________________
2729 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2732 // Look for the maximum by collapsing over the time
2733 // Sum over four pad col and two pad row
2739 Int_t idect = fDetectorPreviousTrack;
2740 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2742 for(Int_t tb = 0; tb < 36; tb++){
2746 //fGoodTracklet = kTRUE;
2747 //fDetectorPreviousTrack = 0;
2750 ///////////////////////////
2752 /////////////////////////
2756 Double_t integralMax = -1;
2758 for (Int_t ir = 1; ir <= 15; ir++)
2760 for (Int_t ic = 2; ic <= 142; ic++)
2762 Double_t integral = 0;
2763 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2765 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2767 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2768 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2771 for(Int_t tb = 0; tb< fTimeMax; tb++){
2772 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2778 if (integralMax < integral)
2782 integralMax = integral;
2783 } // check max integral
2787 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2789 if((imaxRow == 0) || (imaxCol == 0)) {
2793 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2794 //if(!fGoodTracklet) used = 1;;
2796 // /////////////////////////////////////////////////////
2797 // sum ober 2 row and 4 pad cols for each time bins
2798 // ////////////////////////////////////////////////////
2801 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2803 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2805 for(Int_t it = 0; it < fTimeMax; it++){
2806 sum[it] += phvalue[ir][ic][it];
2812 Double_t sumcharge = 0.0;
2813 for(Int_t it = 0; it < fTimeMax; it++){
2814 sumcharge += sum[it];
2815 if(sum[it] > 20.0) nbcl++;
2819 /////////////////////////////////////////////////////////
2821 ////////////////////////////////////////////////////////
2822 if(fDebugLevel > 0){
2823 if ( !fDebugStreamer ) {
2825 TDirectory *backup = gDirectory;
2826 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2827 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2830 Double_t amph0 = sum[0];
2831 Double_t amphlast = sum[fTimeMax-1];
2832 Double_t rms = TMath::RMS(fTimeMax,sum);
2833 Int_t goodtracklet = (Int_t) fGoodTracklet;
2834 for(Int_t it = 0; it < fTimeMax; it++){
2835 Double_t clustera = sum[it];
2837 (* fDebugStreamer) << "FillDAQa"<<
2838 "ampTotal="<<sumcharge<<
2841 "detector="<<idect<<
2843 "amphlast="<<amphlast<<
2844 "goodtracklet="<<goodtracklet<<
2845 "clustera="<<clustera<<
2852 ////////////////////////////////////////////////////////
2854 ///////////////////////////////////////////////////////
2855 if(sum[0] > 100.0) used = 1;
2856 if(nbcl < fNumberClusters) used = 1;
2857 if(nbcl > fNumberClustersf) used = 1;
2859 //if(fDetectorPreviousTrack == 15){
2860 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2862 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2864 for(Int_t it = 0; it < fTimeMax; it++){
2865 UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2869 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2871 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2878 //____________Online trackling in AliTRDtrigger________________________________
2879 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2883 // Fill a simple average pulse height
2887 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2893 //____________Write_____________________________________________________
2894 //_____________________________________________________________________
2895 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2898 // Write infos to file
2902 if ( fDebugStreamer ) {
2903 delete fDebugStreamer;
2904 fDebugStreamer = 0x0;
2907 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2912 ,fNumberUsedPh[1]));
2914 TDirectory *backup = gDirectory;
2920 option = "recreate";
2922 TFile f(filename,option.Data());
2924 TStopwatch stopwatch;
2927 f.WriteTObject(fCalibraVector);
2932 f.WriteTObject(fCH2d);
2937 f.WriteTObject(fPH2d);
2942 f.WriteTObject(fPRF2d);
2945 if(fLinearFitterOn){
2946 AnalyseLinearFitter();
2947 f.WriteTObject(fLinearVdriftFit);
2952 if ( backup ) backup->cd();
2954 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2955 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2957 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2959 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2960 //___________________________________________probe the histos__________________________________________________
2961 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2964 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2965 // debug mode with 2 for TH2I and 3 for TProfile2D
2966 // It gives a pointer to a Double_t[7] with the info following...
2967 // [0] : number of calibration groups with entries
2968 // [1] : minimal number of entries found
2969 // [2] : calibration group number of the min
2970 // [3] : maximal number of entries found
2971 // [4] : calibration group number of the max
2972 // [5] : mean number of entries found
2973 // [6] : mean relative error
2976 Double_t *info = new Double_t[7];
2978 // Number of Xbins (detectors or groups of pads)
2979 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2980 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2983 Double_t nbwe = 0; //number of calibration groups with entries
2984 Double_t minentries = 0; //minimal number of entries found
2985 Double_t maxentries = 0; //maximal number of entries found
2986 Double_t placemin = 0; //calibration group number of the min
2987 Double_t placemax = -1; //calibration group number of the max
2988 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2989 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2991 Double_t counter = 0;
2994 TH1F *nbEntries = 0x0;//distribution of the number of entries
2995 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2996 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2998 // Beginning of the loop over the calibration groups
2999 for (Int_t idect = 0; idect < nbins; idect++) {
3001 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
3002 projch->SetDirectory(0);
3004 // Number of entries for this calibration group
3005 Double_t nentries = 0.0;
3007 for (Int_t k = 0; k < nxbins; k++) {
3008 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
3012 for (Int_t k = 0; k < nxbins; k++) {
3013 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3014 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3015 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3023 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3024 nbEntries = new TH1F("Number of entries","Number of entries"
3025 ,100,(Int_t)nentries/2,nentries*2);
3026 nbEntries->SetDirectory(0);
3027 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3029 nbEntriesPerGroup->SetDirectory(0);
3030 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3031 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3032 nbEntriesPerSp->SetDirectory(0);
3035 if(nentries > 0) nbEntries->Fill(nentries);
3036 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3037 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3042 if(nentries > maxentries){
3043 maxentries = nentries;
3047 minentries = nentries;
3049 if(nentries < minentries){
3050 minentries = nentries;
3056 meanstats += nentries;
3058 }//calibration groups loop
3060 if(nbwe > 0) meanstats /= nbwe;
3061 if(counter > 0) meanrelativerror /= counter;
3063 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3064 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3065 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3066 AliInfo(Form("The mean number of entries is %f",meanstats));
3067 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3070 info[1] = minentries;
3072 info[3] = maxentries;
3074 info[5] = meanstats;
3075 info[6] = meanrelativerror;
3078 gStyle->SetPalette(1);
3079 gStyle->SetOptStat(1111);
3080 gStyle->SetPadBorderMode(0);
3081 gStyle->SetCanvasColor(10);
3082 gStyle->SetPadLeftMargin(0.13);
3083 gStyle->SetPadRightMargin(0.01);
3084 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3087 nbEntries->Draw("");
3089 nbEntriesPerSp->SetStats(0);
3090 nbEntriesPerSp->Draw("");
3091 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3093 nbEntriesPerGroup->SetStats(0);
3094 nbEntriesPerGroup->Draw("");
3100 //____________________________________________________________________________
3101 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3104 // Return a Int_t[4] with:
3105 // 0 Mean number of entries
3106 // 1 median of number of entries
3107 // 2 rms of number of entries
3108 // 3 number of group with entries
3111 Double_t *stat = new Double_t[4];
3114 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3115 Double_t *weight = new Double_t[nbofgroups];
3116 Int_t *nonul = new Int_t[nbofgroups];
3118 for(Int_t k = 0; k < nbofgroups; k++){
3119 if(fEntriesCH[k] > 0) {
3121 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3124 else weight[k] = 0.0;
3126 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3127 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3128 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3133 //____________________________________________________________________________
3134 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3137 // Return a Int_t[4] with:
3138 // 0 Mean number of entries
3139 // 1 median of number of entries
3140 // 2 rms of number of entries
3141 // 3 number of group with entries
3144 Double_t *stat = new Double_t[4];
3147 Int_t nbofgroups = 540;
3148 Double_t *weight = new Double_t[nbofgroups];
3149 Int_t *nonul = new Int_t[nbofgroups];
3151 for(Int_t k = 0; k < nbofgroups; k++){
3152 if(fEntriesLinearFitter[k] > 0) {
3154 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3157 else weight[k] = 0.0;
3159 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3160 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3161 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3166 //////////////////////////////////////////////////////////////////////////////////////
3168 //////////////////////////////////////////////////////////////////////////////////////
3169 //_____________________________________________________________________________
3170 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3173 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3174 // If fNgroupprf is zero then no binning in tnp
3178 name += fCalibraMode->GetNz(2);
3180 name += fCalibraMode->GetNrphi(2);
3184 if(fNgroupprf != 0){
3186 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3187 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3188 fPRF2d->SetYTitle("Det/pad groups");
3189 fPRF2d->SetXTitle("Position x/W [pad width units]");
3190 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3191 fPRF2d->SetStats(0);
3194 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3195 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3196 fPRF2d->SetYTitle("Det/pad groups");
3197 fPRF2d->SetXTitle("Position x/W [pad width units]");
3198 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3199 fPRF2d->SetStats(0);
3204 //_____________________________________________________________________________
3205 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3208 // Create the 2D histos
3212 name += fCalibraMode->GetNz(1);
3214 name += fCalibraMode->GetNrphi(1);
3216 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3217 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3219 fPH2d->SetYTitle("Det/pad groups");
3220 fPH2d->SetXTitle("time [#mus]");
3221 fPH2d->SetZTitle("<PH> [a.u.]");
3225 //_____________________________________________________________________________
3226 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3229 // Create the 2D histos
3233 name += fCalibraMode->GetNz(0);
3235 name += fCalibraMode->GetNrphi(0);
3237 fCH2d = new TH2I("CH2d",(const Char_t *) name
3238 ,fNumberBinCharge,0,300,nn,0,nn);
3239 fCH2d->SetYTitle("Det/pad groups");
3240 fCH2d->SetXTitle("charge deposit [a.u]");
3241 fCH2d->SetZTitle("counts");
3246 //////////////////////////////////////////////////////////////////////////////////
3247 // Set relative scale
3248 /////////////////////////////////////////////////////////////////////////////////
3249 //_____________________________________________________________________________
3250 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3253 // Set the factor that will divide the deposited charge
3254 // to fit in the histo range [0,300]
3257 if (RelativeScale > 0.0) {
3258 fRelativeScale = RelativeScale;
3261 AliInfo("RelativeScale must be strict positif!");
3265 //////////////////////////////////////////////////////////////////////////////////
3266 // Quick way to fill a histo
3267 //////////////////////////////////////////////////////////////////////////////////
3268 //_____________________________________________________________________
3269 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3272 // FillCH2d: Marian style
3275 //skip simply the value out of range
3276 if((y>=300.0) || (y<0.0)) return;
3278 //Calcul the y place
3279 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3280 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3283 fCH2d->GetArray()[place]++;
3287 //////////////////////////////////////////////////////////////////////////////////
3288 // Geometrical functions
3289 ///////////////////////////////////////////////////////////////////////////////////
3290 //_____________________________________________________________________________
3291 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3294 // Reconstruct the plane number from the detector number
3297 return ((Int_t) (d % 6));
3301 //_____________________________________________________________________________
3302 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3305 // Reconstruct the chamber number from the detector number
3309 return ((Int_t) (d % 30) / fgkNplan);
3313 //_____________________________________________________________________________
3314 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3317 // Reconstruct the sector number from the detector number
3321 return ((Int_t) (d / fg));
3324 ///////////////////////////////////////////////////////////////////////////////////
3325 // Getter functions for DAQ of the CH2d and the PH2d
3326 //////////////////////////////////////////////////////////////////////////////////
3327 //_____________________________________________________________________
3328 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3331 // return pointer to fPH2d TProfile2D
3332 // create a new TProfile2D if it doesn't exist allready
3338 fTimeMax = nbtimebin;
3339 fSf = samplefrequency;
3345 //_____________________________________________________________________
3346 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3349 // return pointer to fCH2d TH2I
3350 // create a new TH2I if it doesn't exist allready
3359 ////////////////////////////////////////////////////////////////////////////////////////////
3360 // Drift velocity calibration
3361 ///////////////////////////////////////////////////////////////////////////////////////////
3362 //_____________________________________________________________________
3363 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3366 // return pointer to TLinearFitter Calibration
3367 // if force is true create a new TLinearFitter if it doesn't exist allready
3370 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3371 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3374 // if we are forced and TLinearFitter doesn't yet exist create it
3376 // new TLinearFitter
3377 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3378 fLinearFitterArray.AddAt(linearfitter,detector);
3379 return linearfitter;
3382 //____________________________________________________________________________
3383 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3386 // Analyse array of linear fitter because can not be written
3387 // Store two arrays: one with the param the other one with the error param + number of entries
3390 for(Int_t k = 0; k < 540; k++){
3391 TLinearFitter *linearfitter = GetLinearFitter(k);
3392 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3393 TVectorD *par = new TVectorD(2);
3394 TVectorD pare = TVectorD(2);
3395 TVectorD *parE = new TVectorD(3);
3396 linearfitter->Eval();
3397 linearfitter->GetParameters(*par);
3398 linearfitter->GetErrors(pare);
3399 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3400 (*parE)[0] = pare[0]*ppointError;
3401 (*parE)[1] = pare[1]*ppointError;
3402 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3403 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3404 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);