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, rbailhache@ikf.uni-frankfurt.de)
33 // J. Book (jbook@ikf.uni-frankfurt.de)
35 //////////////////////////////////////////////////////////////////////////////////////
37 #include <TProfile2D.h>
42 #include <TObjArray.h>
47 #include <TStopwatch.h>
49 #include <TDirectory.h>
50 #include <TTreeStream.h>
52 #include <TLinearFitter.h>
56 #include "AliESDtrack.h"
57 #include "AliTRDCalibraFillHisto.h"
58 #include "AliTRDCalibraMode.h"
59 #include "AliTRDCalibraVector.h"
60 #include "AliTRDCalibraVdriftLinearFit.h"
61 #include "AliTRDCalibraExbAltFit.h"
62 #include "AliTRDcalibDB.h"
63 #include "AliTRDCommonParam.h"
64 #include "AliTRDpadPlane.h"
65 #include "AliTRDcluster.h"
66 #include "AliTRDtrackV1.h"
67 #include "AliRawReader.h"
68 #include "AliRawReaderDate.h"
69 #include "AliTRDgeometry.h"
70 #include "AliTRDfeeParam.h"
71 #include "./Cal/AliTRDCalROC.h"
72 #include "./Cal/AliTRDCalPad.h"
73 #include "./Cal/AliTRDCalDet.h"
75 #include "AliTRDdigitsManager.h"
76 #include "AliTRDdigitsParam.h"
77 #include "AliTRDSignalIndex.h"
78 #include "AliTRDarrayADC.h"
80 #include "AliTRDrawStream.h"
82 #include "AliCDBEntry.h"
83 #include "AliCDBManager.h"
90 ClassImp(AliTRDCalibraFillHisto)
92 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
93 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
95 //_____________singleton implementation_________________________________________________
96 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
99 // Singleton implementation
102 if (fgTerminated != kFALSE) {
106 if (fgInstance == 0) {
107 fgInstance = new AliTRDCalibraFillHisto();
114 //______________________________________________________________________________________
115 void AliTRDCalibraFillHisto::Terminate()
118 // Singleton implementation
119 // Deletes the instance of this class
122 fgTerminated = kTRUE;
124 if (fgInstance != 0) {
131 //______________________________________________________________________________________
132 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
142 ,fLinearFitterOn(kFALSE)
143 ,fLinearFitterDebugOn(kFALSE)
144 ,fExbAltFitOn(kFALSE)
146 ,fThresholdClusterPRF2(15.0)
147 ,fLimitChargeIntegration(kFALSE)
148 ,fFillWithZero(kFALSE)
149 ,fNormalizeNbOfCluster(kFALSE)
152 ,fCutWithVdriftCalib(kFALSE)
153 ,fMinNbTRDtracklets(0)
154 ,fMinTRDMomentum(0.0)
157 ,fSubVersionGainUsed(0)
158 ,fFirstRunGainLocal(0)
159 ,fVersionGainLocalUsed(0)
160 ,fSubVersionGainLocalUsed(0)
162 ,fVersionVdriftUsed(0)
163 ,fSubVersionVdriftUsed(0)
166 ,fSubVersionExBUsed(0)
167 ,fCalibraMode(new AliTRDCalibraMode())
170 ,fDetectorPreviousTrack(-1)
174 ,fNumberClustersf(30)
175 ,fNumberClustersProcent(0.5)
176 ,fThresholdClustersDAQ(120.0)
184 ,fNumberBinCharge(50)
190 ,fGoodTracklet(kTRUE)
191 ,fLinearFitterTracklet(0x0)
193 ,fEntriesLinearFitter(0x0)
198 ,fLinearFitterArray(540)
199 ,fLinearVdriftFit(0x0)
205 // Default constructor
209 // Init some default values
212 fNumberUsedCh[0] = 0;
213 fNumberUsedCh[1] = 0;
214 fNumberUsedPh[0] = 0;
215 fNumberUsedPh[1] = 0;
217 fGeo = new AliTRDgeometry();
218 fCalibDB = AliTRDcalibDB::Instance();
221 //______________________________________________________________________________________
222 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
229 ,fPRF2dOn(c.fPRF2dOn)
230 ,fHisto2d(c.fHisto2d)
231 ,fVector2d(c.fVector2d)
232 ,fLinearFitterOn(c.fLinearFitterOn)
233 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
234 ,fExbAltFitOn(c.fExbAltFitOn)
235 ,fRelativeScale(c.fRelativeScale)
236 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
237 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
238 ,fFillWithZero(c.fFillWithZero)
239 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
240 ,fMaxCluster(c.fMaxCluster)
241 ,fNbMaxCluster(c.fNbMaxCluster)
242 ,fCutWithVdriftCalib(c.fCutWithVdriftCalib)
243 ,fMinNbTRDtracklets(c.fMinNbTRDtracklets)
244 ,fMinTRDMomentum(c.fMinTRDMomentum)
245 ,fFirstRunGain(c.fFirstRunGain)
246 ,fVersionGainUsed(c.fVersionGainUsed)
247 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
248 ,fFirstRunGainLocal(c.fFirstRunGainLocal)
249 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
250 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
251 ,fFirstRunVdrift(c.fFirstRunVdrift)
252 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
253 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
254 ,fFirstRunExB(c.fFirstRunExB)
255 ,fVersionExBUsed(c.fVersionExBUsed)
256 ,fSubVersionExBUsed(c.fSubVersionExBUsed)
259 ,fDebugLevel(c.fDebugLevel)
260 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
261 ,fMCMPrevious(c.fMCMPrevious)
262 ,fROBPrevious(c.fROBPrevious)
263 ,fNumberClusters(c.fNumberClusters)
264 ,fNumberClustersf(c.fNumberClustersf)
265 ,fNumberClustersProcent(c.fNumberClustersProcent)
266 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
267 ,fNumberRowDAQ(c.fNumberRowDAQ)
268 ,fNumberColDAQ(c.fNumberColDAQ)
269 ,fProcent(c.fProcent)
270 ,fDifference(c.fDifference)
271 ,fNumberTrack(c.fNumberTrack)
272 ,fTimeMax(c.fTimeMax)
274 ,fNumberBinCharge(c.fNumberBinCharge)
275 ,fNumberBinPRF(c.fNumberBinPRF)
276 ,fNgroupprf(c.fNgroupprf)
280 ,fGoodTracklet(c.fGoodTracklet)
281 ,fLinearFitterTracklet(0x0)
283 ,fEntriesLinearFitter(0x0)
288 ,fLinearFitterArray(540)
289 ,fLinearVdriftFit(0x0)
297 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
298 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
300 fPH2d = new TProfile2D(*c.fPH2d);
301 fPH2d->SetDirectory(0);
304 fPRF2d = new TProfile2D(*c.fPRF2d);
305 fPRF2d->SetDirectory(0);
308 fCH2d = new TH2I(*c.fCH2d);
309 fCH2d->SetDirectory(0);
311 if(c.fLinearVdriftFit){
312 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
315 fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
318 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
319 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
324 fGeo = new AliTRDgeometry();
325 fCalibDB = AliTRDcalibDB::Instance();
327 fNumberUsedCh[0] = 0;
328 fNumberUsedCh[1] = 0;
329 fNumberUsedPh[0] = 0;
330 fNumberUsedPh[1] = 0;
334 //____________________________________________________________________________________
335 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
338 // AliTRDCalibraFillHisto destructor
342 if ( fDebugStreamer ) delete fDebugStreamer;
344 if ( fCalDetGain ) delete fCalDetGain;
345 if ( fCalROCGain ) delete fCalROCGain;
347 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
351 delete [] fEntriesCH;
352 delete [] fEntriesLinearFitter;
355 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
356 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
359 if(fLinearVdriftFit) delete fLinearVdriftFit;
360 if(fExbAltFit) delete fExbAltFit;
366 //_____________________________________________________________________________
367 void AliTRDCalibraFillHisto::Destroy()
378 //_____________________________________________________________________________
379 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
382 // Delete DebugStreamer
385 if ( fDebugStreamer ) delete fDebugStreamer;
386 fDebugStreamer = 0x0;
389 //_____________________________________________________________________________
390 void AliTRDCalibraFillHisto::ClearHistos()
410 //////////////////////////////////////////////////////////////////////////////////
411 // calibration with AliTRDtrackV1: Init, Update
412 //////////////////////////////////////////////////////////////////////////////////
413 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
414 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
417 // Init the histograms and stuff to be filled
422 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
424 AliInfo("Could not get calibDB");
427 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
429 AliInfo("Could not get CommonParam");
434 if(nboftimebin > 0) fTimeMax = nboftimebin;
435 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
436 if(fTimeMax <= 0) fTimeMax = 30;
437 printf("////////////////////////////////////////////\n");
438 printf("Number of time bins in calibration component %d\n",fTimeMax);
439 printf("////////////////////////////////////////////\n");
440 fSf = parCom->GetSamplingFrequency();
441 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
442 else fRelativeScale = 1.18;
443 fNumberClustersf = fTimeMax;
444 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
446 // Init linear fitter
447 if(!fLinearFitterTracklet) {
448 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
449 fLinearFitterTracklet->StoreData(kTRUE);
452 // Calcul Xbins Chambd0, Chamb2
453 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
454 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
455 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
457 // If vector method On initialised all the stuff
459 fCalibraVector = new AliTRDCalibraVector();
460 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
461 fCalibraVector->SetTimeMax(fTimeMax);
462 if(fNgroupprf != 0) {
463 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
464 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
467 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
468 fCalibraVector->SetPRFRange(1.5);
470 for(Int_t k = 0; k < 3; k++){
471 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
472 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
474 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
475 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
476 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
477 fCalibraVector->SetNbGroupPRF(fNgroupprf);
480 // Create the 2D histos corresponding to the pad groupCalibration mode
483 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
484 ,fCalibraMode->GetNz(0)
485 ,fCalibraMode->GetNrphi(0)));
487 // Create the 2D histo
492 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
493 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
497 fEntriesCH = new Int_t[ntotal0];
498 for(Int_t k = 0; k < ntotal0; k++){
505 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
506 ,fCalibraMode->GetNz(1)
507 ,fCalibraMode->GetNrphi(1)));
509 // Create the 2D histo
514 fPHPlace = new Short_t[fTimeMax];
515 for (Int_t k = 0; k < fTimeMax; k++) {
518 fPHValue = new Float_t[fTimeMax];
519 for (Int_t k = 0; k < fTimeMax; k++) {
523 if (fLinearFitterOn) {
524 if(fLinearFitterDebugOn) {
525 fLinearFitterArray.SetName("ArrayLinearFitters");
526 fEntriesLinearFitter = new Int_t[540];
527 for(Int_t k = 0; k < 540; k++){
528 fEntriesLinearFitter[k] = 0;
531 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
532 TString nameee("Ver");
533 nameee += fVersionExBUsed;
535 nameee += fSubVersionExBUsed;
536 nameee += "FirstRun";
537 nameee += fFirstRunExB;
539 fLinearVdriftFit->SetNameCalibUsed(nameee);
542 fExbAltFit = new AliTRDCalibraExbAltFit();
547 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
548 ,fCalibraMode->GetNz(2)
549 ,fCalibraMode->GetNrphi(2)));
550 // Create the 2D histo
552 CreatePRF2d(ntotal2);
559 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
560 Bool_t AliTRDCalibraFillHisto::InitCalDet()
563 // Init the Gain Cal Det
568 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
570 AliError("No gain det calibration entry found");
573 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
575 AliError("No calDet gain found");
581 fCalDetGain->~AliTRDCalDet();
582 new(fCalDetGain) AliTRDCalDet(*(calDet));
583 }else fCalDetGain = new AliTRDCalDet(*(calDet));
588 name += fVersionGainUsed;
590 name += fSubVersionGainUsed;
592 name += fFirstRunGain;
594 name += fCalibraMode->GetNz(0);
596 name += fCalibraMode->GetNrphi(0);
598 fCH2d->SetTitle(name);
601 TString namee("Ver");
602 namee += fVersionVdriftUsed;
604 namee += fSubVersionVdriftUsed;
606 namee += fFirstRunVdrift;
608 namee += fCalibraMode->GetNz(1);
610 namee += fCalibraMode->GetNrphi(1);
612 fPH2d->SetTitle(namee);
614 // title AliTRDCalibraVdriftLinearFit
615 TString nameee("Ver");
616 nameee += fVersionExBUsed;
618 nameee += fSubVersionExBUsed;
619 nameee += "FirstRun";
620 nameee += fFirstRunExB;
624 fLinearVdriftFit->SetNameCalibUsed(nameee);
631 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
632 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
635 // Init the Gain Cal Pad
640 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainLocalUsed,fSubVersionGainLocalUsed);
642 AliError("No gain pad calibration entry found");
645 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
647 AliError("No calPad gain found");
650 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
652 AliError("No calRoc gain found");
657 fCalROCGain->~AliTRDCalROC();
658 new(fCalROCGain) AliTRDCalROC(*(calRoc));
659 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
668 //____________Offline tracking in the AliTRDtracker____________________________
669 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t,const AliESDtrack *esdtrack)
672 // Use AliTRDtrackV1 for the calibration
676 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
677 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
678 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
679 Bool_t newtr = kTRUE; // new track
683 // Cut on the number of TRD tracklets
685 Int_t numberoftrdtracklets = t->GetNumberOfTracklets();
686 if(numberoftrdtracklets < fMinNbTRDtracklets) return kFALSE;
690 AliInfo("Could not get calibDB");
695 ///////////////////////////
696 // loop over the tracklet
697 ///////////////////////////
698 for(Int_t itr = 0; itr < 6; itr++){
700 if(!(tracklet = t->GetTracklet(itr))) continue;
701 if(!tracklet->IsOK()) continue;
703 ResetfVariablestracklet();
704 Float_t momentum = t->GetMomentum(itr);
705 if(TMath::Abs(momentum) < fMinTRDMomentum) continue;
708 //////////////////////////////////////////
709 // localisation of the tracklet and dqdl
710 //////////////////////////////////////////
711 Int_t layer = tracklet->GetPlane();
713 while(!(cl = tracklet->GetClusters(ic++))) continue;
714 Int_t detector = cl->GetDetector();
715 if (detector != fDetectorPreviousTrack) {
716 // if not a new track
718 // don't use the rest of this track if in the same plane
719 if (layer == GetLayer(fDetectorPreviousTrack)) {
720 //printf("bad tracklet, same layer for detector %d\n",detector);
724 //Localise the detector bin
725 LocalisationDetectorXbins(detector);
727 if(!fIsHLT) InitCalPad(detector);
730 fDetectorPreviousTrack = detector;
734 ////////////////////////////
735 // loop over the clusters
736 ////////////////////////////
737 Double_t chargeQ = 0.0;
738 Int_t nbclusters = 0;
739 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
740 if(!(cl = tracklet->GetClusters(jc))) continue;
743 // Store the info bis of the tracklet
744 Int_t row = cl->GetPadRow();
745 Int_t col = cl->GetPadCol();
746 CheckGoodTrackletV1(cl);
747 Int_t group[2] = {0,0};
748 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
749 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
750 // Add the charge if shared cluster
751 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
753 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
756 ////////////////////////////////////////
757 // Fill the stuffs if a good tracklet
758 ////////////////////////////////////////
761 // drift velocity unables to cut bad tracklets
762 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
764 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
768 if(fCutWithVdriftCalib) {
769 if(pass) FillTheInfoOfTheTrackCH(nbclusters);
771 FillTheInfoOfTheTrackCH(nbclusters);
777 if(fCutWithVdriftCalib) {
778 if(pass) FillTheInfoOfTheTrackPH();
781 FillTheInfoOfTheTrackPH();
785 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
788 /////////////////////////////////////////////////////////
790 ////////////////////////////////////////////////////////
793 if ( !fDebugStreamer ) {
795 TDirectory *backup = gDirectory;
796 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
797 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
800 Int_t stacke = AliTRDgeometry::GetStack(detector);
801 Int_t sme = AliTRDgeometry::GetSector(detector);
802 Int_t layere = AliTRDgeometry::GetLayer(detector);
804 Float_t b[2] = {0.0,0.0};
805 Float_t bCov[3] = {0.0,0.0,0.0};
806 if(esdtrack) esdtrack->GetImpactParameters(b,bCov);
807 if (bCov[0]<=0 || bCov[2]<=0) {
808 bCov[0]=0; bCov[2]=0;
810 Float_t dcaxy = b[0];
812 Int_t tpcnbclusters = 0;
813 if(esdtrack) tpcnbclusters = esdtrack->GetTPCclusters(0);
814 Double_t tpcsignal = 0.0;
815 if(esdtrack) tpcsignal = esdtrack->GetTPCsignal();
816 Int_t cutvdriftlinear = 0;
817 if(!pass) cutvdriftlinear = 1;
819 (* fDebugStreamer) << "FillCharge"<<
820 "detector="<<detector<<
826 "nbtpccls="<<tpcnbclusters<<
827 "tpcsignal="<<tpcsignal<<
828 "cutvdriftlinear="<<cutvdriftlinear<<
830 "nbtrdtracklet="<<numberoftrdtracklets<<
836 } // if a good tracklet
842 ///////////////////////////////////////////////////////////////////////////////////
843 // Routine inside the update with AliTRDtrack
844 ///////////////////////////////////////////////////////////////////////////////////
845 //____________Offine tracking in the AliTRDtracker_____________________________
846 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
849 // Drift velocity calibration:
850 // Fit the clusters with a straight line
851 // From the slope find the drift velocity
854 ////////////////////////////////////////////////
855 //Number of points: if less than 3 return kFALSE
856 /////////////////////////////////////////////////
857 if(nbclusters <= 2) return kFALSE;
862 // results of the linear fit
863 Double_t dydt = 0.0; // dydt tracklet after straight line fit
864 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
865 Double_t pointError = 0.0; // error after straight line fit
866 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
867 Int_t crossrow = 0; // if it crosses a pad row
868 Int_t rowp = -1; // if it crosses a pad row
869 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
870 fLinearFitterTracklet->ClearPoints();
873 ///////////////////////////////////////////
874 // Take the parameters of the track
875 //////////////////////////////////////////
876 // take now the snp, tnp and tgl from the track
877 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
878 Double_t tnp = 0.0; // dy/dx at the end of the chamber
879 if( TMath::Abs(snp) < 1.){
880 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
882 Double_t tgl = tracklet->GetTgl(); // dz/dl
883 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
885 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
886 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
887 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
888 // at the end with correction due to linear fit
889 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
890 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
893 ////////////////////////////
894 // loop over the clusters
895 ////////////////////////////
897 AliTRDcluster *cl = 0x0;
898 //////////////////////////////
899 // Check no shared clusters
900 //////////////////////////////
901 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
902 cl = tracklet->GetClusters(icc);
905 //////////////////////////////////
907 //////////////////////////////////
909 Float_t sigArr[AliTRDfeeParam::GetNcol()];
910 memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
911 Int_t ncl=0, tbf=0, tbl=0;
913 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
914 if(!(cl = tracklet->GetClusters(ic))) continue;
919 Int_t col = cl->GetPadCol();
920 for(int ip=-1, jp=2; jp<5; ip++, jp++){
922 if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
923 sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
926 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
928 Double_t ycluster = cl->GetY();
929 Int_t time = cl->GetPadTime();
930 Double_t timeis = time/fSf;
931 //See if cross two pad rows
932 Int_t row = cl->GetPadRow();
933 if(rowp==-1) rowp = row;
934 if(row != rowp) crossrow = 1;
936 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
942 ////////////////////////////////////
943 // Do the straight line fit now
944 ///////////////////////////////////
946 fLinearFitterTracklet->ClearPoints();
950 fLinearFitterTracklet->Eval();
951 fLinearFitterTracklet->GetParameters(pars);
952 pointError = TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
953 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
955 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
956 fLinearFitterTracklet->ClearPoints();
958 ////////////////////////////////////
959 // Calc the projection of the clusters on the y direction
960 ///////////////////////////////////
962 Float_t signalSum(0.);
963 Float_t mean = 0.0, rms = 0.0;
964 Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
965 Float_t dz = dzdx*(tbl-tbf)/10;
967 for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
968 signalSum+=sigArr[ip];
971 if(signalSum > 0.0) mean/=signalSum;
973 for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
974 rms+=sigArr[ip]*(ip-mean)*(ip-mean);
976 if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
978 rms -= TMath::Abs(dz*tilt);
982 ////////////////////////////////
984 ///////////////////////////////
988 if ( !fDebugStreamer ) {
990 TDirectory *backup = gDirectory;
991 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
992 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
995 float xcoord = tnp-dzdx*tnt;
996 float pt = tracklet->GetPt();
997 Int_t layer = GetLayer(fDetectorPreviousTrack);
999 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1000 //"snpright="<<snpright<<
1002 "nbclusters="<<nbclusters<<
1003 "detector="<<fDetectorPreviousTrack<<
1011 "crossrow="<<crossrow<<
1012 "errorpar="<<errorpar<<
1013 "pointError="<<pointError<<
1025 /////////////////////////
1027 ////////////////////////
1029 if(nbclusters < fNumberClusters) return kFALSE;
1030 if(nbclusters > fNumberClustersf) return kFALSE;
1031 if(pointError >= 0.3) return kFALSE;
1032 if(crossrow == 1) return kTRUE;
1034 ///////////////////////
1036 //////////////////////
1038 if(fLinearFitterOn){
1039 //Add to the linear fitter of the detector
1040 if( TMath::Abs(snp) < 1.){
1041 Double_t x = tnp-dzdx*tnt;
1042 if(fLinearFitterDebugOn) {
1043 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1044 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1046 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1050 fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
1055 //____________Offine tracking in the AliTRDtracker_____________________________
1056 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1059 // PRF width calibration
1060 // Assume a Gaussian shape: determinate the position of the three pad clusters
1061 // Fit with a straight line
1062 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1063 // Fill the PRF as function of angle of the track
1067 //printf("begin\n");
1068 ///////////////////////////////////////////
1069 // Take the parameters of the track
1070 //////////////////////////////////////////
1071 // take now the snp, tnp and tgl from the track
1072 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1073 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1074 if( TMath::Abs(snp) < 1.){
1075 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1077 Double_t tgl = tracklet->GetTgl(); // dz/dl
1078 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1080 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1081 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1082 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1083 // at the end with correction due to linear fit
1084 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1085 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1087 ///////////////////////////////
1088 // Calculate tnp group shift
1089 ///////////////////////////////
1090 Bool_t echec = kFALSE;
1091 Double_t shift = 0.0;
1092 //Calculate the shift in x coresponding to this tnp
1093 if(fNgroupprf != 0.0){
1094 shift = -3.0*(fNgroupprf-1)-1.5;
1095 Double_t limithigh = -0.2*(fNgroupprf-1);
1096 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1098 while(tnp > limithigh){
1104 // do nothing if out of tnp range
1105 //printf("echec %d\n",(Int_t)echec);
1106 if(echec) return kFALSE;
1108 ///////////////////////
1110 //////////////////////
1112 Int_t nb3pc = 0; // number of three pads clusters used for fit
1113 // to see the difference between the fit and the 3 pad clusters position
1114 Double_t padPositions[AliTRDseedV1::kNtb];
1115 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1116 fLinearFitterTracklet->ClearPoints();
1118 //printf("loop clusters \n");
1119 ////////////////////////////
1120 // loop over the clusters
1121 ////////////////////////////
1122 AliTRDcluster *cl = 0x0;
1123 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1124 // reject shared clusters on pad row
1125 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1126 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1129 cl = tracklet->GetClusters(ic);
1131 Double_t time = cl->GetPadTime();
1132 if((time<=7) || (time>=21)) continue;
1133 Short_t *signals = cl->GetSignals();
1134 Float_t xcenter = 0.0;
1135 Bool_t echec1 = kTRUE;
1137 /////////////////////////////////////////////////////////////
1138 // Center 3 balanced: position with the center of the pad
1139 /////////////////////////////////////////////////////////////
1140 if ((((Float_t) signals[3]) > 0.0) &&
1141 (((Float_t) signals[2]) > 0.0) &&
1142 (((Float_t) signals[4]) > 0.0)) {
1144 // Security if the denomiateur is 0
1145 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1146 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1147 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1148 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1149 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1155 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1156 if(echec1) continue;
1158 ////////////////////////////////////////////////////////
1159 //if no echec1: calculate with the position of the pad
1160 // Position of the cluster
1161 // fill the linear fitter
1162 ///////////////////////////////////////////////////////
1163 Double_t padPosition = xcenter + cl->GetPadCol();
1164 padPositions[ic] = padPosition;
1166 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1171 //printf("Fin loop clusters \n");
1172 //////////////////////////////
1173 // fit with a straight line
1174 /////////////////////////////
1176 fLinearFitterTracklet->ClearPoints();
1179 fLinearFitterTracklet->Eval();
1181 fLinearFitterTracklet->GetParameters(line);
1182 Float_t pointError = -1.0;
1183 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1184 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1186 fLinearFitterTracklet->ClearPoints();
1188 //printf("PRF second loop \n");
1189 ////////////////////////////////////////////////
1190 // Fill the PRF: Second loop over clusters
1191 //////////////////////////////////////////////
1192 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1193 // reject shared clusters on pad row
1194 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1195 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1197 cl = tracklet->GetClusters(ic);
1200 Short_t *signals = cl->GetSignals(); // signal
1201 Double_t time = cl->GetPadTime(); // time bin
1202 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1203 Float_t padPos = cl->GetPadCol(); // middle pad
1204 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1205 Float_t ycenter = 0.0; // relative center charge
1206 Float_t ymin = 0.0; // relative left charge
1207 Float_t ymax = 0.0; // relative right charge
1209 ////////////////////////////////////////////////////////////////
1210 // Calculate ycenter, ymin and ymax even for two pad clusters
1211 ////////////////////////////////////////////////////////////////
1212 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1213 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1214 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1215 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1216 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1217 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1220 /////////////////////////
1221 // Calibration group
1222 ////////////////////////
1223 Int_t rowcl = cl->GetPadRow(); // row of cluster
1224 Int_t colcl = cl->GetPadCol(); // col of cluster
1225 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1226 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1227 Float_t xcl = cl->GetY(); // y cluster
1228 Float_t qcl = cl->GetQ(); // charge cluster
1229 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1230 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1231 Double_t xdiff = dpad; // reconstructed position constant
1232 Double_t x = dpad; // reconstructed position moved
1233 Float_t ep = pointError; // error of fit
1234 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1235 Float_t signal3 = (Float_t)signals[3]; // signal
1236 Float_t signal2 = (Float_t)signals[2]; // signal
1237 Float_t signal4 = (Float_t)signals[4]; // signal
1238 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1242 /////////////////////
1244 ////////////////////
1246 if(fDebugLevel > 1){
1247 if ( !fDebugStreamer ) {
1249 TDirectory *backup = gDirectory;
1250 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1251 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1256 Float_t y = ycenter;
1257 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1258 "caligroup="<<caligroup<<
1259 "detector="<<fDetectorPreviousTrack<<
1262 "npoints="<<nbclusters<<
1271 "padPosition="<<padPositions[ic]<<
1272 "padPosTracklet="<<padPosTracklet<<
1277 "signal1="<<signal1<<
1278 "signal2="<<signal2<<
1279 "signal3="<<signal3<<
1280 "signal4="<<signal4<<
1281 "signal5="<<signal5<<
1287 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1288 "caligroup="<<caligroup<<
1289 "detector="<<fDetectorPreviousTrack<<
1292 "npoints="<<nbclusters<<
1301 "padPosition="<<padPositions[ic]<<
1302 "padPosTracklet="<<padPosTracklet<<
1307 "signal1="<<signal1<<
1308 "signal2="<<signal2<<
1309 "signal3="<<signal3<<
1310 "signal4="<<signal4<<
1311 "signal5="<<signal5<<
1317 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1318 "caligroup="<<caligroup<<
1319 "detector="<<fDetectorPreviousTrack<<
1322 "npoints="<<nbclusters<<
1331 "padPosition="<<padPositions[ic]<<
1332 "padPosTracklet="<<padPosTracklet<<
1337 "signal1="<<signal1<<
1338 "signal2="<<signal2<<
1339 "signal3="<<signal3<<
1340 "signal4="<<signal4<<
1341 "signal5="<<signal5<<
1347 /////////////////////
1349 /////////////////////
1350 if(nbclusters < fNumberClusters) continue;
1351 if(nbclusters > fNumberClustersf) continue;
1352 if(nb3pc <= 5) continue;
1353 if((time >= 21) || (time < 7)) continue;
1354 if(TMath::Abs(qcl) < 80) continue;
1355 if( TMath::Abs(snp) > 1.) continue;
1358 ////////////////////////
1360 ///////////////////////
1362 if(TMath::Abs(dpad) < 1.5) {
1363 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1364 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1365 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1367 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1368 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1369 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1371 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1372 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1373 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1378 if(TMath::Abs(dpad) < 1.5) {
1379 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1380 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1382 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1383 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1384 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1386 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1387 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1388 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1391 } // second loop over clusters
1396 ///////////////////////////////////////////////////////////////////////////////////////
1397 // Pad row col stuff: see if masked or not
1398 ///////////////////////////////////////////////////////////////////////////////////////
1399 //_____________________________________________________________________________
1400 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1403 // See if we are not near a masked pad
1406 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1410 //_____________________________________________________________________________
1411 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1414 // See if we are not near a masked pad
1417 if (!IsPadOn(detector, col, row)) {
1418 fGoodTracklet = kFALSE;
1422 if (!IsPadOn(detector, col-1, row)) {
1423 fGoodTracklet = kFALSE;
1428 if (!IsPadOn(detector, col+1, row)) {
1429 fGoodTracklet = kFALSE;
1434 //_____________________________________________________________________________
1435 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1438 // Look in the choosen database if the pad is On.
1439 // If no the track will be "not good"
1442 // Get the parameter object
1443 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1445 AliInfo("Could not get calibDB");
1449 if (!cal->IsChamberGood(detector) ||
1450 cal->IsChamberNoData(detector) ||
1451 cal->IsPadMasked(detector,col,row)) {
1459 ///////////////////////////////////////////////////////////////////////////////////////
1460 // Calibration groups: calculate the number of groups, localise...
1461 ////////////////////////////////////////////////////////////////////////////////////////
1462 //_____________________________________________________________________________
1463 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1466 // Calculate the calibration group number for i
1469 // Row of the cluster and position in the pad groups
1471 if (fCalibraMode->GetNnZ(i) != 0) {
1472 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1476 // Col of the cluster and position in the pad groups
1478 if (fCalibraMode->GetNnRphi(i) != 0) {
1479 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1482 return posc*fCalibraMode->GetNfragZ(i)+posr;
1485 //____________________________________________________________________________________
1486 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1489 // Calculate the total number of calibration groups
1495 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1497 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1502 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1504 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1509 fCalibraMode->ModePadCalibration(2,i);
1510 fCalibraMode->ModePadFragmentation(0,2,0,i);
1511 fCalibraMode->SetDetChamb2(i);
1512 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1513 fCalibraMode->ModePadCalibration(0,i);
1514 fCalibraMode->ModePadFragmentation(0,0,0,i);
1515 fCalibraMode->SetDetChamb0(i);
1516 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1517 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1521 //_____________________________________________________________________________
1522 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1525 // Set the mode of calibration group in the z direction for the parameter i
1530 fCalibraMode->SetNz(i, Nz);
1533 AliInfo("You have to choose between 0 and 4");
1537 //_____________________________________________________________________________
1538 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1541 // Set the mode of calibration group in the rphi direction for the parameter i
1546 fCalibraMode->SetNrphi(i ,Nrphi);
1549 AliInfo("You have to choose between 0 and 6");
1554 //_____________________________________________________________________________
1555 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1558 // Set the mode of calibration group all together
1560 if(fVector2d == kTRUE) {
1561 AliInfo("Can not work with the vector method");
1564 fCalibraMode->SetAllTogether(i);
1568 //_____________________________________________________________________________
1569 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1572 // Set the mode of calibration group per supermodule
1574 if(fVector2d == kTRUE) {
1575 AliInfo("Can not work with the vector method");
1578 fCalibraMode->SetPerSuperModule(i);
1582 //____________Set the pad calibration variables for the detector_______________
1583 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1586 // For the detector calcul the first Xbins and set the number of row
1587 // and col pads per calibration groups, the number of calibration
1588 // groups in the detector.
1591 // first Xbins of the detector
1593 fCalibraMode->CalculXBins(detector,0);
1596 fCalibraMode->CalculXBins(detector,1);
1599 fCalibraMode->CalculXBins(detector,2);
1602 // fragmentation of idect
1603 for (Int_t i = 0; i < 3; i++) {
1604 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1605 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1606 , (Int_t) GetStack(detector)
1607 , (Int_t) GetSector(detector),i);
1613 //_____________________________________________________________________________
1614 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1617 // Should be between 0 and 6
1620 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1621 AliInfo("The number of groups must be between 0 and 6!");
1624 fNgroupprf = numberGroupsPRF;
1628 ///////////////////////////////////////////////////////////////////////////////////////////
1629 // Per tracklet: store or reset the info, fill the histos with the info
1630 //////////////////////////////////////////////////////////////////////////////////////////
1631 //_____________________________________________________________________________
1632 Float_t AliTRDCalibraFillHisto::StoreInfoCHPHtrack(const AliTRDcluster *cl,const Double_t dqdl,const Int_t *group,const Int_t row,const Int_t col,const AliTRDcluster *cls)
1635 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1636 // Correct from the gain correction before
1637 // cls is shared cluster if any
1638 // Return the charge
1642 //printf("StoreInfoCHPHtrack\n");
1644 // time bin of the cluster not corrected
1645 Int_t time = cl->GetPadTime();
1646 Float_t charge = TMath::Abs(cl->GetQ());
1648 charge += TMath::Abs(cls->GetQ());
1649 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1652 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1654 //Correct for the gain coefficient used in the database for reconstruction
1655 Float_t correctthegain = 1.0;
1656 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1657 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1658 Float_t correction = 1.0;
1659 Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1660 // we divide with gain in AliTRDclusterizer::Transform...
1661 if( correctthegain > 0 ) normalisation /= correctthegain;
1664 // take dd/dl corrected from the angle of the track
1665 correction = dqdl / (normalisation);
1668 // Fill the fAmpTotal with the charge
1670 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1671 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1672 fAmpTotal[(Int_t) group[0]] += correction;
1676 // Fill the fPHPlace and value
1678 fPHPlace[time] = group[1];
1679 fPHValue[time] = charge;
1685 //____________Offine tracking in the AliTRDtracker_____________________________
1686 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1689 // Reset values per tracklet
1692 //Reset good tracklet
1693 fGoodTracklet = kTRUE;
1695 // Reset the fPHValue
1697 //Reset the fPHValue and fPHPlace
1698 for (Int_t k = 0; k < fTimeMax; k++) {
1704 // Reset the fAmpTotal where we put value
1706 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1711 //____________Offine tracking in the AliTRDtracker_____________________________
1712 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1715 // For the offline tracking or mcm tracklets
1716 // This function will be called in the functions UpdateHistogram...
1717 // to fill the info of a track for the relativ gain calibration
1720 Int_t nb = 0; // Nombre de zones traversees
1721 Int_t fd = -1; // Premiere zone non nulle
1722 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1724 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1726 if(nbclusters < fNumberClusters) return;
1727 if(nbclusters > fNumberClustersf) return;
1730 // Normalize with the number of clusters
1731 Double_t normalizeCst = fRelativeScale;
1732 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1734 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1736 // See if the track goes through different zones
1737 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1738 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1739 if (fAmpTotal[k] > 0.0) {
1740 totalcharge += fAmpTotal[k];
1748 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1754 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1756 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1757 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1760 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1764 if ((fAmpTotal[fd] > 0.0) &&
1765 (fAmpTotal[fd+1] > 0.0)) {
1766 // One of the two very big
1767 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1769 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1770 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1773 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1776 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1778 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1780 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1781 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1784 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1787 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1790 if (fCalibraMode->GetNfragZ(0) > 1) {
1791 if (fAmpTotal[fd] > 0.0) {
1792 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1793 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1794 // One of the two very big
1795 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1797 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1798 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1801 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1804 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1806 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1808 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1809 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1812 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1814 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1825 //____________Offine tracking in the AliTRDtracker_____________________________
1826 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1829 // For the offline tracking or mcm tracklets
1830 // This function will be called in the functions UpdateHistogram...
1831 // to fill the info of a track for the drift velocity calibration
1834 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1835 Int_t fd1 = -1; // Premiere zone non nulle
1836 Int_t fd2 = -1; // Deuxieme zone non nulle
1837 Int_t k1 = -1; // Debut de la premiere zone
1838 Int_t k2 = -1; // Debut de la seconde zone
1839 Int_t nbclusters = 0; // number of clusters
1843 // See if the track goes through different zones
1844 for (Int_t k = 0; k < fTimeMax; k++) {
1845 if (fPHValue[k] > 0.0) {
1851 if (fPHPlace[k] != fd1) {
1857 if (fPHPlace[k] != fd2) {
1864 // See if noise before and after
1865 if(fMaxCluster > 0) {
1866 if(fPHValue[0] > fMaxCluster) return;
1867 if(fTimeMax > fNbMaxCluster) {
1868 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1869 if(fPHValue[k] > fMaxCluster) return;
1874 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1876 if(nbclusters < fNumberClusters) return;
1877 if(nbclusters > fNumberClustersf) return;
1883 for (Int_t i = 0; i < fTimeMax; i++) {
1885 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1887 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1889 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1892 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1894 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1900 if ((fd1 == fd2+1) ||
1902 // One of the two fast all the think
1903 if (k2 > (k1+fDifference)) {
1904 //we choose to fill the fd1 with all the values
1906 for (Int_t i = 0; i < fTimeMax; i++) {
1908 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1910 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1914 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1916 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1921 if ((k2+fDifference) < fTimeMax) {
1922 //we choose to fill the fd2 with all the values
1924 for (Int_t i = 0; i < fTimeMax; i++) {
1926 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1928 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1932 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1934 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1940 // Two zones voisines sinon rien!
1941 if (fCalibraMode->GetNfragZ(1) > 1) {
1943 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1944 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1945 // One of the two fast all the think
1946 if (k2 > (k1+fDifference)) {
1947 //we choose to fill the fd1 with all the values
1949 for (Int_t i = 0; i < fTimeMax; i++) {
1951 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1953 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1957 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1959 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1964 if ((k2+fDifference) < fTimeMax) {
1965 //we choose to fill the fd2 with all the values
1967 for (Int_t i = 0; i < fTimeMax; i++) {
1969 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1971 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1975 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1977 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1984 // Two zones voisines sinon rien!
1986 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1987 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1988 // One of the two fast all the think
1989 if (k2 > (k1 + fDifference)) {
1990 //we choose to fill the fd1 with all the values
1992 for (Int_t i = 0; i < fTimeMax; i++) {
1994 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1996 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2000 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2002 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2007 if ((k2+fDifference) < fTimeMax) {
2008 //we choose to fill the fd2 with all the values
2010 for (Int_t i = 0; i < fTimeMax; i++) {
2012 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2014 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2018 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2020 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2032 //////////////////////////////////////////////////////////////////////////////////////////
2033 // DAQ process functions
2034 /////////////////////////////////////////////////////////////////////////////////////////
2035 //_____________________________________________________________________
2036 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2039 // Event Processing loop - AliTRDrawStream
2041 // 0 timebin problem
2044 // Same algorithm as TestBeam but different reader
2047 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2049 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2050 digitsManager->CreateArrays();
2052 Int_t withInput = 1;
2054 Double_t phvalue[16][144][36];
2055 for(Int_t k = 0; k < 36; k++){
2056 for(Int_t j = 0; j < 16; j++){
2057 for(Int_t c = 0; c < 144; c++){
2058 phvalue[j][c][k] = 0.0;
2063 fDetectorPreviousTrack = -1;
2067 Int_t nbtimebin = 0;
2068 Int_t baseline = 10;
2074 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2076 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2077 // printf("there is ADC data on this chamber!\n");
2079 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2080 if (digits->HasData()) { //array
2082 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2083 if (indexes->IsAllocated() == kFALSE) {
2084 AliError("Indexes do not exist!");
2089 indexes->ResetCounters();
2091 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2092 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2093 //while (rawStream->Next()) {
2095 Int_t idetector = det; // current detector
2096 //Int_t imcm = rawStream->GetMCM(); // current MCM
2097 //Int_t irob = rawStream->GetROB(); // current ROB
2100 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2102 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2105 for(Int_t k = 0; k < 36; k++){
2106 for(Int_t j = 0; j < 16; j++){
2107 for(Int_t c = 0; c < 144; c++){
2108 phvalue[j][c][k] = 0.0;
2114 fDetectorPreviousTrack = idetector;
2115 //fMCMPrevious = imcm;
2116 //fROBPrevious = irob;
2118 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2119 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2120 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2121 baseline = digitParam->GetADCbaseline(det); // baseline
2123 if(nbtimebin == 0) return 0;
2124 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2125 fTimeMax = nbtimebin;
2127 fNumberClustersf = fTimeMax;
2128 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2131 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2132 // phvalue[row][col][itime] = signal[itime]-baseline;
2133 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2134 /*if(phvalue[iRow][iCol][itime] >= 20) {
2135 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2139 (Short_t)(digits->GetData(iRow,iCol,itime)),
2146 // fill the last one
2147 if(fDetectorPreviousTrack != -1){
2150 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2151 // printf("\n ---> withinput %d\n\n",withInput);
2153 for(Int_t k = 0; k < 36; k++){
2154 for(Int_t j = 0; j < 16; j++){
2155 for(Int_t c = 0; c < 144; c++){
2156 phvalue[j][c][k] = 0.0;
2164 digitsManager->ClearArrays(det);
2166 delete digitsManager;
2171 //_____________________________________________________________________
2172 //////////////////////////////////////////////////////////////////////////////
2173 // Routine inside the DAQ process
2174 /////////////////////////////////////////////////////////////////////////////
2175 //_______________________________________________________________________
2176 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2179 // Look for the maximum by collapsing over the time
2180 // Sum over four pad col and two pad row
2186 Int_t idect = fDetectorPreviousTrack;
2187 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2189 for(Int_t tb = 0; tb < 36; tb++){
2193 //fGoodTracklet = kTRUE;
2194 //fDetectorPreviousTrack = 0;
2197 ///////////////////////////
2199 /////////////////////////
2203 Double_t integralMax = -1;
2205 for (Int_t ir = 1; ir <= 15; ir++)
2207 for (Int_t ic = 2; ic <= 142; ic++)
2209 Double_t integral = 0;
2210 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2212 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2214 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2215 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2218 for(Int_t tb = 0; tb< fTimeMax; tb++){
2219 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2224 if (integralMax < integral)
2228 integralMax = integral;
2230 } // check max integral
2234 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2235 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2240 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2244 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2245 //if(!fGoodTracklet) used = 1;;
2247 // /////////////////////////////////////////////////////
2248 // sum ober 2 row and 4 pad cols for each time bins
2249 // ////////////////////////////////////////////////////
2253 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2255 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2257 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2258 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2260 for(Int_t it = 0; it < fTimeMax; it++){
2261 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2268 Double_t sumcharge = 0.0;
2269 for(Int_t it = 0; it < fTimeMax; it++){
2270 sumcharge += sum[it];
2271 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2275 /////////////////////////////////////////////////////////
2277 ////////////////////////////////////////////////////////
2278 if(fDebugLevel > 1){
2279 if ( !fDebugStreamer ) {
2281 TDirectory *backup = gDirectory;
2282 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2283 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2286 Double_t amph0 = sum[0];
2287 Double_t amphlast = sum[fTimeMax-1];
2288 Double_t rms = TMath::RMS(fTimeMax,sum);
2289 Int_t goodtracklet = (Int_t) fGoodTracklet;
2290 for(Int_t it = 0; it < fTimeMax; it++){
2291 Double_t clustera = sum[it];
2293 (* fDebugStreamer) << "FillDAQa"<<
2294 "ampTotal="<<sumcharge<<
2297 "detector="<<idect<<
2299 "amphlast="<<amphlast<<
2300 "goodtracklet="<<goodtracklet<<
2301 "clustera="<<clustera<<
2309 ////////////////////////////////////////////////////////
2311 ///////////////////////////////////////////////////////
2312 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2313 if(sum[0] > 100.0) used = 1;
2314 if(nbcl < fNumberClusters) used = 1;
2315 if(nbcl > fNumberClustersf) used = 1;
2317 //if(fDetectorPreviousTrack == 15){
2318 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2320 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2322 for(Int_t it = 0; it < fTimeMax; it++){
2323 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2325 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2327 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2329 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2334 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2336 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2343 //____________Online trackling in AliTRDtrigger________________________________
2344 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2348 // Fill a simple average pulse height
2352 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2358 //____________Write_____________________________________________________
2359 //_____________________________________________________________________
2360 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2363 // Write infos to file
2367 if ( fDebugStreamer ) {
2368 delete fDebugStreamer;
2369 fDebugStreamer = 0x0;
2372 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2377 ,fNumberUsedPh[1]));
2379 TDirectory *backup = gDirectory;
2385 option = "recreate";
2387 TFile f(filename,option.Data());
2389 TStopwatch stopwatch;
2392 f.WriteTObject(fCalibraVector);
2397 f.WriteTObject(fCH2d);
2402 f.WriteTObject(fPH2d);
2407 f.WriteTObject(fPRF2d);
2410 if(fLinearFitterOn){
2411 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2412 f.WriteTObject(fLinearVdriftFit);
2417 if ( backup ) backup->cd();
2419 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2420 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2422 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2424 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2425 //___________________________________________probe the histos__________________________________________________
2426 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2429 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2430 // debug mode with 2 for TH2I and 3 for TProfile2D
2431 // It gives a pointer to a Double_t[7] with the info following...
2432 // [0] : number of calibration groups with entries
2433 // [1] : minimal number of entries found
2434 // [2] : calibration group number of the min
2435 // [3] : maximal number of entries found
2436 // [4] : calibration group number of the max
2437 // [5] : mean number of entries found
2438 // [6] : mean relative error
2441 Double_t *info = new Double_t[7];
2443 // Number of Xbins (detectors or groups of pads)
2444 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2445 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2448 Double_t nbwe = 0; //number of calibration groups with entries
2449 Double_t minentries = 0; //minimal number of entries found
2450 Double_t maxentries = 0; //maximal number of entries found
2451 Double_t placemin = 0; //calibration group number of the min
2452 Double_t placemax = -1; //calibration group number of the max
2453 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2454 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2456 Double_t counter = 0;
2459 TH1F *nbEntries = 0x0;//distribution of the number of entries
2460 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2461 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2463 // Beginning of the loop over the calibration groups
2464 for (Int_t idect = 0; idect < nbins; idect++) {
2466 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2467 projch->SetDirectory(0);
2469 // Number of entries for this calibration group
2470 Double_t nentries = 0.0;
2472 for (Int_t k = 0; k < nxbins; k++) {
2473 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2477 for (Int_t k = 0; k < nxbins; k++) {
2478 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2479 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2480 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2489 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2490 nbEntries->SetDirectory(0);
2491 nbEntries->Fill(nentries);
2492 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2493 nbEntriesPerGroup->SetDirectory(0);
2494 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2495 if(!((Bool_t)nbEntriesPerSp)) nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule",(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2496 nbEntriesPerSp->SetDirectory(0);
2497 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2502 if(nentries > maxentries){
2503 maxentries = nentries;
2507 minentries = nentries;
2509 if(nentries < minentries){
2510 minentries = nentries;
2516 meanstats += nentries;
2518 }//calibration groups loop
2520 if(nbwe > 0) meanstats /= nbwe;
2521 if(counter > 0) meanrelativerror /= counter;
2523 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2524 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2525 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2526 AliInfo(Form("The mean number of entries is %f",meanstats));
2527 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2530 info[1] = minentries;
2532 info[3] = maxentries;
2534 info[5] = meanstats;
2535 info[6] = meanrelativerror;
2537 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2538 gStyle->SetPalette(1);
2539 gStyle->SetOptStat(1111);
2540 gStyle->SetPadBorderMode(0);
2541 gStyle->SetCanvasColor(10);
2542 gStyle->SetPadLeftMargin(0.13);
2543 gStyle->SetPadRightMargin(0.01);
2544 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2547 nbEntries->Draw("");
2549 nbEntriesPerSp->SetStats(0);
2550 nbEntriesPerSp->Draw("");
2551 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2553 nbEntriesPerGroup->SetStats(0);
2554 nbEntriesPerGroup->Draw("");
2560 //____________________________________________________________________________
2561 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2564 // Return a Int_t[4] with:
2565 // 0 Mean number of entries
2566 // 1 median of number of entries
2567 // 2 rms of number of entries
2568 // 3 number of group with entries
2571 Double_t *stat = new Double_t[4];
2574 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2576 Double_t *weight = new Double_t[nbofgroups];
2577 Double_t *nonul = new Double_t[nbofgroups];
2579 for(Int_t k = 0; k < nbofgroups; k++){
2580 if(fEntriesCH[k] > 0) {
2582 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2585 else weight[k] = 0.0;
2587 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2588 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2589 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2597 //____________________________________________________________________________
2598 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2601 // Return a Int_t[4] with:
2602 // 0 Mean number of entries
2603 // 1 median of number of entries
2604 // 2 rms of number of entries
2605 // 3 number of group with entries
2608 Double_t *stat = new Double_t[4];
2611 Int_t nbofgroups = 540;
2612 Double_t *weight = new Double_t[nbofgroups];
2613 Int_t *nonul = new Int_t[nbofgroups];
2615 for(Int_t k = 0; k < nbofgroups; k++){
2616 if(fEntriesLinearFitter[k] > 0) {
2618 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2621 else weight[k] = 0.0;
2623 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2624 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2625 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2633 //////////////////////////////////////////////////////////////////////////////////////
2635 //////////////////////////////////////////////////////////////////////////////////////
2636 //_____________________________________________________________________________
2637 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2640 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2641 // If fNgroupprf is zero then no binning in tnp
2645 name += fCalibraMode->GetNz(2);
2647 name += fCalibraMode->GetNrphi(2);
2651 if(fNgroupprf != 0){
2653 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2654 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2655 fPRF2d->SetYTitle("Det/pad groups");
2656 fPRF2d->SetXTitle("Position x/W [pad width units]");
2657 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2658 fPRF2d->SetStats(0);
2661 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2662 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2663 fPRF2d->SetYTitle("Det/pad groups");
2664 fPRF2d->SetXTitle("Position x/W [pad width units]");
2665 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2666 fPRF2d->SetStats(0);
2671 //_____________________________________________________________________________
2672 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2675 // Create the 2D histos
2678 TString name("Ver");
2679 name += fVersionVdriftUsed;
2681 name += fSubVersionVdriftUsed;
2683 name += fFirstRunVdrift;
2685 name += fCalibraMode->GetNz(1);
2687 name += fCalibraMode->GetNrphi(1);
2689 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2690 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2692 fPH2d->SetYTitle("Det/pad groups");
2693 fPH2d->SetXTitle("time [#mus]");
2694 fPH2d->SetZTitle("<PH> [a.u.]");
2698 //_____________________________________________________________________________
2699 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2702 // Create the 2D histos
2705 TString name("Ver");
2706 name += fVersionGainUsed;
2708 name += fSubVersionGainUsed;
2710 name += fFirstRunGain;
2712 name += fCalibraMode->GetNz(0);
2714 name += fCalibraMode->GetNrphi(0);
2716 fCH2d = new TH2I("CH2d",(const Char_t *) name
2717 ,fNumberBinCharge,0,300,nn,0,nn);
2718 fCH2d->SetYTitle("Det/pad groups");
2719 fCH2d->SetXTitle("charge deposit [a.u]");
2720 fCH2d->SetZTitle("counts");
2725 //////////////////////////////////////////////////////////////////////////////////
2726 // Set relative scale
2727 /////////////////////////////////////////////////////////////////////////////////
2728 //_____________________________________________________________________________
2729 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2732 // Set the factor that will divide the deposited charge
2733 // to fit in the histo range [0,300]
2736 if (RelativeScale > 0.0) {
2737 fRelativeScale = RelativeScale;
2740 AliInfo("RelativeScale must be strict positif!");
2744 //////////////////////////////////////////////////////////////////////////////////
2745 // Quick way to fill a histo
2746 //////////////////////////////////////////////////////////////////////////////////
2747 //_____________________________________________________________________
2748 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2751 // FillCH2d: Marian style
2754 //skip simply the value out of range
2755 if((y>=300.0) || (y<0.0)) return;
2757 //Calcul the y place
2758 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2759 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2762 fCH2d->GetArray()[place]++;
2766 //////////////////////////////////////////////////////////////////////////////////
2767 // Geometrical functions
2768 ///////////////////////////////////////////////////////////////////////////////////
2769 //_____________________________________________________________________________
2770 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2773 // Reconstruct the layer number from the detector number
2776 return ((Int_t) (d % 6));
2780 //_____________________________________________________________________________
2781 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2784 // Reconstruct the stack number from the detector number
2786 const Int_t kNlayer = 6;
2788 return ((Int_t) (d % 30) / kNlayer);
2792 //_____________________________________________________________________________
2793 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2796 // Reconstruct the sector number from the detector number
2800 return ((Int_t) (d / fg));
2803 ///////////////////////////////////////////////////////////////////////////////////
2804 // Getter functions for DAQ of the CH2d and the PH2d
2805 //////////////////////////////////////////////////////////////////////////////////
2806 //_____________________________________________________________________
2807 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2810 // return pointer to fPH2d TProfile2D
2811 // create a new TProfile2D if it doesn't exist allready
2817 fTimeMax = nbtimebin;
2818 fSf = samplefrequency;
2824 //_____________________________________________________________________
2825 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2828 // return pointer to fCH2d TH2I
2829 // create a new TH2I if it doesn't exist allready
2838 ////////////////////////////////////////////////////////////////////////////////////////////
2839 // Drift velocity calibration
2840 ///////////////////////////////////////////////////////////////////////////////////////////
2841 //_____________________________________________________________________
2842 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2845 // return pointer to TLinearFitter Calibration
2846 // if force is true create a new TLinearFitter if it doesn't exist allready
2849 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2850 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2853 // if we are forced and TLinearFitter doesn't yet exist create it
2855 // new TLinearFitter
2856 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2857 fLinearFitterArray.AddAt(linearfitter,detector);
2858 return linearfitter;
2861 //____________________________________________________________________________
2862 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2865 // Analyse array of linear fitter because can not be written
2866 // Store two arrays: one with the param the other one with the error param + number of entries
2869 for(Int_t k = 0; k < 540; k++){
2870 TLinearFitter *linearfitter = GetLinearFitter(k);
2871 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2872 TVectorD *par = new TVectorD(2);
2873 TVectorD pare = TVectorD(2);
2874 TVectorD *parE = new TVectorD(3);
2875 if((linearfitter->EvalRobust(0.8)==0)) {
2876 //linearfitter->Eval();
2877 linearfitter->GetParameters(*par);
2878 //linearfitter->GetErrors(pare);
2879 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2880 //(*parE)[0] = pare[0]*ppointError;
2881 //(*parE)[1] = pare[1]*ppointError;
2885 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2886 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2887 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);