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 "AliTRDtrack.h"
66 #include "AliTRDtrackV1.h"
67 #include "AliTRDrawStreamBase.h"
68 #include "AliRawReader.h"
69 #include "AliRawReaderDate.h"
70 #include "AliTRDgeometry.h"
71 #include "./Cal/AliTRDCalROC.h"
72 #include "./Cal/AliTRDCalDet.h"
79 ClassImp(AliTRDCalibraFillHisto)
81 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
82 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
84 //_____________singleton implementation_________________________________________________
85 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
88 // Singleton implementation
91 if (fgTerminated != kFALSE) {
95 if (fgInstance == 0) {
96 fgInstance = new AliTRDCalibraFillHisto();
103 //______________________________________________________________________________________
104 void AliTRDCalibraFillHisto::Terminate()
107 // Singleton implementation
108 // Deletes the instance of this class
111 fgTerminated = kTRUE;
113 if (fgInstance != 0) {
120 //______________________________________________________________________________________
121 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
124 ,fMcmCorrectAngle(kFALSE)
130 ,fLinearFitterOn(kFALSE)
131 ,fLinearFitterDebugOn(kFALSE)
133 ,fThresholdClusterPRF2(15.0)
134 ,fLimitChargeIntegration(kFALSE)
135 ,fFillWithZero(kFALSE)
136 ,fCalibraMode(new AliTRDCalibraMode())
139 ,fDetectorPreviousTrack(-1)
143 ,fNumberClustersf(30)
149 ,fNumberBinCharge(100)
155 ,fGoodTracklet(kTRUE)
156 ,fLinearFitterTracklet(0x0)
158 ,fEntriesLinearFitter(0x0)
163 ,fLinearFitterArray(540)
164 ,fLinearVdriftFit(0x0)
169 // Default constructor
173 // Init some default values
176 fNumberUsedCh[0] = 0;
177 fNumberUsedCh[1] = 0;
178 fNumberUsedPh[0] = 0;
179 fNumberUsedPh[1] = 0;
181 fGeo = new AliTRDgeometry();
185 //______________________________________________________________________________________
186 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
189 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
192 ,fPRF2dOn(c.fPRF2dOn)
193 ,fHisto2d(c.fHisto2d)
194 ,fVector2d(c.fVector2d)
195 ,fLinearFitterOn(c.fLinearFitterOn)
196 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
197 ,fRelativeScale(c.fRelativeScale)
198 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
199 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
200 ,fFillWithZero(c.fFillWithZero)
203 ,fDebugLevel(c.fDebugLevel)
204 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
205 ,fMCMPrevious(c.fMCMPrevious)
206 ,fROBPrevious(c.fROBPrevious)
207 ,fNumberClusters(c.fNumberClusters)
208 ,fNumberClustersf(c.fNumberClustersf)
209 ,fProcent(c.fProcent)
210 ,fDifference(c.fDifference)
211 ,fNumberTrack(c.fNumberTrack)
212 ,fTimeMax(c.fTimeMax)
214 ,fNumberBinCharge(c.fNumberBinCharge)
215 ,fNumberBinPRF(c.fNumberBinPRF)
216 ,fNgroupprf(c.fNgroupprf)
220 ,fGoodTracklet(c.fGoodTracklet)
221 ,fLinearFitterTracklet(0x0)
223 ,fEntriesLinearFitter(0x0)
228 ,fLinearFitterArray(540)
229 ,fLinearVdriftFit(0x0)
236 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
237 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
239 fPH2d = new TProfile2D(*c.fPH2d);
240 fPH2d->SetDirectory(0);
243 fPRF2d = new TProfile2D(*c.fPRF2d);
244 fPRF2d->SetDirectory(0);
247 fCH2d = new TH2I(*c.fCH2d);
248 fCH2d->SetDirectory(0);
250 if(c.fLinearVdriftFit){
251 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
254 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
255 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
260 fGeo = new AliTRDgeometry();
263 //____________________________________________________________________________________
264 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
267 // AliTRDCalibraFillHisto destructor
271 if ( fDebugStreamer ) delete fDebugStreamer;
273 if ( fCalDetGain ) delete fCalDetGain;
274 if ( fCalROCGain ) delete fCalROCGain;
276 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
280 delete [] fEntriesCH;
281 delete [] fEntriesLinearFitter;
284 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
285 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
293 //_____________________________________________________________________________
294 void AliTRDCalibraFillHisto::Destroy()
305 //_____________________________________________________________________________
306 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
309 // Delete DebugStreamer
312 if ( fDebugStreamer ) delete fDebugStreamer;
315 //_____________________________________________________________________________
316 void AliTRDCalibraFillHisto::ClearHistos()
336 //////////////////////////////////////////////////////////////////////////////////
337 // calibration with AliTRDtrackV1: Init, Update
338 //////////////////////////////////////////////////////////////////////////////////
339 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
340 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
343 // Init the histograms and stuff to be filled
348 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
350 AliInfo("Could not get calibDB");
353 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
355 AliInfo("Could not get CommonParam");
360 fTimeMax = cal->GetNumberOfTimeBins();
361 fSf = parCom->GetSamplingFrequency();
363 fNumberClustersf = fTimeMax;
364 fNumberClusters = (Int_t)(0.5*fTimeMax);
366 // Init linear fitter
367 if(!fLinearFitterTracklet) {
368 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
369 fLinearFitterTracklet->StoreData(kFALSE);
372 //calib object from database used for reconstruction
374 fCalDetGain->~AliTRDCalDet();
375 new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
376 }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
378 // Calcul Xbins Chambd0, Chamb2
379 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
380 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
381 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
383 // If vector method On initialised all the stuff
385 fCalibraVector = new AliTRDCalibraVector();
386 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
387 fCalibraVector->SetTimeMax(fTimeMax);
388 if(fNgroupprf != 0) {
389 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
390 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
393 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
394 fCalibraVector->SetPRFRange(1.5);
396 for(Int_t k = 0; k < 3; k++){
397 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
398 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
400 TString namech("Nz");
401 namech += fCalibraMode->GetNz(0);
403 namech += fCalibraMode->GetNrphi(0);
404 fCalibraVector->SetNameCH((const char* ) namech);
405 TString nameph("Nz");
406 nameph += fCalibraMode->GetNz(1);
408 nameph += fCalibraMode->GetNrphi(1);
409 fCalibraVector->SetNamePH((const char* ) nameph);
410 TString nameprf("Nz");
411 nameprf += fCalibraMode->GetNz(2);
413 nameprf += fCalibraMode->GetNrphi(2);
415 nameprf += fNgroupprf;
416 fCalibraVector->SetNamePRF((const char* ) nameprf);
419 // Create the 2D histos corresponding to the pad groupCalibration mode
422 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
423 ,fCalibraMode->GetNz(0)
424 ,fCalibraMode->GetNrphi(0)));
426 // Create the 2D histo
431 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
432 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
436 fEntriesCH = new Int_t[ntotal0];
437 for(Int_t k = 0; k < ntotal0; k++){
444 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
445 ,fCalibraMode->GetNz(1)
446 ,fCalibraMode->GetNrphi(1)));
448 // Create the 2D histo
453 fPHPlace = new Short_t[fTimeMax];
454 for (Int_t k = 0; k < fTimeMax; k++) {
457 fPHValue = new Float_t[fTimeMax];
458 for (Int_t k = 0; k < fTimeMax; k++) {
462 if (fLinearFitterOn) {
463 //fLinearFitterArray.Expand(540);
464 fLinearFitterArray.SetName("ArrayLinearFitters");
465 fEntriesLinearFitter = new Int_t[540];
466 for(Int_t k = 0; k < 540; k++){
467 fEntriesLinearFitter[k] = 0;
469 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
474 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
475 ,fCalibraMode->GetNz(2)
476 ,fCalibraMode->GetNrphi(2)));
477 // Create the 2D histo
479 CreatePRF2d(ntotal2);
486 //____________Offline tracking in the AliTRDtracker____________________________
487 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
490 // Use AliTRDtrack for the calibration
494 AliTRDcluster *cl = 0x0; // pointeur to cluster
495 Int_t index0 = 0; // index of the first cluster in the current chamber
496 Int_t index1 = 0; // index of the current cluster in the current chamber
498 //////////////////////////////
499 // loop over the clusters
500 ///////////////////////////////
501 while((cl = t->GetCluster(index1))){
503 /////////////////////////////////////////////////////////////////////////
504 // Fill the infos for the previous clusters if not the same detector
505 ////////////////////////////////////////////////////////////////////////
506 Int_t detector = cl->GetDetector();
507 if ((detector != fDetectorPreviousTrack) &&
508 (index0 != index1)) {
512 //If the same track, then look if the previous detector is in
513 //the same plane, if yes: not a good track
514 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
518 // Fill only if the track doesn't touch a masked pad or doesn't
521 // drift velocity unables to cut bad tracklets
522 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
526 FillTheInfoOfTheTrackCH(index1-index0);
531 FillTheInfoOfTheTrackPH();
534 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
537 } // if a good tracklet
540 ResetfVariablestracklet();
543 } // Fill at the end the charge
546 /////////////////////////////////////////////////////////////
547 // Localise and take the calib gain object for the detector
548 ////////////////////////////////////////////////////////////
549 if (detector != fDetectorPreviousTrack) {
551 //Localise the detector
552 LocalisationDetectorXbins(detector);
555 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
557 AliInfo("Could not get calibDB");
563 fCalROCGain->~AliTRDCalROC();
564 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
565 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
569 // Reset the detectbjobsor
570 fDetectorPreviousTrack = detector;
572 ///////////////////////////////////////
573 // Store the info of the cluster
574 ///////////////////////////////////////
575 Int_t row = cl->GetPadRow();
576 Int_t col = cl->GetPadCol();
577 CheckGoodTrackletV1(cl);
578 Int_t group[2] = {0,0};
579 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
580 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
581 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
585 } // while on clusters
587 ///////////////////////////
588 // Fill the last plane
589 //////////////////////////
590 if( index0 != index1 ){
596 // drift velocity unables to cut bad tracklets
597 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
601 FillTheInfoOfTheTrackCH(index1-index0);
606 FillTheInfoOfTheTrackPH();
609 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
611 } // if a good tracklet
616 ResetfVariablestracklet();
621 //____________Offline tracking in the AliTRDtracker____________________________
622 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
625 // Use AliTRDtrackV1 for the calibration
629 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
630 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
631 Bool_t newtr = kTRUE; // new track
634 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
636 AliInfo("Could not get calibDB");
640 ///////////////////////////
641 // loop over the tracklet
642 ///////////////////////////
643 for(Int_t itr = 0; itr < 6; itr++){
645 if(!(tracklet = t->GetTracklet(itr))) continue;
646 if(!tracklet->IsOK()) continue;
648 ResetfVariablestracklet();
650 //////////////////////////////////////////
651 // localisation of the tracklet and dqdl
652 //////////////////////////////////////////
653 Int_t layer = tracklet->GetPlane();
655 while(!(cl = tracklet->GetClusters(ic++))) continue;
656 Int_t detector = cl->GetDetector();
657 if (detector != fDetectorPreviousTrack) {
658 // if not a new track
660 // don't use the rest of this track if in the same plane
661 if (layer == GetLayer(fDetectorPreviousTrack)) {
662 //printf("bad tracklet, same layer for detector %d\n",detector);
666 //Localise the detector bin
667 LocalisationDetectorXbins(detector);
670 fCalROCGain->~AliTRDCalROC();
671 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
672 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
675 fDetectorPreviousTrack = detector;
679 ////////////////////////////
680 // loop over the clusters
681 ////////////////////////////
682 Int_t nbclusters = 0;
683 for(int jc=0; jc<AliTRDseed::knTimebins; jc++){
684 if(!(cl = tracklet->GetClusters(jc))) continue;
687 // Store the info bis of the tracklet
688 Int_t row = cl->GetPadRow();
689 Int_t col = cl->GetPadCol();
690 CheckGoodTrackletV1(cl);
691 Int_t group[2] = {0,0};
692 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
693 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
694 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col);
697 ////////////////////////////////////////
698 // Fill the stuffs if a good tracklet
699 ////////////////////////////////////////
702 // drift velocity unables to cut bad tracklets
703 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
707 FillTheInfoOfTheTrackCH(nbclusters);
712 FillTheInfoOfTheTrackPH();
715 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
717 } // if a good tracklet
723 ///////////////////////////////////////////////////////////////////////////////////
724 // Routine inside the update with AliTRDtrack
725 ///////////////////////////////////////////////////////////////////////////////////
726 //____________Offine tracking in the AliTRDtracker_____________________________
727 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
730 // Drift velocity calibration:
731 // Fit the clusters with a straight line
732 // From the slope find the drift velocity
735 //Number of points: if less than 3 return kFALSE
736 Int_t npoints = index1-index0;
737 if(npoints <= 2) return kFALSE;
742 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
743 // parameters of the track
744 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
745 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
746 Double_t tnp = 0.0; // tan angle in the plan xy track
747 if( TMath::Abs(snp) < 1.){
748 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
750 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
751 // tilting pad and cross row
752 Int_t crossrow = 0; // if it crosses a pad row
753 Int_t rowp = -1; // if it crosses a pad row
754 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
755 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
756 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
758 fLinearFitterTracklet->ClearPoints();
759 Double_t dydt = 0.0; // dydt tracklet after straight line fit
760 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
761 Double_t pointError = 0.0; // error after straight line fit
762 Int_t nbli = 0; // number linear fitter points
764 //////////////////////////////
765 // loop over clusters
766 ////////////////////////////
767 for(Int_t k = 0; k < npoints; k++){
769 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
770 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
771 Double_t ycluster = cl->GetY();
772 Int_t time = cl->GetPadTime();
773 Double_t timeis = time/fSf;
774 //Double_t q = cl->GetQ();
775 //See if cross two pad rows
776 Int_t row = cl->GetPadRow();
778 if(row != rowp) crossrow = 1;
780 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
785 //////////////////////////////
787 //////////////////////////////
789 fLinearFitterTracklet->ClearPoints();
793 fLinearFitterTracklet->Eval();
794 fLinearFitterTracklet->GetParameters(pars);
795 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
796 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
798 fLinearFitterTracklet->ClearPoints();
800 /////////////////////////////
802 ////////////////////////////
804 if ( !fDebugStreamer ) {
806 TDirectory *backup = gDirectory;
807 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
808 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
812 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
813 //"snpright="<<snpright<<
814 "npoints="<<npoints<<
816 "detector="<<detector<<
823 "crossrow="<<crossrow<<
824 "errorpar="<<errorpar<<
825 "pointError="<<pointError<<
829 Int_t nbclusters = index1-index0;
830 Int_t layer = GetLayer(fDetectorPreviousTrack);
832 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
833 //"snpright="<<snpright<<
834 "nbclusters="<<nbclusters<<
835 "detector="<<fDetectorPreviousTrack<<
841 ///////////////////////////
843 ///////////////////////////
844 if(npoints < fNumberClusters) return kFALSE;
845 if(npoints > fNumberClustersf) return kFALSE;
846 if(pointError >= 0.1) return kFALSE;
847 if(crossrow == 1) return kFALSE;
849 ////////////////////////////
851 ////////////////////////////
853 //Add to the linear fitter of the detector
854 if( TMath::Abs(snp) < 1.){
855 Double_t x = tnp-dzdx*tnt;
856 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
857 if(fLinearFitterDebugOn) {
858 fLinearVdriftFit->Update(detector,x,pars[1]);
860 fEntriesLinearFitter[detector]++;
866 //____________Offine tracking in the AliTRDtracker_____________________________
867 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
870 // Drift velocity calibration:
871 // Fit the clusters with a straight line
872 // From the slope find the drift velocity
875 ////////////////////////////////////////////////
876 //Number of points: if less than 3 return kFALSE
877 /////////////////////////////////////////////////
878 if(nbclusters <= 2) return kFALSE;
883 // results of the linear fit
884 Double_t dydt = 0.0; // dydt tracklet after straight line fit
885 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
886 Double_t pointError = 0.0; // error after straight line fit
887 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
888 Int_t crossrow = 0; // if it crosses a pad row
889 Int_t rowp = -1; // if it crosses a pad row
890 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
891 fLinearFitterTracklet->ClearPoints();
894 ///////////////////////////////////////////
895 // Take the parameters of the track
896 //////////////////////////////////////////
897 // take now the snp, tnp and tgl from the track
898 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
899 Double_t tnp = 0.0; // dy/dx at the end of the chamber
900 if( TMath::Abs(snp) < 1.){
901 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
903 Double_t tgl = tracklet->GetTgl(); // dz/dl
904 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
906 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
907 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
908 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
909 // at the end with correction due to linear fit
910 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
911 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
914 ////////////////////////////
915 // loop over the clusters
916 ////////////////////////////
918 AliTRDcluster *cl = 0x0;
919 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
920 if(!(cl = tracklet->GetClusters(ic))) continue;
921 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
923 Double_t ycluster = cl->GetY();
924 Int_t time = cl->GetPadTime();
925 Double_t timeis = time/fSf;
926 //See if cross two pad rows
927 Int_t row = cl->GetPadRow();
928 if(rowp==-1) rowp = row;
929 if(row != rowp) crossrow = 1;
931 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
936 ////////////////////////////////////
937 // Do the straight line fit now
938 ///////////////////////////////////
940 fLinearFitterTracklet->ClearPoints();
944 fLinearFitterTracklet->Eval();
945 fLinearFitterTracklet->GetParameters(pars);
946 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
947 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
949 fLinearFitterTracklet->ClearPoints();
951 ////////////////////////////////
953 ///////////////////////////////
957 if ( !fDebugStreamer ) {
959 TDirectory *backup = gDirectory;
960 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
961 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
965 Int_t layer = GetLayer(fDetectorPreviousTrack);
967 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
968 //"snpright="<<snpright<<
969 "nbclusters="<<nbclusters<<
970 "detector="<<fDetectorPreviousTrack<<
978 "crossrow="<<crossrow<<
979 "errorpar="<<errorpar<<
980 "pointError="<<pointError<<
985 /////////////////////////
987 ////////////////////////
989 if(nbclusters < fNumberClusters) return kFALSE;
990 if(nbclusters > fNumberClustersf) return kFALSE;
991 if(pointError >= 0.1) return kFALSE;
992 if(crossrow == 1) return kFALSE;
994 ///////////////////////
996 //////////////////////
999 //Add to the linear fitter of the detector
1000 if( TMath::Abs(snp) < 1.){
1001 Double_t x = tnp-dzdx*tnt;
1002 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1003 if(fLinearFitterDebugOn) {
1004 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1006 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1012 //____________Offine tracking in the AliTRDtracker_____________________________
1013 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1016 // PRF width calibration
1017 // Assume a Gaussian shape: determinate the position of the three pad clusters
1018 // Fit with a straight line
1019 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1020 // Fill the PRF as function of angle of the track
1025 //////////////////////////
1027 /////////////////////////
1028 Int_t npoints = index1-index0; // number of total points
1029 Int_t nb3pc = 0; // number of three pads clusters used for fit
1030 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1031 // To see the difference due to the fit
1032 Double_t *padPositions;
1033 padPositions = new Double_t[npoints];
1034 for(Int_t k = 0; k < npoints; k++){
1035 padPositions[k] = 0.0;
1037 // Take the tgl and snp with the track t now
1038 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1039 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1040 Float_t dzdx = 0.0; // dzdx
1042 if(TMath::Abs(snp) < 1.0){
1043 tnp = snp / (TMath::Sqrt(1-snp*snp));
1044 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1047 fLinearFitterTracklet->ClearPoints();
1049 ///////////////////////////
1050 // calcul the tnp group
1051 ///////////////////////////
1052 Bool_t echec = kFALSE;
1053 Double_t shift = 0.0;
1054 //Calculate the shift in x coresponding to this tnp
1055 if(fNgroupprf != 0.0){
1056 shift = -3.0*(fNgroupprf-1)-1.5;
1057 Double_t limithigh = -0.2*(fNgroupprf-1);
1058 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1060 while(tnp > limithigh){
1067 delete [] padPositions;
1071 //////////////////////
1073 /////////////////////
1074 for(Int_t k = 0; k < npoints; k++){
1076 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1077 Short_t *signals = cl->GetSignals();
1078 Double_t time = cl->GetLocalTimeBin();
1079 //Calculate x if possible
1080 Float_t xcenter = 0.0;
1081 Bool_t echec1 = kTRUE;
1082 if((time<=7) || (time>=21)) continue;
1083 // Center 3 balanced: position with the center of the pad
1084 if ((((Float_t) signals[3]) > 0.0) &&
1085 (((Float_t) signals[2]) > 0.0) &&
1086 (((Float_t) signals[4]) > 0.0)) {
1088 // Security if the denomiateur is 0
1089 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1090 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1091 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1092 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1093 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1099 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1101 //if no echec: calculate with the position of the pad
1102 // Position of the cluster
1103 Double_t padPosition = xcenter + cl->GetPadCol();
1104 padPositions[k] = padPosition;
1106 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1110 /////////////////////////////
1112 ////////////////////////////
1114 delete [] padPositions;
1115 fLinearFitterTracklet->ClearPoints();
1118 fLinearFitterTracklet->Eval();
1120 fLinearFitterTracklet->GetParameters(line);
1121 Float_t pointError = -1.0;
1122 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1123 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1125 fLinearFitterTracklet->ClearPoints();
1127 /////////////////////////////////////////////////////
1128 // Now fill the PRF: second loop over clusters
1129 /////////////////////////////////////////////////////
1130 for(Int_t k = 0; k < npoints; k++){
1132 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1133 Short_t *signals = cl->GetSignals(); // signal
1134 Double_t time = cl->GetLocalTimeBin(); // time bin
1135 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1136 Float_t padPos = cl->GetPadCol(); // middle pad
1137 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1138 Float_t ycenter = 0.0; // relative center charge
1139 Float_t ymin = 0.0; // relative left charge
1140 Float_t ymax = 0.0; // relative right charge
1142 //Requiere simply two pads clusters at least
1143 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1144 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1145 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1146 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1147 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1148 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1152 Int_t rowcl = cl->GetPadRow(); // row of cluster
1153 Int_t colcl = cl->GetPadCol(); // col of cluster
1154 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1155 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1156 Float_t xcl = cl->GetY(); // y cluster
1157 Float_t qcl = cl->GetQ(); // charge cluster
1158 Int_t layer = GetLayer(detector); // layer
1159 Int_t stack = GetStack(detector); // stack
1160 Double_t xdiff = dpad; // reconstructed position constant
1161 Double_t x = dpad; // reconstructed position moved
1162 Float_t ep = pointError; // error of fit
1163 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1164 Float_t signal3 = (Float_t)signals[3]; // signal
1165 Float_t signal2 = (Float_t)signals[2]; // signal
1166 Float_t signal4 = (Float_t)signals[4]; // signal
1167 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1169 //////////////////////////////
1171 /////////////////////////////
1172 if(fDebugLevel > 0){
1173 if ( !fDebugStreamer ) {
1175 TDirectory *backup = gDirectory;
1176 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1177 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1183 Float_t y = ycenter;
1184 (* fDebugStreamer) << "HandlePRFtracklet"<<
1185 "caligroup="<<caligroup<<
1186 "detector="<<detector<<
1189 "npoints="<<npoints<<
1198 "padPosition="<<padPositions[k]<<
1199 "padPosTracklet="<<padPosTracklet<<
1204 "signal1="<<signal1<<
1205 "signal2="<<signal2<<
1206 "signal3="<<signal3<<
1207 "signal4="<<signal4<<
1208 "signal5="<<signal5<<
1214 (* fDebugStreamer) << "HandlePRFtracklet"<<
1215 "caligroup="<<caligroup<<
1216 "detector="<<detector<<
1219 "npoints="<<npoints<<
1228 "padPosition="<<padPositions[k]<<
1229 "padPosTracklet="<<padPosTracklet<<
1234 "signal1="<<signal1<<
1235 "signal2="<<signal2<<
1236 "signal3="<<signal3<<
1237 "signal4="<<signal4<<
1238 "signal5="<<signal5<<
1244 (* fDebugStreamer) << "HandlePRFtracklet"<<
1245 "caligroup="<<caligroup<<
1246 "detector="<<detector<<
1249 "npoints="<<npoints<<
1258 "padPosition="<<padPositions[k]<<
1259 "padPosTracklet="<<padPosTracklet<<
1264 "signal1="<<signal1<<
1265 "signal2="<<signal2<<
1266 "signal3="<<signal3<<
1267 "signal4="<<signal4<<
1268 "signal5="<<signal5<<
1274 ////////////////////////////
1276 ///////////////////////////
1277 if(npoints < fNumberClusters) continue;
1278 if(npoints > fNumberClustersf) continue;
1279 if(nb3pc <= 5) continue;
1280 if((time >= 21) || (time < 7)) continue;
1281 if(TMath::Abs(snp) >= 1.0) continue;
1282 if(TMath::Abs(qcl) < 80) continue;
1284 ////////////////////////////
1286 ///////////////////////////
1288 if(TMath::Abs(dpad) < 1.5) {
1289 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1290 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1292 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1293 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1294 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1296 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1297 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1298 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1302 if(TMath::Abs(dpad) < 1.5) {
1303 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1304 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1306 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1307 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1308 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1310 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1311 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1312 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1316 delete [] padPositions;
1320 //____________Offine tracking in the AliTRDtracker_____________________________
1321 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1324 // PRF width calibration
1325 // Assume a Gaussian shape: determinate the position of the three pad clusters
1326 // Fit with a straight line
1327 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1328 // Fill the PRF as function of angle of the track
1332 //printf("begin\n");
1333 ///////////////////////////////////////////
1334 // Take the parameters of the track
1335 //////////////////////////////////////////
1336 // take now the snp, tnp and tgl from the track
1337 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1338 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1339 if( TMath::Abs(snp) < 1.){
1340 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1342 Double_t tgl = tracklet->GetTgl(); // dz/dl
1343 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1345 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1346 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1347 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1348 // at the end with correction due to linear fit
1349 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1350 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1352 ///////////////////////////////
1353 // Calculate tnp group shift
1354 ///////////////////////////////
1355 Bool_t echec = kFALSE;
1356 Double_t shift = 0.0;
1357 //Calculate the shift in x coresponding to this tnp
1358 if(fNgroupprf != 0.0){
1359 shift = -3.0*(fNgroupprf-1)-1.5;
1360 Double_t limithigh = -0.2*(fNgroupprf-1);
1361 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1363 while(tnp > limithigh){
1369 // do nothing if out of tnp range
1370 //printf("echec %d\n",(Int_t)echec);
1371 if(echec) return kFALSE;
1373 ///////////////////////
1375 //////////////////////
1377 Int_t nb3pc = 0; // number of three pads clusters used for fit
1378 Double_t *padPositions; // to see the difference between the fit and the 3 pad clusters position
1379 padPositions = new Double_t[AliTRDseed::knTimebins];
1380 for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1381 padPositions[k] = 0.0;
1383 fLinearFitterTracklet->ClearPoints();
1385 //printf("loop clusters \n");
1386 ////////////////////////////
1387 // loop over the clusters
1388 ////////////////////////////
1389 AliTRDcluster *cl = 0x0;
1390 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1391 if(!(cl = tracklet->GetClusters(ic))) continue;
1393 Double_t time = cl->GetLocalTimeBin();
1394 if((time<=7) || (time>=21)) continue;
1395 Short_t *signals = cl->GetSignals();
1396 Float_t xcenter = 0.0;
1397 Bool_t echec1 = kTRUE;
1399 /////////////////////////////////////////////////////////////
1400 // Center 3 balanced: position with the center of the pad
1401 /////////////////////////////////////////////////////////////
1402 if ((((Float_t) signals[3]) > 0.0) &&
1403 (((Float_t) signals[2]) > 0.0) &&
1404 (((Float_t) signals[4]) > 0.0)) {
1406 // Security if the denomiateur is 0
1407 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1408 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1409 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1410 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1411 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1417 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1418 if(echec1) continue;
1420 ////////////////////////////////////////////////////////
1421 //if no echec1: calculate with the position of the pad
1422 // Position of the cluster
1423 // fill the linear fitter
1424 ///////////////////////////////////////////////////////
1425 Double_t padPosition = xcenter + cl->GetPadCol();
1426 padPositions[ic] = padPosition;
1428 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1433 //printf("Fin loop clusters \n");
1434 //////////////////////////////
1435 // fit with a straight line
1436 /////////////////////////////
1438 delete [] padPositions;
1439 fLinearFitterTracklet->ClearPoints();
1440 delete [] padPositions;
1443 fLinearFitterTracklet->Eval();
1445 fLinearFitterTracklet->GetParameters(line);
1446 Float_t pointError = -1.0;
1447 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1448 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1450 fLinearFitterTracklet->ClearPoints();
1452 //printf("PRF second loop \n");
1453 ////////////////////////////////////////////////
1454 // Fill the PRF: Second loop over clusters
1455 //////////////////////////////////////////////
1456 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1457 if(!(cl = tracklet->GetClusters(ic))) continue;
1459 Short_t *signals = cl->GetSignals(); // signal
1460 Double_t time = cl->GetLocalTimeBin(); // time bin
1461 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1462 Float_t padPos = cl->GetPadCol(); // middle pad
1463 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1464 Float_t ycenter = 0.0; // relative center charge
1465 Float_t ymin = 0.0; // relative left charge
1466 Float_t ymax = 0.0; // relative right charge
1468 ////////////////////////////////////////////////////////////////
1469 // Calculate ycenter, ymin and ymax even for two pad clusters
1470 ////////////////////////////////////////////////////////////////
1471 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1472 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1473 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1474 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1475 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1476 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1479 /////////////////////////
1480 // Calibration group
1481 ////////////////////////
1482 Int_t rowcl = cl->GetPadRow(); // row of cluster
1483 Int_t colcl = cl->GetPadCol(); // col of cluster
1484 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1485 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1486 Float_t xcl = cl->GetY(); // y cluster
1487 Float_t qcl = cl->GetQ(); // charge cluster
1488 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1489 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1490 Double_t xdiff = dpad; // reconstructed position constant
1491 Double_t x = dpad; // reconstructed position moved
1492 Float_t ep = pointError; // error of fit
1493 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1494 Float_t signal3 = (Float_t)signals[3]; // signal
1495 Float_t signal2 = (Float_t)signals[2]; // signal
1496 Float_t signal4 = (Float_t)signals[4]; // signal
1497 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1501 /////////////////////
1503 ////////////////////
1505 if(fDebugLevel > 0){
1506 if ( !fDebugStreamer ) {
1508 TDirectory *backup = gDirectory;
1509 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1510 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1515 Float_t y = ycenter;
1516 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1517 "caligroup="<<caligroup<<
1518 "detector="<<fDetectorPreviousTrack<<
1521 "npoints="<<nbclusters<<
1530 "padPosition="<<padPositions[ic]<<
1531 "padPosTracklet="<<padPosTracklet<<
1536 "signal1="<<signal1<<
1537 "signal2="<<signal2<<
1538 "signal3="<<signal3<<
1539 "signal4="<<signal4<<
1540 "signal5="<<signal5<<
1546 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1547 "caligroup="<<caligroup<<
1548 "detector="<<fDetectorPreviousTrack<<
1551 "npoints="<<nbclusters<<
1560 "padPosition="<<padPositions[ic]<<
1561 "padPosTracklet="<<padPosTracklet<<
1566 "signal1="<<signal1<<
1567 "signal2="<<signal2<<
1568 "signal3="<<signal3<<
1569 "signal4="<<signal4<<
1570 "signal5="<<signal5<<
1576 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1577 "caligroup="<<caligroup<<
1578 "detector="<<fDetectorPreviousTrack<<
1581 "npoints="<<nbclusters<<
1590 "padPosition="<<padPositions[ic]<<
1591 "padPosTracklet="<<padPosTracklet<<
1596 "signal1="<<signal1<<
1597 "signal2="<<signal2<<
1598 "signal3="<<signal3<<
1599 "signal4="<<signal4<<
1600 "signal5="<<signal5<<
1606 /////////////////////
1608 /////////////////////
1609 if(nbclusters < fNumberClusters) continue;
1610 if(nbclusters > fNumberClustersf) continue;
1611 if(nb3pc <= 5) continue;
1612 if((time >= 21) || (time < 7)) continue;
1613 if(TMath::Abs(qcl) < 80) continue;
1614 if( TMath::Abs(snp) > 1.) continue;
1617 ////////////////////////
1619 ///////////////////////
1621 if(TMath::Abs(dpad) < 1.5) {
1622 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1623 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1624 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1626 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1627 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1628 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1630 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1631 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1632 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1637 if(TMath::Abs(dpad) < 1.5) {
1638 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1639 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1641 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1642 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1643 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1645 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1646 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1647 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1650 } // second loop over clusters
1653 delete [] padPositions;
1657 ///////////////////////////////////////////////////////////////////////////////////////
1658 // Pad row col stuff: see if masked or not
1659 ///////////////////////////////////////////////////////////////////////////////////////
1660 //_____________________________________________________________________________
1661 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1664 // See if we are not near a masked pad
1667 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1671 //_____________________________________________________________________________
1672 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1675 // See if we are not near a masked pad
1678 if (!IsPadOn(detector, col, row)) {
1679 fGoodTracklet = kFALSE;
1683 if (!IsPadOn(detector, col-1, row)) {
1684 fGoodTracklet = kFALSE;
1689 if (!IsPadOn(detector, col+1, row)) {
1690 fGoodTracklet = kFALSE;
1695 //_____________________________________________________________________________
1696 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1699 // Look in the choosen database if the pad is On.
1700 // If no the track will be "not good"
1703 // Get the parameter object
1704 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1706 AliInfo("Could not get calibDB");
1710 if (!cal->IsChamberInstalled(detector) ||
1711 cal->IsChamberMasked(detector) ||
1712 cal->IsPadMasked(detector,col,row)) {
1720 ///////////////////////////////////////////////////////////////////////////////////////
1721 // Calibration groups: calculate the number of groups, localise...
1722 ////////////////////////////////////////////////////////////////////////////////////////
1723 //_____________________________________________________________________________
1724 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1727 // Calculate the calibration group number for i
1730 // Row of the cluster and position in the pad groups
1732 if (fCalibraMode->GetNnZ(i) != 0) {
1733 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1737 // Col of the cluster and position in the pad groups
1739 if (fCalibraMode->GetNnRphi(i) != 0) {
1740 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1743 return posc*fCalibraMode->GetNfragZ(i)+posr;
1746 //____________________________________________________________________________________
1747 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1750 // Calculate the total number of calibration groups
1754 fCalibraMode->ModePadCalibration(2,i);
1755 fCalibraMode->ModePadFragmentation(0,2,0,i);
1756 fCalibraMode->SetDetChamb2(i);
1757 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1758 fCalibraMode->ModePadCalibration(0,i);
1759 fCalibraMode->ModePadFragmentation(0,0,0,i);
1760 fCalibraMode->SetDetChamb0(i);
1761 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1762 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1766 //_____________________________________________________________________________
1767 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1770 // Set the mode of calibration group in the z direction for the parameter i
1775 fCalibraMode->SetNz(i, Nz);
1778 AliInfo("You have to choose between 0 and 4");
1782 //_____________________________________________________________________________
1783 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1786 // Set the mode of calibration group in the rphi direction for the parameter i
1791 fCalibraMode->SetNrphi(i ,Nrphi);
1794 AliInfo("You have to choose between 0 and 6");
1798 //____________Set the pad calibration variables for the detector_______________
1799 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1802 // For the detector calcul the first Xbins and set the number of row
1803 // and col pads per calibration groups, the number of calibration
1804 // groups in the detector.
1807 // first Xbins of the detector
1809 fCalibraMode->CalculXBins(detector,0);
1812 fCalibraMode->CalculXBins(detector,1);
1815 fCalibraMode->CalculXBins(detector,2);
1818 // fragmentation of idect
1819 for (Int_t i = 0; i < 3; i++) {
1820 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1821 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1822 , (Int_t) GetStack(detector)
1823 , (Int_t) GetSector(detector),i);
1829 //_____________________________________________________________________________
1830 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1833 // Should be between 0 and 6
1836 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1837 AliInfo("The number of groups must be between 0 and 6!");
1840 fNgroupprf = numberGroupsPRF;
1844 ///////////////////////////////////////////////////////////////////////////////////////////
1845 // Per tracklet: store or reset the info, fill the histos with the info
1846 //////////////////////////////////////////////////////////////////////////////////////////
1847 //_____________________________________________________________________________
1848 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1851 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1852 // Correct from the gain correction before
1855 // time bin of the cluster not corrected
1856 Int_t time = cl->GetPadTime();
1858 //Correct for the gain coefficient used in the database for reconstruction
1859 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1860 Float_t correction = 1.0;
1861 Float_t normalisation = 6.67;
1862 // we divide with gain in AliTRDclusterizer::Transform...
1863 if( correctthegain > 0 ) normalisation /= correctthegain;
1866 // take dd/dl corrected from the angle of the track
1867 correction = dqdl / (normalisation);
1870 // Fill the fAmpTotal with the charge
1872 if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1875 // Fill the fPHPlace and value
1877 fPHPlace[time] = group[1];
1878 fPHValue[time] = correction;
1882 //____________Offine tracking in the AliTRDtracker_____________________________
1883 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1886 // Reset values per tracklet
1889 //Reset good tracklet
1890 fGoodTracklet = kTRUE;
1892 // Reset the fPHValue
1894 //Reset the fPHValue and fPHPlace
1895 for (Int_t k = 0; k < fTimeMax; k++) {
1901 // Reset the fAmpTotal where we put value
1903 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1908 //____________Offine tracking in the AliTRDtracker_____________________________
1909 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1912 // For the offline tracking or mcm tracklets
1913 // This function will be called in the functions UpdateHistogram...
1914 // to fill the info of a track for the relativ gain calibration
1917 Int_t nb = 0; // Nombre de zones traversees
1918 Int_t fd = -1; // Premiere zone non nulle
1919 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1921 if(nbclusters < fNumberClusters) return;
1922 if(nbclusters > fNumberClustersf) return;
1925 // See if the track goes through different zones
1926 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1927 if (fAmpTotal[k] > 0.0) {
1928 totalcharge += fAmpTotal[k];
1941 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1943 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1944 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1947 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1951 if ((fAmpTotal[fd] > 0.0) &&
1952 (fAmpTotal[fd+1] > 0.0)) {
1953 // One of the two very big
1954 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1956 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1957 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1960 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1963 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1965 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1967 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1968 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1971 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1974 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1977 if (fCalibraMode->GetNfragZ(0) > 1) {
1978 if (fAmpTotal[fd] > 0.0) {
1979 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1980 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1981 // One of the two very big
1982 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1984 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1985 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1988 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1991 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1993 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1995 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1996 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1999 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2001 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
2012 //____________Offine tracking in the AliTRDtracker_____________________________
2013 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2016 // For the offline tracking or mcm tracklets
2017 // This function will be called in the functions UpdateHistogram...
2018 // to fill the info of a track for the drift velocity calibration
2021 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2022 Int_t fd1 = -1; // Premiere zone non nulle
2023 Int_t fd2 = -1; // Deuxieme zone non nulle
2024 Int_t k1 = -1; // Debut de la premiere zone
2025 Int_t k2 = -1; // Debut de la seconde zone
2026 Int_t nbclusters = 0; // number of clusters
2030 // See if the track goes through different zones
2031 for (Int_t k = 0; k < fTimeMax; k++) {
2032 if (fPHValue[k] > 0.0) {
2038 if (fPHPlace[k] != fd1) {
2044 if (fPHPlace[k] != fd2) {
2051 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2053 if(nbclusters < fNumberClusters) return;
2054 if(nbclusters > fNumberClustersf) return;
2060 for (Int_t i = 0; i < fTimeMax; i++) {
2062 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2064 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2066 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2069 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2071 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2077 if ((fd1 == fd2+1) ||
2079 // One of the two fast all the think
2080 if (k2 > (k1+fDifference)) {
2081 //we choose to fill the fd1 with all the values
2083 for (Int_t i = 0; i < fTimeMax; i++) {
2085 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2087 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2091 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2093 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2098 if ((k2+fDifference) < fTimeMax) {
2099 //we choose to fill the fd2 with all the values
2101 for (Int_t i = 0; i < fTimeMax; i++) {
2103 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2105 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2109 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2111 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2117 // Two zones voisines sinon rien!
2118 if (fCalibraMode->GetNfragZ(1) > 1) {
2120 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2121 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2122 // One of the two fast all the think
2123 if (k2 > (k1+fDifference)) {
2124 //we choose to fill the fd1 with all the values
2126 for (Int_t i = 0; i < fTimeMax; i++) {
2128 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2130 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2134 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2136 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2141 if ((k2+fDifference) < fTimeMax) {
2142 //we choose to fill the fd2 with all the values
2144 for (Int_t i = 0; i < fTimeMax; i++) {
2146 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2148 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2152 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2154 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2161 // Two zones voisines sinon rien!
2163 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2164 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2165 // One of the two fast all the think
2166 if (k2 > (k1 + fDifference)) {
2167 //we choose to fill the fd1 with all the values
2169 for (Int_t i = 0; i < fTimeMax; i++) {
2171 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2173 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2177 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2179 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2184 if ((k2+fDifference) < fTimeMax) {
2185 //we choose to fill the fd2 with all the values
2187 for (Int_t i = 0; i < fTimeMax; i++) {
2189 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2191 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2195 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2197 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2209 //////////////////////////////////////////////////////////////////////////////////////////
2210 // DAQ process functions
2211 /////////////////////////////////////////////////////////////////////////////////////////
2212 //_____________________________________________________________________
2213 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2216 // Event Processing loop - AliTRDrawStreamBase
2217 // TestBeam 2007 version
2218 // 0 timebin problem
2222 // Same algorithm as TestBeam but different reader
2225 Int_t withInput = 1;
2227 Double_t phvalue[16][144][36];
2228 for(Int_t k = 0; k < 36; k++){
2229 for(Int_t j = 0; j < 16; j++){
2230 for(Int_t c = 0; c < 144; c++){
2231 phvalue[j][c][k] = 0.0;
2236 fDetectorPreviousTrack = -1;
2239 Int_t nbtimebin = 0;
2240 Int_t baseline = 10;
2247 while (rawStream->Next()) {
2249 Int_t idetector = rawStream->GetDet(); // current detector
2250 Int_t imcm = rawStream->GetMCM(); // current MCM
2251 Int_t irob = rawStream->GetROB(); // current ROB
2253 //printf("Detector %d\n",idetector);
2255 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2258 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2262 for(Int_t k = 0; k < 36; k++){
2263 for(Int_t j = 0; j < 16; j++){
2264 for(Int_t c = 0; c < 144; c++){
2265 phvalue[j][c][k] = 0.0;
2271 fDetectorPreviousTrack = idetector;
2272 fMCMPrevious = imcm;
2273 fROBPrevious = irob;
2275 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2276 if(nbtimebin == 0) return 0;
2277 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2278 fTimeMax = nbtimebin;
2280 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2281 fNumberClustersf = fTimeMax;
2282 fNumberClusters = (Int_t)(0.6*fTimeMax);
2285 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2286 Int_t col = rawStream->GetCol();
2287 Int_t row = rawStream->GetRow();
2290 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2293 for(Int_t itime = 0; itime < nbtimebin; itime++){
2294 phvalue[row][col][itime] = signal[itime]-baseline;
2298 // fill the last one
2299 if(fDetectorPreviousTrack != -1){
2302 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2305 for(Int_t k = 0; k < 36; k++){
2306 for(Int_t j = 0; j < 16; j++){
2307 for(Int_t c = 0; c < 144; c++){
2308 phvalue[j][c][k] = 0.0;
2317 while (rawStream->Next()) {
2319 Int_t idetector = rawStream->GetDet(); // current detector
2320 Int_t imcm = rawStream->GetMCM(); // current MCM
2321 Int_t irob = rawStream->GetROB(); // current ROB
2323 //printf("Detector %d\n",idetector);
2325 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2328 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2331 for(Int_t k = 0; k < 36; k++){
2332 for(Int_t j = 0; j < 16; j++){
2333 for(Int_t c = 0; c < 144; c++){
2334 phvalue[j][c][k] = 0.0;
2340 fDetectorPreviousTrack = idetector;
2341 fMCMPrevious = imcm;
2342 fROBPrevious = irob;
2344 //baseline = rawStream->GetCommonAdditive(); // common baseline
2346 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2347 fNumberClustersf = fTimeMax;
2348 fNumberClusters = (Int_t)(0.6*fTimeMax);
2349 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2350 Int_t col = rawStream->GetCol();
2351 Int_t row = rawStream->GetRow();
2354 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2356 for(Int_t itime = 0; itime < fTimeMax; itime++){
2357 phvalue[row][col][itime] = signal[itime]-baseline;
2361 // fill the last one
2362 if(fDetectorPreviousTrack != -1){
2365 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2368 for(Int_t k = 0; k < 36; k++){
2369 for(Int_t j = 0; j < 16; j++){
2370 for(Int_t c = 0; c < 144; c++){
2371 phvalue[j][c][k] = 0.0;
2381 //_____________________________________________________________________
2382 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2385 // Event Processing loop - AliTRDrawStreamBase
2386 // Use old AliTRDmcmtracklet code
2387 // 0 timebin problem
2391 // Algorithm with mcm tracklet
2394 Int_t withInput = 1;
2396 AliTRDmcm mcm = AliTRDmcm(0);
2397 AliTRDtrigParam *param = AliTRDtrigParam::Instance();
2398 rawStream->SetSharedPadReadout(kTRUE);
2400 fDetectorPreviousTrack = -1;
2404 Int_t nbtimebin = 0;
2405 Int_t baseline = 10;
2412 while (rawStream->Next()) {
2414 Int_t idetector = rawStream->GetDet(); // current detector
2415 Int_t imcm = rawStream->GetMCM(); // current MCM
2416 Int_t irob = rawStream->GetROB(); // current ROB
2417 row = rawStream->GetRow();
2420 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2423 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2429 mcm.SetChaId(idetector);
2431 mcm.SetColRange(0,21);
2434 if(fDetectorPreviousTrack == -1){
2437 mcm.SetChaId(idetector);
2439 mcm.SetColRange(0,21);
2443 fDetectorPreviousTrack = idetector;
2444 fMCMPrevious = imcm;
2445 fROBPrevious = irob;
2447 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2448 if(nbtimebin == 0) return 0;
2449 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2450 fTimeMax = nbtimebin;
2451 fNumberClustersf = fTimeMax;
2452 fNumberClusters = (Int_t)(0.6*fTimeMax);
2453 param->SetTimeRange(0,fTimeMax);
2455 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2457 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2458 Int_t adc = rawStream->GetADC();
2461 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2463 for(Int_t itime = 0; itime < nbtimebin; itime++){
2464 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2468 // fill the last one
2469 if(fDetectorPreviousTrack != -1){
2472 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2477 mcm.SetRobId(fROBPrevious);
2478 mcm.SetChaId(fDetectorPreviousTrack);
2480 mcm.SetColRange(0,21);
2487 while (rawStream->Next()) {
2489 Int_t idetector = rawStream->GetDet(); // current detector
2490 Int_t imcm = rawStream->GetMCM(); // current MCM
2491 Int_t irob = rawStream->GetROB(); // current ROB
2492 row = rawStream->GetRow();
2494 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2497 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2503 mcm.SetChaId(idetector);
2505 mcm.SetColRange(0,21);
2509 fDetectorPreviousTrack = idetector;
2510 fMCMPrevious = imcm;
2511 fROBPrevious = irob;
2513 //baseline = rawStream->GetCommonAdditive(); // common baseline
2515 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2516 fNumberClustersf = fTimeMax;
2517 fNumberClusters = (Int_t)(0.6*fTimeMax);
2518 param->SetTimeRange(0,fTimeMax);
2519 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2520 Int_t adc = rawStream->GetADC();
2523 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2525 for(Int_t itime = 0; itime < fTimeMax; itime++){
2526 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2530 // fill the last one
2531 if(fDetectorPreviousTrack != -1){
2534 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2538 mcm.SetRobId(fROBPrevious);
2539 mcm.SetChaId(fDetectorPreviousTrack);
2541 mcm.SetColRange(0,21);
2549 //_____________________________________________________________________
2550 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2553 // Event processing loop - AliRawReader
2554 // Testbeam 2007 version
2557 AliTRDrawStreamBase rawStream(rawReader);
2559 rawReader->Select("TRD");
2561 return ProcessEventDAQ(&rawStream, nocheck);
2564 //_________________________________________________________________________
2565 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2567 eventHeaderStruct *event,
2570 eventHeaderStruct* /*event*/,
2577 // process date event
2578 // Testbeam 2007 version
2581 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2582 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2586 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2591 //_____________________________________________________________________
2592 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliRawReader *rawReader, Bool_t nocheck)
2595 // Event processing loop - AliRawReader
2596 // use the old mcm traklet code
2599 AliTRDrawStreamBase rawStream(rawReader);
2601 rawReader->Select("TRD");
2603 return ProcessEventDAQV1(&rawStream, nocheck);
2605 //_________________________________________________________________________
2606 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(
2608 eventHeaderStruct *event,
2611 eventHeaderStruct* /*event*/,
2618 // process date event
2619 // use the old mcm tracklet code
2622 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2623 Int_t result=ProcessEventDAQV1(rawReader, nocheck);
2627 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2632 //////////////////////////////////////////////////////////////////////////////
2633 // Routine inside the DAQ process
2634 /////////////////////////////////////////////////////////////////////////////
2635 //_______________________________________________________________________
2636 Int_t AliTRDCalibraFillHisto::FillDAQ(AliTRDmcm *mcm){
2639 // Return 2 if some tracklets are found and used, 1 if nothing
2646 for (Int_t iSeed = 0; iSeed < 4; iSeed++) {
2648 if (mcm->GetSeedCol()[iSeed] < 0) {
2652 nbev += TestTracklet(mcm->GetChaId(),mcm->GetRow(),iSeed,mcm);
2657 if(nbev > 0) nbev = 2;
2663 //__________________________________________________________________________
2664 Int_t AliTRDCalibraFillHisto::TestTracklet( Int_t idet, Int_t row, Int_t iSeed, AliTRDmcm *mcm){
2667 // Build the tracklet and return if the tracklet if finally used or not (1/0)
2672 AliTRDmcmTracklet mcmtracklet = AliTRDmcmTracklet();
2673 //mcmtracklet.Reset();
2674 mcmtracklet.SetDetector(idet);
2675 mcmtracklet.SetRow(row);
2676 mcmtracklet.SetN(0);
2678 Int_t iCol, iCol1, iCol2, track[3];
2679 iCol = mcm->GetSeedCol()[iSeed]; // 0....20 (MCM)
2680 mcm->GetColRange(iCol1,iCol2); // range in the pad plane
2683 for (Int_t iTime = 0; iTime < fTimeMax; iTime++) {
2685 amp[0] = mcm->GetADC(iCol-1,iTime);
2686 amp[1] = mcm->GetADC(iCol ,iTime);
2687 amp[2] = mcm->GetADC(iCol+1,iTime);
2689 if(mcm->IsCluster(iCol,iTime)) {
2691 mcmtracklet.AddCluster(iCol+iCol1,iTime,amp,track);
2694 else if ((iCol+1+1) < 21) {
2696 amp[0] = mcm->GetADC(iCol-1+1,iTime);
2697 amp[1] = mcm->GetADC(iCol +1,iTime);
2698 amp[2] = mcm->GetADC(iCol+1+1,iTime);
2700 if(mcm->IsCluster(iCol+1,iTime)) {
2702 mcmtracklet.AddCluster(iCol+1+iCol1,iTime,amp,track);
2711 nbev = UpdateHistogramcm(&mcmtracklet);
2716 //____________Online trackling in AliTRDtrigger________________________________
2717 Int_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2720 // Return if the tracklet is finally used or not (1/0) for calibration
2725 //fGoodTracklet = kTRUE;
2727 // Localisation of the Xbins involved
2728 Int_t idect = trk->GetDetector();
2729 Int_t idectrue = trk->GetDetector();
2732 Int_t nbclusters = trk->GetNclusters();
2734 // Eventuelle correction due to track angle in z direction
2735 Float_t correction = 1.0;
2736 if (fMcmCorrectAngle) {
2737 Float_t z = trk->GetRowz();
2738 Float_t r = trk->GetTime0();
2739 correction = r / TMath::Sqrt((r*r+z*z));
2743 Int_t row = trk->GetRow();
2746 // Boucle sur les clusters
2747 // Condition on number of cluster: don't come from the middle of the detector
2750 for(Int_t k =0; k < 36; k++) amph[k]=0.0;
2751 Double_t ampTotal = 0.0;
2753 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2755 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2756 Int_t time = trk->GetClusterTime(icl);
2757 Int_t col = trk->GetClusterCol(icl);
2759 //CheckGoodTrackletV0(idect,row,col);
2761 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2762 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2763 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2765 ampTotal += (Float_t) (amp[0]+amp[1]+amp[2]);
2766 amph[time]=amp[0]+amp[1]+amp[2];
2768 if(fDebugLevel > 0){
2769 if ( !fDebugStreamer ) {
2771 TDirectory *backup = gDirectory;
2772 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2773 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2776 Double_t amp0 = amp[0];
2777 Double_t amp1 = amp[1];
2778 Double_t amp2 = amp[2];
2780 (* fDebugStreamer) << "UpdateHistogramcm0"<<
2781 "nbclusters="<<nbclusters<<
2788 "detector="<<idectrue<<
2792 } // Boucle clusters
2794 if((amph[0] > 100.0) || (!fGoodTracklet) || (trk->GetNclusters() < fNumberClusters) || (trk->GetNclusters() > fNumberClustersf)) used = 0;
2797 for(Int_t k = 0; k < fTimeMax; k++) UpdateDAQ(idect,0,0,k,amph[k],fTimeMax);
2798 //((TH2I *)GetCH2d()->Fill(ampTotal/30.0,idect));
2802 if(fDebugLevel > 0){
2803 if ( !fDebugStreamer ) {
2805 TDirectory *backup = gDirectory;
2806 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2807 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2810 Double_t amph0 = amph[0];
2811 Double_t amphlast = amph[fTimeMax-1];
2812 Double_t rms = TMath::RMS(fTimeMax,amph);
2813 Int_t goodtracklet = (Int_t) fGoodTracklet;
2815 (* fDebugStreamer) << "UpdateHistogramcm1"<<
2816 "nbclusters="<<nbclusters<<
2817 "ampTotal="<<ampTotal<<
2819 "detector="<<idectrue<<
2821 "amphlast="<<amphlast<<
2822 "goodtracklet="<<goodtracklet<<
2830 //_______________________________________________________________________
2831 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2834 // Look for the maximum by collapsing over the time
2835 // Sum over four pad col and two pad row
2841 Int_t idect = fDetectorPreviousTrack;
2842 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2844 for(Int_t tb = 0; tb < 36; tb++){
2848 //fGoodTracklet = kTRUE;
2849 //fDetectorPreviousTrack = 0;
2852 ///////////////////////////
2854 /////////////////////////
2858 Double_t integralMax = -1;
2860 for (Int_t ir = 1; ir <= 15; ir++)
2862 for (Int_t ic = 2; ic <= 142; ic++)
2864 Double_t integral = 0;
2865 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2867 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2869 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2870 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2873 for(Int_t tb = 0; tb< fTimeMax; tb++){
2874 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2880 if (integralMax < integral)
2884 integralMax = integral;
2885 } // check max integral
2889 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2891 if((imaxRow == 0) || (imaxCol == 0)) {
2895 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2896 //if(!fGoodTracklet) used = 1;;
2898 // /////////////////////////////////////////////////////
2899 // sum ober 2 row and 4 pad cols for each time bins
2900 // ////////////////////////////////////////////////////
2903 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2905 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2907 for(Int_t it = 0; it < fTimeMax; it++){
2908 sum[it] += phvalue[ir][ic][it];
2914 Double_t sumcharge = 0.0;
2915 for(Int_t it = 0; it < fTimeMax; it++){
2916 sumcharge += sum[it];
2917 if(sum[it] > 20.0) nbcl++;
2921 /////////////////////////////////////////////////////////
2923 ////////////////////////////////////////////////////////
2924 if(fDebugLevel > 0){
2925 if ( !fDebugStreamer ) {
2927 TDirectory *backup = gDirectory;
2928 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2929 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2932 Double_t amph0 = sum[0];
2933 Double_t amphlast = sum[fTimeMax-1];
2934 Double_t rms = TMath::RMS(fTimeMax,sum);
2935 Int_t goodtracklet = (Int_t) fGoodTracklet;
2936 for(Int_t it = 0; it < fTimeMax; it++){
2937 Double_t clustera = sum[it];
2939 (* fDebugStreamer) << "FillDAQa"<<
2940 "ampTotal="<<sumcharge<<
2943 "detector="<<idect<<
2945 "amphlast="<<amphlast<<
2946 "goodtracklet="<<goodtracklet<<
2947 "clustera="<<clustera<<
2954 ////////////////////////////////////////////////////////
2956 ///////////////////////////////////////////////////////
2957 if(sum[0] > 100.0) used = 1;
2958 if(nbcl < fNumberClusters) used = 1;
2959 if(nbcl > fNumberClustersf) used = 1;
2961 //if(fDetectorPreviousTrack == 15){
2962 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2964 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2966 for(Int_t it = 0; it < fTimeMax; it++){
2967 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2969 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2974 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2976 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2983 //____________Online trackling in AliTRDtrigger________________________________
2984 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2988 // Fill a simple average pulse height
2992 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2998 //____________Write_____________________________________________________
2999 //_____________________________________________________________________
3000 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
3003 // Write infos to file
3007 if ( fDebugStreamer ) {
3008 delete fDebugStreamer;
3009 fDebugStreamer = 0x0;
3012 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
3017 ,fNumberUsedPh[1]));
3019 TDirectory *backup = gDirectory;
3025 option = "recreate";
3027 TFile f(filename,option.Data());
3029 TStopwatch stopwatch;
3032 f.WriteTObject(fCalibraVector);
3037 f.WriteTObject(fCH2d);
3042 f.WriteTObject(fPH2d);
3047 f.WriteTObject(fPRF2d);
3050 if(fLinearFitterOn){
3051 AnalyseLinearFitter();
3052 f.WriteTObject(fLinearVdriftFit);
3057 if ( backup ) backup->cd();
3059 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
3060 ,stopwatch.RealTime(),stopwatch.CpuTime()));
3062 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3064 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3065 //___________________________________________probe the histos__________________________________________________
3066 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
3069 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
3070 // debug mode with 2 for TH2I and 3 for TProfile2D
3071 // It gives a pointer to a Double_t[7] with the info following...
3072 // [0] : number of calibration groups with entries
3073 // [1] : minimal number of entries found
3074 // [2] : calibration group number of the min
3075 // [3] : maximal number of entries found
3076 // [4] : calibration group number of the max
3077 // [5] : mean number of entries found
3078 // [6] : mean relative error
3081 Double_t *info = new Double_t[7];
3083 // Number of Xbins (detectors or groups of pads)
3084 Int_t nbins = h->GetNbinsY(); //number of calibration groups
3085 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
3088 Double_t nbwe = 0; //number of calibration groups with entries
3089 Double_t minentries = 0; //minimal number of entries found
3090 Double_t maxentries = 0; //maximal number of entries found
3091 Double_t placemin = 0; //calibration group number of the min
3092 Double_t placemax = -1; //calibration group number of the max
3093 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
3094 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
3096 Double_t counter = 0;
3099 TH1F *nbEntries = 0x0;//distribution of the number of entries
3100 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
3101 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
3103 // Beginning of the loop over the calibration groups
3104 for (Int_t idect = 0; idect < nbins; idect++) {
3106 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
3107 projch->SetDirectory(0);
3109 // Number of entries for this calibration group
3110 Double_t nentries = 0.0;
3112 for (Int_t k = 0; k < nxbins; k++) {
3113 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
3117 for (Int_t k = 0; k < nxbins; k++) {
3118 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3119 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3120 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3128 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3129 nbEntries = new TH1F("Number of entries","Number of entries"
3130 ,100,(Int_t)nentries/2,nentries*2);
3131 nbEntries->SetDirectory(0);
3132 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3134 nbEntriesPerGroup->SetDirectory(0);
3135 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3136 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3137 nbEntriesPerSp->SetDirectory(0);
3140 if(nentries > 0) nbEntries->Fill(nentries);
3141 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3142 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3147 if(nentries > maxentries){
3148 maxentries = nentries;
3152 minentries = nentries;
3154 if(nentries < minentries){
3155 minentries = nentries;
3161 meanstats += nentries;
3163 }//calibration groups loop
3165 if(nbwe > 0) meanstats /= nbwe;
3166 if(counter > 0) meanrelativerror /= counter;
3168 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3169 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3170 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3171 AliInfo(Form("The mean number of entries is %f",meanstats));
3172 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3175 info[1] = minentries;
3177 info[3] = maxentries;
3179 info[5] = meanstats;
3180 info[6] = meanrelativerror;
3183 gStyle->SetPalette(1);
3184 gStyle->SetOptStat(1111);
3185 gStyle->SetPadBorderMode(0);
3186 gStyle->SetCanvasColor(10);
3187 gStyle->SetPadLeftMargin(0.13);
3188 gStyle->SetPadRightMargin(0.01);
3189 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3192 nbEntries->Draw("");
3194 nbEntriesPerSp->SetStats(0);
3195 nbEntriesPerSp->Draw("");
3196 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3198 nbEntriesPerGroup->SetStats(0);
3199 nbEntriesPerGroup->Draw("");
3205 //____________________________________________________________________________
3206 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3209 // Return a Int_t[4] with:
3210 // 0 Mean number of entries
3211 // 1 median of number of entries
3212 // 2 rms of number of entries
3213 // 3 number of group with entries
3216 Double_t *stat = new Double_t[4];
3219 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3220 Double_t *weight = new Double_t[nbofgroups];
3221 Int_t *nonul = new Int_t[nbofgroups];
3223 for(Int_t k = 0; k < nbofgroups; k++){
3224 if(fEntriesCH[k] > 0) {
3226 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3229 else weight[k] = 0.0;
3231 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3232 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3233 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3238 //____________________________________________________________________________
3239 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3242 // Return a Int_t[4] with:
3243 // 0 Mean number of entries
3244 // 1 median of number of entries
3245 // 2 rms of number of entries
3246 // 3 number of group with entries
3249 Double_t *stat = new Double_t[4];
3252 Int_t nbofgroups = 540;
3253 Double_t *weight = new Double_t[nbofgroups];
3254 Int_t *nonul = new Int_t[nbofgroups];
3256 for(Int_t k = 0; k < nbofgroups; k++){
3257 if(fEntriesLinearFitter[k] > 0) {
3259 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3262 else weight[k] = 0.0;
3264 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3265 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3266 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3271 //////////////////////////////////////////////////////////////////////////////////////
3273 //////////////////////////////////////////////////////////////////////////////////////
3274 //_____________________________________________________________________________
3275 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3278 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3279 // If fNgroupprf is zero then no binning in tnp
3283 name += fCalibraMode->GetNz(2);
3285 name += fCalibraMode->GetNrphi(2);
3289 if(fNgroupprf != 0){
3291 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3292 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3293 fPRF2d->SetYTitle("Det/pad groups");
3294 fPRF2d->SetXTitle("Position x/W [pad width units]");
3295 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3296 fPRF2d->SetStats(0);
3299 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3300 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3301 fPRF2d->SetYTitle("Det/pad groups");
3302 fPRF2d->SetXTitle("Position x/W [pad width units]");
3303 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3304 fPRF2d->SetStats(0);
3309 //_____________________________________________________________________________
3310 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3313 // Create the 2D histos
3317 name += fCalibraMode->GetNz(1);
3319 name += fCalibraMode->GetNrphi(1);
3321 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3322 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3324 fPH2d->SetYTitle("Det/pad groups");
3325 fPH2d->SetXTitle("time [#mus]");
3326 fPH2d->SetZTitle("<PH> [a.u.]");
3330 //_____________________________________________________________________________
3331 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3334 // Create the 2D histos
3338 name += fCalibraMode->GetNz(0);
3340 name += fCalibraMode->GetNrphi(0);
3342 fCH2d = new TH2I("CH2d",(const Char_t *) name
3343 ,fNumberBinCharge,0,300,nn,0,nn);
3344 fCH2d->SetYTitle("Det/pad groups");
3345 fCH2d->SetXTitle("charge deposit [a.u]");
3346 fCH2d->SetZTitle("counts");
3351 //////////////////////////////////////////////////////////////////////////////////
3352 // Set relative scale
3353 /////////////////////////////////////////////////////////////////////////////////
3354 //_____________________________________________________________________________
3355 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3358 // Set the factor that will divide the deposited charge
3359 // to fit in the histo range [0,300]
3362 if (RelativeScale > 0.0) {
3363 fRelativeScale = RelativeScale;
3366 AliInfo("RelativeScale must be strict positif!");
3370 //////////////////////////////////////////////////////////////////////////////////
3371 // Quick way to fill a histo
3372 //////////////////////////////////////////////////////////////////////////////////
3373 //_____________________________________________________________________
3374 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3377 // FillCH2d: Marian style
3380 //skip simply the value out of range
3381 if((y>=300.0) || (y<0.0)) return;
3383 //Calcul the y place
3384 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3385 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3388 fCH2d->GetArray()[place]++;
3392 //////////////////////////////////////////////////////////////////////////////////
3393 // Geometrical functions
3394 ///////////////////////////////////////////////////////////////////////////////////
3395 //_____________________________________________________________________________
3396 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3399 // Reconstruct the layer number from the detector number
3402 return ((Int_t) (d % 6));
3406 //_____________________________________________________________________________
3407 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3410 // Reconstruct the stack number from the detector number
3412 const Int_t kNlayer = 6;
3414 return ((Int_t) (d % 30) / kNlayer);
3418 //_____________________________________________________________________________
3419 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3422 // Reconstruct the sector number from the detector number
3426 return ((Int_t) (d / fg));
3429 ///////////////////////////////////////////////////////////////////////////////////
3430 // Getter functions for DAQ of the CH2d and the PH2d
3431 //////////////////////////////////////////////////////////////////////////////////
3432 //_____________________________________________________________________
3433 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3436 // return pointer to fPH2d TProfile2D
3437 // create a new TProfile2D if it doesn't exist allready
3443 fTimeMax = nbtimebin;
3444 fSf = samplefrequency;
3450 //_____________________________________________________________________
3451 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3454 // return pointer to fCH2d TH2I
3455 // create a new TH2I if it doesn't exist allready
3464 ////////////////////////////////////////////////////////////////////////////////////////////
3465 // Drift velocity calibration
3466 ///////////////////////////////////////////////////////////////////////////////////////////
3467 //_____________________________________________________________________
3468 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3471 // return pointer to TLinearFitter Calibration
3472 // if force is true create a new TLinearFitter if it doesn't exist allready
3475 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3476 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3479 // if we are forced and TLinearFitter doesn't yet exist create it
3481 // new TLinearFitter
3482 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3483 fLinearFitterArray.AddAt(linearfitter,detector);
3484 return linearfitter;
3487 //____________________________________________________________________________
3488 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3491 // Analyse array of linear fitter because can not be written
3492 // Store two arrays: one with the param the other one with the error param + number of entries
3495 for(Int_t k = 0; k < 540; k++){
3496 TLinearFitter *linearfitter = GetLinearFitter(k);
3497 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3498 TVectorD *par = new TVectorD(2);
3499 TVectorD pare = TVectorD(2);
3500 TVectorD *parE = new TVectorD(3);
3501 linearfitter->Eval();
3502 linearfitter->GetParameters(*par);
3503 linearfitter->GetErrors(pare);
3504 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3505 (*parE)[0] = pare[0]*ppointError;
3506 (*parE)[1] = pare[1]*ppointError;
3507 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3508 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3509 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);