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)
145 ,fScaleWithTPCSignal(kFALSE)
147 ,fThresholdClusterPRF2(15.0)
148 ,fLimitChargeIntegration(kFALSE)
149 ,fFillWithZero(kFALSE)
150 ,fNormalizeNbOfCluster(kFALSE)
153 ,fCutWithVdriftCalib(kFALSE)
154 ,fMinNbTRDtracklets(0)
155 ,fMinTRDMomentum(0.0)
158 ,fSubVersionGainUsed(0)
159 ,fFirstRunGainLocal(0)
160 ,fVersionGainLocalUsed(0)
161 ,fSubVersionGainLocalUsed(0)
163 ,fVersionVdriftUsed(0)
164 ,fSubVersionVdriftUsed(0)
167 ,fSubVersionExBUsed(0)
168 ,fCalibraMode(new AliTRDCalibraMode())
171 ,fDetectorPreviousTrack(-1)
175 ,fNumberClustersf(30)
176 ,fNumberClustersProcent(0.5)
177 ,fThresholdClustersDAQ(120.0)
185 ,fRangeHistoCharge(150)
186 ,fNumberBinCharge(50)
192 ,fGoodTracklet(kTRUE)
193 ,fLinearFitterTracklet(0x0)
195 ,fEntriesLinearFitter(0x0)
200 ,fLinearFitterArray(540)
201 ,fLinearVdriftFit(0x0)
207 // Default constructor
211 // Init some default values
214 fNumberUsedCh[0] = 0;
215 fNumberUsedCh[1] = 0;
216 fNumberUsedPh[0] = 0;
217 fNumberUsedPh[1] = 0;
219 fGeo = new AliTRDgeometry();
220 fCalibDB = AliTRDcalibDB::Instance();
223 //______________________________________________________________________________________
224 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
231 ,fPRF2dOn(c.fPRF2dOn)
232 ,fHisto2d(c.fHisto2d)
233 ,fVector2d(c.fVector2d)
234 ,fLinearFitterOn(c.fLinearFitterOn)
235 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
236 ,fExbAltFitOn(c.fExbAltFitOn)
237 ,fScaleWithTPCSignal(c.fScaleWithTPCSignal)
238 ,fRelativeScale(c.fRelativeScale)
239 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
240 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
241 ,fFillWithZero(c.fFillWithZero)
242 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
243 ,fMaxCluster(c.fMaxCluster)
244 ,fNbMaxCluster(c.fNbMaxCluster)
245 ,fCutWithVdriftCalib(c.fCutWithVdriftCalib)
246 ,fMinNbTRDtracklets(c.fMinNbTRDtracklets)
247 ,fMinTRDMomentum(c.fMinTRDMomentum)
248 ,fFirstRunGain(c.fFirstRunGain)
249 ,fVersionGainUsed(c.fVersionGainUsed)
250 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
251 ,fFirstRunGainLocal(c.fFirstRunGainLocal)
252 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
253 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
254 ,fFirstRunVdrift(c.fFirstRunVdrift)
255 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
256 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
257 ,fFirstRunExB(c.fFirstRunExB)
258 ,fVersionExBUsed(c.fVersionExBUsed)
259 ,fSubVersionExBUsed(c.fSubVersionExBUsed)
262 ,fDebugLevel(c.fDebugLevel)
263 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
264 ,fMCMPrevious(c.fMCMPrevious)
265 ,fROBPrevious(c.fROBPrevious)
266 ,fNumberClusters(c.fNumberClusters)
267 ,fNumberClustersf(c.fNumberClustersf)
268 ,fNumberClustersProcent(c.fNumberClustersProcent)
269 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
270 ,fNumberRowDAQ(c.fNumberRowDAQ)
271 ,fNumberColDAQ(c.fNumberColDAQ)
272 ,fProcent(c.fProcent)
273 ,fDifference(c.fDifference)
274 ,fNumberTrack(c.fNumberTrack)
275 ,fTimeMax(c.fTimeMax)
277 ,fRangeHistoCharge(c.fRangeHistoCharge)
278 ,fNumberBinCharge(c.fNumberBinCharge)
279 ,fNumberBinPRF(c.fNumberBinPRF)
280 ,fNgroupprf(c.fNgroupprf)
284 ,fGoodTracklet(c.fGoodTracklet)
285 ,fLinearFitterTracklet(0x0)
287 ,fEntriesLinearFitter(0x0)
292 ,fLinearFitterArray(540)
293 ,fLinearVdriftFit(0x0)
301 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
302 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
304 fPH2d = new TProfile2D(*c.fPH2d);
305 fPH2d->SetDirectory(0);
308 fPRF2d = new TProfile2D(*c.fPRF2d);
309 fPRF2d->SetDirectory(0);
312 fCH2d = new TH2I(*c.fCH2d);
313 fCH2d->SetDirectory(0);
315 if(c.fLinearVdriftFit){
316 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
319 fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
322 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
323 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
328 fGeo = new AliTRDgeometry();
329 fCalibDB = AliTRDcalibDB::Instance();
331 fNumberUsedCh[0] = 0;
332 fNumberUsedCh[1] = 0;
333 fNumberUsedPh[0] = 0;
334 fNumberUsedPh[1] = 0;
338 //____________________________________________________________________________________
339 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
342 // AliTRDCalibraFillHisto destructor
346 if ( fDebugStreamer ) delete fDebugStreamer;
348 if ( fCalDetGain ) delete fCalDetGain;
349 if ( fCalROCGain ) delete fCalROCGain;
351 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
355 delete [] fEntriesCH;
356 delete [] fEntriesLinearFitter;
359 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
360 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
363 if(fLinearVdriftFit) delete fLinearVdriftFit;
364 if(fExbAltFit) delete fExbAltFit;
370 //_____________________________________________________________________________
371 void AliTRDCalibraFillHisto::Destroy()
382 //_____________________________________________________________________________
383 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
386 // Delete DebugStreamer
389 if ( fDebugStreamer ) delete fDebugStreamer;
390 fDebugStreamer = 0x0;
393 //_____________________________________________________________________________
394 void AliTRDCalibraFillHisto::ClearHistos()
414 //////////////////////////////////////////////////////////////////////////////////
415 // calibration with AliTRDtrackV1: Init, Update
416 //////////////////////////////////////////////////////////////////////////////////
417 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
418 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
421 // Init the histograms and stuff to be filled
426 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
428 AliInfo("Could not get calibDB");
431 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
433 AliInfo("Could not get CommonParam");
438 if(nboftimebin > 0) fTimeMax = nboftimebin;
439 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
440 if(fTimeMax <= 0) fTimeMax = 30;
441 printf("////////////////////////////////////////////\n");
442 printf("Number of time bins in calibration component %d\n",fTimeMax);
443 printf("////////////////////////////////////////////\n");
444 fSf = parCom->GetSamplingFrequency();
445 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
446 else fRelativeScale = 1.18;
447 fNumberClustersf = fTimeMax;
448 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
450 // Init linear fitter
451 if(!fLinearFitterTracklet) {
452 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
453 fLinearFitterTracklet->StoreData(kTRUE);
456 // Calcul Xbins Chambd0, Chamb2
457 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
458 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
459 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
461 // If vector method On initialised all the stuff
463 fCalibraVector = new AliTRDCalibraVector();
464 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
465 fCalibraVector->SetTimeMax(fTimeMax);
466 if(fNgroupprf != 0) {
467 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
468 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
471 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
472 fCalibraVector->SetPRFRange(1.5);
474 for(Int_t k = 0; k < 3; k++){
475 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
476 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
478 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
479 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
480 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
481 fCalibraVector->SetNbGroupPRF(fNgroupprf);
484 // Create the 2D histos corresponding to the pad groupCalibration mode
487 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
488 ,fCalibraMode->GetNz(0)
489 ,fCalibraMode->GetNrphi(0)));
491 // Create the 2D histo
496 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
497 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
501 fEntriesCH = new Int_t[ntotal0];
502 for(Int_t k = 0; k < ntotal0; k++){
509 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
510 ,fCalibraMode->GetNz(1)
511 ,fCalibraMode->GetNrphi(1)));
513 // Create the 2D histo
518 fPHPlace = new Short_t[fTimeMax];
519 for (Int_t k = 0; k < fTimeMax; k++) {
522 fPHValue = new Float_t[fTimeMax];
523 for (Int_t k = 0; k < fTimeMax; k++) {
527 if (fLinearFitterOn) {
528 if(fLinearFitterDebugOn) {
529 fLinearFitterArray.SetName("ArrayLinearFitters");
530 fEntriesLinearFitter = new Int_t[540];
531 for(Int_t k = 0; k < 540; k++){
532 fEntriesLinearFitter[k] = 0;
535 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
536 TString nameee("Ver");
537 nameee += fVersionExBUsed;
539 nameee += fSubVersionExBUsed;
540 nameee += "FirstRun";
541 nameee += fFirstRunExB;
543 fLinearVdriftFit->SetNameCalibUsed(nameee);
546 fExbAltFit = new AliTRDCalibraExbAltFit();
551 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
552 ,fCalibraMode->GetNz(2)
553 ,fCalibraMode->GetNrphi(2)));
554 // Create the 2D histo
556 CreatePRF2d(ntotal2);
563 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
564 Bool_t AliTRDCalibraFillHisto::InitCalDet()
567 // Init the Gain Cal Det
572 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
574 AliError("No gain det calibration entry found");
577 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
579 AliError("No calDet gain found");
585 fCalDetGain->~AliTRDCalDet();
586 new(fCalDetGain) AliTRDCalDet(*(calDet));
587 }else fCalDetGain = new AliTRDCalDet(*(calDet));
592 name += fVersionGainUsed;
594 name += fSubVersionGainUsed;
596 name += fFirstRunGain;
598 name += fCalibraMode->GetNz(0);
600 name += fCalibraMode->GetNrphi(0);
602 fCH2d->SetTitle(name);
605 TString namee("Ver");
606 namee += fVersionVdriftUsed;
608 namee += fSubVersionVdriftUsed;
610 namee += fFirstRunVdrift;
612 namee += fCalibraMode->GetNz(1);
614 namee += fCalibraMode->GetNrphi(1);
616 fPH2d->SetTitle(namee);
618 // title AliTRDCalibraVdriftLinearFit
619 TString nameee("Ver");
620 nameee += fVersionExBUsed;
622 nameee += fSubVersionExBUsed;
623 nameee += "FirstRun";
624 nameee += fFirstRunExB;
628 fLinearVdriftFit->SetNameCalibUsed(nameee);
635 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
636 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
639 // Init the Gain Cal Pad
644 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainLocalUsed,fSubVersionGainLocalUsed);
646 AliError("No gain pad calibration entry found");
649 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
651 AliError("No calPad gain found");
654 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
656 AliError("No calRoc gain found");
661 fCalROCGain->~AliTRDCalROC();
662 new(fCalROCGain) AliTRDCalROC(*(calRoc));
663 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
672 //____________Offline tracking in the AliTRDtracker____________________________
673 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t,const AliESDtrack *esdtrack)
676 // Use AliTRDtrackV1 for the calibration
680 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
681 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
682 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
683 Bool_t newtr = kTRUE; // new track
687 // Cut on the number of TRD tracklets
689 Int_t numberoftrdtracklets = t->GetNumberOfTracklets();
690 if(numberoftrdtracklets < fMinNbTRDtracklets) return kFALSE;
692 Double_t tpcsignal = 1.0;
693 if(esdtrack) tpcsignal = esdtrack->GetTPCsignal()/50.0;
694 if(fScaleWithTPCSignal && tpcsignal <0.00001) return kFALSE;
698 AliInfo("Could not get calibDB");
703 ///////////////////////////
704 // loop over the tracklet
705 ///////////////////////////
706 for(Int_t itr = 0; itr < 6; itr++){
708 if(!(tracklet = t->GetTracklet(itr))) continue;
709 if(!tracklet->IsOK()) continue;
711 ResetfVariablestracklet();
712 Float_t momentum = t->GetMomentum(itr);
713 if(TMath::Abs(momentum) < fMinTRDMomentum) continue;
716 //////////////////////////////////////////
717 // localisation of the tracklet and dqdl
718 //////////////////////////////////////////
719 Int_t layer = tracklet->GetPlane();
721 while(!(cl = tracklet->GetClusters(ic++))) continue;
722 Int_t detector = cl->GetDetector();
723 if (detector != fDetectorPreviousTrack) {
724 // if not a new track
726 // don't use the rest of this track if in the same plane
727 if (layer == GetLayer(fDetectorPreviousTrack)) {
728 //printf("bad tracklet, same layer for detector %d\n",detector);
732 //Localise the detector bin
733 LocalisationDetectorXbins(detector);
735 if(!fIsHLT) InitCalPad(detector);
738 fDetectorPreviousTrack = detector;
742 ////////////////////////////
743 // loop over the clusters
744 ////////////////////////////
745 Double_t chargeQ = 0.0;
746 Int_t nbclusters = 0;
747 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
748 if(!(cl = tracklet->GetClusters(jc))) continue;
751 // Store the info bis of the tracklet
752 Int_t row = cl->GetPadRow();
753 Int_t col = cl->GetPadCol();
754 CheckGoodTrackletV1(cl);
755 Int_t group[2] = {0,0};
756 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
757 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
758 // Add the charge if shared cluster
759 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
761 //Scale with TPC signal or not
762 if(!fScaleWithTPCSignal) {
763 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
764 //printf("Do not scale now\n");
767 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc)/tpcsignal,group,row,col,cls); //tracklet->GetdQdl(jc)
772 ////////////////////////////////////////
773 // Fill the stuffs if a good tracklet
774 ////////////////////////////////////////
777 // drift velocity unables to cut bad tracklets
778 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
780 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
784 if(fCutWithVdriftCalib) {
785 if(pass) FillTheInfoOfTheTrackCH(nbclusters);
787 FillTheInfoOfTheTrackCH(nbclusters);
793 if(fCutWithVdriftCalib) {
794 if(pass) FillTheInfoOfTheTrackPH();
797 FillTheInfoOfTheTrackPH();
801 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
804 /////////////////////////////////////////////////////////
806 ////////////////////////////////////////////////////////
809 if ( !fDebugStreamer ) {
811 TDirectory *backup = gDirectory;
812 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
813 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
816 Int_t stacke = AliTRDgeometry::GetStack(detector);
817 Int_t sme = AliTRDgeometry::GetSector(detector);
818 Int_t layere = AliTRDgeometry::GetLayer(detector);
820 Float_t b[2] = {0.0,0.0};
821 Float_t bCov[3] = {0.0,0.0,0.0};
822 if(esdtrack) esdtrack->GetImpactParameters(b,bCov);
823 if (bCov[0]<=0 || bCov[2]<=0) {
824 bCov[0]=0; bCov[2]=0;
826 Float_t dcaxy = b[0];
828 Int_t tpcnbclusters = 0;
829 if(esdtrack) tpcnbclusters = esdtrack->GetTPCclusters(0);
830 Double_t ttpcsignal = 0.0;
831 if(esdtrack) ttpcsignal = esdtrack->GetTPCsignal();
832 Int_t cutvdriftlinear = 0;
833 if(!pass) cutvdriftlinear = 1;
835 (* fDebugStreamer) << "FillCharge"<<
836 "detector="<<detector<<
842 "nbtpccls="<<tpcnbclusters<<
843 "tpcsignal="<<ttpcsignal<<
844 "cutvdriftlinear="<<cutvdriftlinear<<
846 "nbtrdtracklet="<<numberoftrdtracklets<<
852 } // if a good tracklet
858 ///////////////////////////////////////////////////////////////////////////////////
859 // Routine inside the update with AliTRDtrack
860 ///////////////////////////////////////////////////////////////////////////////////
861 //____________Offine tracking in the AliTRDtracker_____________________________
862 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
865 // Drift velocity calibration:
866 // Fit the clusters with a straight line
867 // From the slope find the drift velocity
870 ////////////////////////////////////////////////
871 //Number of points: if less than 3 return kFALSE
872 /////////////////////////////////////////////////
873 if(nbclusters <= 2) return kFALSE;
878 // results of the linear fit
879 Double_t dydt = 0.0; // dydt tracklet after straight line fit
880 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
881 Double_t pointError = 0.0; // error after straight line fit
882 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
883 Int_t crossrow = 0; // if it crosses a pad row
884 Int_t rowp = -1; // if it crosses a pad row
885 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
886 fLinearFitterTracklet->ClearPoints();
889 ///////////////////////////////////////////
890 // Take the parameters of the track
891 //////////////////////////////////////////
892 // take now the snp, tnp and tgl from the track
893 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
894 Double_t tnp = 0.0; // dy/dx at the end of the chamber
895 if( TMath::Abs(snp) < 1.){
896 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
898 Double_t tgl = tracklet->GetTgl(); // dz/dl
899 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
901 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
902 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
903 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
904 // at the end with correction due to linear fit
905 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
906 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
909 ////////////////////////////
910 // loop over the clusters
911 ////////////////////////////
913 AliTRDcluster *cl = 0x0;
914 //////////////////////////////
915 // Check no shared clusters
916 //////////////////////////////
917 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
918 cl = tracklet->GetClusters(icc);
921 //////////////////////////////////
923 //////////////////////////////////
925 Float_t sigArr[AliTRDfeeParam::GetNcol()];
926 memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
927 Int_t ncl=0, tbf=0, tbl=0;
929 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
930 if(!(cl = tracklet->GetClusters(ic))) continue;
935 Int_t col = cl->GetPadCol();
936 for(int ip=-1, jp=2; jp<5; ip++, jp++){
938 if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
939 sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
942 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
944 Double_t ycluster = cl->GetY();
945 Int_t time = cl->GetPadTime();
946 Double_t timeis = time/fSf;
947 //See if cross two pad rows
948 Int_t row = cl->GetPadRow();
949 if(rowp==-1) rowp = row;
950 if(row != rowp) crossrow = 1;
952 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
958 ////////////////////////////////////
959 // Do the straight line fit now
960 ///////////////////////////////////
962 fLinearFitterTracklet->ClearPoints();
966 fLinearFitterTracklet->Eval();
967 fLinearFitterTracklet->GetParameters(pars);
968 pointError = TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
969 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
971 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
972 fLinearFitterTracklet->ClearPoints();
974 ////////////////////////////////////
975 // Calc the projection of the clusters on the y direction
976 ///////////////////////////////////
978 Float_t signalSum(0.);
979 Float_t mean = 0.0, rms = 0.0;
980 Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
981 Float_t dz = dzdx*(tbl-tbf)/10;
983 for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
984 signalSum+=sigArr[ip];
987 if(signalSum > 0.0) mean/=signalSum;
989 for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
990 rms+=sigArr[ip]*(ip-mean)*(ip-mean);
992 if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
994 rms -= TMath::Abs(dz*tilt);
998 ////////////////////////////////
1000 ///////////////////////////////
1003 if(fDebugLevel > 1){
1004 if ( !fDebugStreamer ) {
1006 TDirectory *backup = gDirectory;
1007 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1008 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1011 float xcoord = tnp-dzdx*tnt;
1012 float pt = tracklet->GetPt();
1013 Int_t layer = GetLayer(fDetectorPreviousTrack);
1015 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1016 //"snpright="<<snpright<<
1018 "nbclusters="<<nbclusters<<
1019 "detector="<<fDetectorPreviousTrack<<
1027 "crossrow="<<crossrow<<
1028 "errorpar="<<errorpar<<
1029 "pointError="<<pointError<<
1041 /////////////////////////
1043 ////////////////////////
1045 if(nbclusters < fNumberClusters) return kFALSE;
1046 if(nbclusters > fNumberClustersf) return kFALSE;
1047 if(pointError >= 0.3) return kFALSE;
1048 if(crossrow == 1) return kTRUE;
1050 ///////////////////////
1052 //////////////////////
1054 if(fLinearFitterOn){
1055 //Add to the linear fitter of the detector
1056 if( TMath::Abs(snp) < 1.){
1057 Double_t x = tnp-dzdx*tnt;
1058 if(fLinearFitterDebugOn) {
1059 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1060 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1062 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1066 fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
1071 //____________Offine tracking in the AliTRDtracker_____________________________
1072 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1075 // PRF width calibration
1076 // Assume a Gaussian shape: determinate the position of the three pad clusters
1077 // Fit with a straight line
1078 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1079 // Fill the PRF as function of angle of the track
1083 //printf("begin\n");
1084 ///////////////////////////////////////////
1085 // Take the parameters of the track
1086 //////////////////////////////////////////
1087 // take now the snp, tnp and tgl from the track
1088 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1089 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1090 if( TMath::Abs(snp) < 1.){
1091 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1093 Double_t tgl = tracklet->GetTgl(); // dz/dl
1094 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1096 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1097 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1098 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1099 // at the end with correction due to linear fit
1100 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1101 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1103 ///////////////////////////////
1104 // Calculate tnp group shift
1105 ///////////////////////////////
1106 Bool_t echec = kFALSE;
1107 Double_t shift = 0.0;
1108 //Calculate the shift in x coresponding to this tnp
1109 if(fNgroupprf != 0.0){
1110 shift = -3.0*(fNgroupprf-1)-1.5;
1111 Double_t limithigh = -0.2*(fNgroupprf-1);
1112 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1114 while(tnp > limithigh){
1120 // do nothing if out of tnp range
1121 //printf("echec %d\n",(Int_t)echec);
1122 if(echec) return kFALSE;
1124 ///////////////////////
1126 //////////////////////
1128 Int_t nb3pc = 0; // number of three pads clusters used for fit
1129 // to see the difference between the fit and the 3 pad clusters position
1130 Double_t padPositions[AliTRDseedV1::kNtb];
1131 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1132 fLinearFitterTracklet->ClearPoints();
1134 //printf("loop clusters \n");
1135 ////////////////////////////
1136 // loop over the clusters
1137 ////////////////////////////
1138 AliTRDcluster *cl = 0x0;
1139 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1140 // reject shared clusters on pad row
1141 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1142 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1145 cl = tracklet->GetClusters(ic);
1147 Double_t time = cl->GetPadTime();
1148 if((time<=7) || (time>=21)) continue;
1149 Short_t *signals = cl->GetSignals();
1150 Float_t xcenter = 0.0;
1151 Bool_t echec1 = kTRUE;
1153 /////////////////////////////////////////////////////////////
1154 // Center 3 balanced: position with the center of the pad
1155 /////////////////////////////////////////////////////////////
1156 if ((((Float_t) signals[3]) > 0.0) &&
1157 (((Float_t) signals[2]) > 0.0) &&
1158 (((Float_t) signals[4]) > 0.0)) {
1160 // Security if the denomiateur is 0
1161 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1162 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1163 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1164 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1165 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1171 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1172 if(echec1) continue;
1174 ////////////////////////////////////////////////////////
1175 //if no echec1: calculate with the position of the pad
1176 // Position of the cluster
1177 // fill the linear fitter
1178 ///////////////////////////////////////////////////////
1179 Double_t padPosition = xcenter + cl->GetPadCol();
1180 padPositions[ic] = padPosition;
1182 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1187 //printf("Fin loop clusters \n");
1188 //////////////////////////////
1189 // fit with a straight line
1190 /////////////////////////////
1192 fLinearFitterTracklet->ClearPoints();
1195 fLinearFitterTracklet->Eval();
1197 fLinearFitterTracklet->GetParameters(line);
1198 Float_t pointError = -1.0;
1199 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1200 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1202 fLinearFitterTracklet->ClearPoints();
1204 //printf("PRF second loop \n");
1205 ////////////////////////////////////////////////
1206 // Fill the PRF: Second loop over clusters
1207 //////////////////////////////////////////////
1208 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1209 // reject shared clusters on pad row
1210 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1211 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1213 cl = tracklet->GetClusters(ic);
1216 Short_t *signals = cl->GetSignals(); // signal
1217 Double_t time = cl->GetPadTime(); // time bin
1218 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1219 Float_t padPos = cl->GetPadCol(); // middle pad
1220 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1221 Float_t ycenter = 0.0; // relative center charge
1222 Float_t ymin = 0.0; // relative left charge
1223 Float_t ymax = 0.0; // relative right charge
1225 ////////////////////////////////////////////////////////////////
1226 // Calculate ycenter, ymin and ymax even for two pad clusters
1227 ////////////////////////////////////////////////////////////////
1228 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1229 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1230 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1231 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1232 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1233 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1236 /////////////////////////
1237 // Calibration group
1238 ////////////////////////
1239 Int_t rowcl = cl->GetPadRow(); // row of cluster
1240 Int_t colcl = cl->GetPadCol(); // col of cluster
1241 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1242 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1243 Float_t xcl = cl->GetY(); // y cluster
1244 Float_t qcl = cl->GetQ(); // charge cluster
1245 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1246 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1247 Double_t xdiff = dpad; // reconstructed position constant
1248 Double_t x = dpad; // reconstructed position moved
1249 Float_t ep = pointError; // error of fit
1250 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1251 Float_t signal3 = (Float_t)signals[3]; // signal
1252 Float_t signal2 = (Float_t)signals[2]; // signal
1253 Float_t signal4 = (Float_t)signals[4]; // signal
1254 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1258 /////////////////////
1260 ////////////////////
1262 if(fDebugLevel > 1){
1263 if ( !fDebugStreamer ) {
1265 TDirectory *backup = gDirectory;
1266 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1267 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1272 Float_t y = ycenter;
1273 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1274 "caligroup="<<caligroup<<
1275 "detector="<<fDetectorPreviousTrack<<
1278 "npoints="<<nbclusters<<
1287 "padPosition="<<padPositions[ic]<<
1288 "padPosTracklet="<<padPosTracklet<<
1293 "signal1="<<signal1<<
1294 "signal2="<<signal2<<
1295 "signal3="<<signal3<<
1296 "signal4="<<signal4<<
1297 "signal5="<<signal5<<
1303 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1304 "caligroup="<<caligroup<<
1305 "detector="<<fDetectorPreviousTrack<<
1308 "npoints="<<nbclusters<<
1317 "padPosition="<<padPositions[ic]<<
1318 "padPosTracklet="<<padPosTracklet<<
1323 "signal1="<<signal1<<
1324 "signal2="<<signal2<<
1325 "signal3="<<signal3<<
1326 "signal4="<<signal4<<
1327 "signal5="<<signal5<<
1333 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1334 "caligroup="<<caligroup<<
1335 "detector="<<fDetectorPreviousTrack<<
1338 "npoints="<<nbclusters<<
1347 "padPosition="<<padPositions[ic]<<
1348 "padPosTracklet="<<padPosTracklet<<
1353 "signal1="<<signal1<<
1354 "signal2="<<signal2<<
1355 "signal3="<<signal3<<
1356 "signal4="<<signal4<<
1357 "signal5="<<signal5<<
1363 /////////////////////
1365 /////////////////////
1366 if(nbclusters < fNumberClusters) continue;
1367 if(nbclusters > fNumberClustersf) continue;
1368 if(nb3pc <= 5) continue;
1369 if((time >= 21) || (time < 7)) continue;
1370 if(TMath::Abs(qcl) < 80) continue;
1371 if( TMath::Abs(snp) > 1.) continue;
1374 ////////////////////////
1376 ///////////////////////
1378 if(TMath::Abs(dpad) < 1.5) {
1379 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1380 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1381 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1383 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1384 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1385 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1387 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1388 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1389 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1394 if(TMath::Abs(dpad) < 1.5) {
1395 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1396 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1398 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1399 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1400 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1402 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1403 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1404 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1407 } // second loop over clusters
1412 ///////////////////////////////////////////////////////////////////////////////////////
1413 // Pad row col stuff: see if masked or not
1414 ///////////////////////////////////////////////////////////////////////////////////////
1415 //_____________________________________________________________________________
1416 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1419 // See if we are not near a masked pad
1422 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1426 //_____________________________________________________________________________
1427 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1430 // See if we are not near a masked pad
1433 if (!IsPadOn(detector, col, row)) {
1434 fGoodTracklet = kFALSE;
1438 if (!IsPadOn(detector, col-1, row)) {
1439 fGoodTracklet = kFALSE;
1444 if (!IsPadOn(detector, col+1, row)) {
1445 fGoodTracklet = kFALSE;
1450 //_____________________________________________________________________________
1451 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1454 // Look in the choosen database if the pad is On.
1455 // If no the track will be "not good"
1458 // Get the parameter object
1459 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1461 AliInfo("Could not get calibDB");
1465 if (!cal->IsChamberGood(detector) ||
1466 cal->IsChamberNoData(detector) ||
1467 cal->IsPadMasked(detector,col,row)) {
1475 ///////////////////////////////////////////////////////////////////////////////////////
1476 // Calibration groups: calculate the number of groups, localise...
1477 ////////////////////////////////////////////////////////////////////////////////////////
1478 //_____________________________________________________________________________
1479 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1482 // Calculate the calibration group number for i
1485 // Row of the cluster and position in the pad groups
1487 if (fCalibraMode->GetNnZ(i) != 0) {
1488 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1492 // Col of the cluster and position in the pad groups
1494 if (fCalibraMode->GetNnRphi(i) != 0) {
1495 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1498 return posc*fCalibraMode->GetNfragZ(i)+posr;
1501 //____________________________________________________________________________________
1502 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1505 // Calculate the total number of calibration groups
1511 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1513 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1518 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1520 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1525 fCalibraMode->ModePadCalibration(2,i);
1526 fCalibraMode->ModePadFragmentation(0,2,0,i);
1527 fCalibraMode->SetDetChamb2(i);
1528 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1529 fCalibraMode->ModePadCalibration(0,i);
1530 fCalibraMode->ModePadFragmentation(0,0,0,i);
1531 fCalibraMode->SetDetChamb0(i);
1532 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1533 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1537 //_____________________________________________________________________________
1538 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1541 // Set the mode of calibration group in the z direction for the parameter i
1546 fCalibraMode->SetNz(i, Nz);
1549 AliInfo("You have to choose between 0 and 4");
1553 //_____________________________________________________________________________
1554 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1557 // Set the mode of calibration group in the rphi direction for the parameter i
1562 fCalibraMode->SetNrphi(i ,Nrphi);
1565 AliInfo("You have to choose between 0 and 6");
1570 //_____________________________________________________________________________
1571 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1574 // Set the mode of calibration group all together
1576 if(fVector2d == kTRUE) {
1577 AliInfo("Can not work with the vector method");
1580 fCalibraMode->SetAllTogether(i);
1584 //_____________________________________________________________________________
1585 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1588 // Set the mode of calibration group per supermodule
1590 if(fVector2d == kTRUE) {
1591 AliInfo("Can not work with the vector method");
1594 fCalibraMode->SetPerSuperModule(i);
1598 //____________Set the pad calibration variables for the detector_______________
1599 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1602 // For the detector calcul the first Xbins and set the number of row
1603 // and col pads per calibration groups, the number of calibration
1604 // groups in the detector.
1607 // first Xbins of the detector
1609 fCalibraMode->CalculXBins(detector,0);
1612 fCalibraMode->CalculXBins(detector,1);
1615 fCalibraMode->CalculXBins(detector,2);
1618 // fragmentation of idect
1619 for (Int_t i = 0; i < 3; i++) {
1620 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1621 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1622 , (Int_t) GetStack(detector)
1623 , (Int_t) GetSector(detector),i);
1629 //_____________________________________________________________________________
1630 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1633 // Should be between 0 and 6
1636 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1637 AliInfo("The number of groups must be between 0 and 6!");
1640 fNgroupprf = numberGroupsPRF;
1644 ///////////////////////////////////////////////////////////////////////////////////////////
1645 // Per tracklet: store or reset the info, fill the histos with the info
1646 //////////////////////////////////////////////////////////////////////////////////////////
1647 //_____________________________________________________________________________
1648 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)
1651 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1652 // Correct from the gain correction before
1653 // cls is shared cluster if any
1654 // Return the charge
1658 //printf("StoreInfoCHPHtrack\n");
1660 // time bin of the cluster not corrected
1661 Int_t time = cl->GetPadTime();
1662 Float_t charge = TMath::Abs(cl->GetQ());
1664 charge += TMath::Abs(cls->GetQ());
1665 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1668 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1670 //Correct for the gain coefficient used in the database for reconstruction
1671 Float_t correctthegain = 1.0;
1672 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1673 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1674 Float_t correction = 1.0;
1675 Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1676 // we divide with gain in AliTRDclusterizer::Transform...
1677 if( correctthegain > 0 ) normalisation /= correctthegain;
1680 // take dd/dl corrected from the angle of the track
1681 correction = dqdl / (normalisation);
1684 // Fill the fAmpTotal with the charge
1686 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1687 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1688 fAmpTotal[(Int_t) group[0]] += correction;
1692 // Fill the fPHPlace and value
1694 fPHPlace[time] = group[1];
1695 fPHValue[time] = charge;
1701 //____________Offine tracking in the AliTRDtracker_____________________________
1702 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1705 // Reset values per tracklet
1708 //Reset good tracklet
1709 fGoodTracklet = kTRUE;
1711 // Reset the fPHValue
1713 //Reset the fPHValue and fPHPlace
1714 for (Int_t k = 0; k < fTimeMax; k++) {
1720 // Reset the fAmpTotal where we put value
1722 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1727 //____________Offine tracking in the AliTRDtracker_____________________________
1728 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1731 // For the offline tracking or mcm tracklets
1732 // This function will be called in the functions UpdateHistogram...
1733 // to fill the info of a track for the relativ gain calibration
1736 Int_t nb = 0; // Nombre de zones traversees
1737 Int_t fd = -1; // Premiere zone non nulle
1738 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1740 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1742 if(nbclusters < fNumberClusters) return;
1743 if(nbclusters > fNumberClustersf) return;
1746 // Normalize with the number of clusters
1747 Double_t normalizeCst = fRelativeScale;
1748 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1750 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1752 // See if the track goes through different zones
1753 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1754 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1755 if (fAmpTotal[k] > 0.0) {
1756 totalcharge += fAmpTotal[k];
1764 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1770 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1772 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1773 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1776 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1780 if ((fAmpTotal[fd] > 0.0) &&
1781 (fAmpTotal[fd+1] > 0.0)) {
1782 // One of the two very big
1783 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1785 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1786 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1789 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1792 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1794 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1796 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1797 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1800 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1803 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1806 if (fCalibraMode->GetNfragZ(0) > 1) {
1807 if (fAmpTotal[fd] > 0.0) {
1808 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1809 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1810 // One of the two very big
1811 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1813 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1814 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1817 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1820 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1822 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1824 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1825 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1828 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1830 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1841 //____________Offine tracking in the AliTRDtracker_____________________________
1842 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1845 // For the offline tracking or mcm tracklets
1846 // This function will be called in the functions UpdateHistogram...
1847 // to fill the info of a track for the drift velocity calibration
1850 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1851 Int_t fd1 = -1; // Premiere zone non nulle
1852 Int_t fd2 = -1; // Deuxieme zone non nulle
1853 Int_t k1 = -1; // Debut de la premiere zone
1854 Int_t k2 = -1; // Debut de la seconde zone
1855 Int_t nbclusters = 0; // number of clusters
1859 // See if the track goes through different zones
1860 for (Int_t k = 0; k < fTimeMax; k++) {
1861 if (fPHValue[k] > 0.0) {
1867 if (fPHPlace[k] != fd1) {
1873 if (fPHPlace[k] != fd2) {
1880 // See if noise before and after
1881 if(fMaxCluster > 0) {
1882 if(fPHValue[0] > fMaxCluster) return;
1883 if(fTimeMax > fNbMaxCluster) {
1884 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1885 if(fPHValue[k] > fMaxCluster) return;
1890 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1892 if(nbclusters < fNumberClusters) return;
1893 if(nbclusters > fNumberClustersf) return;
1899 for (Int_t i = 0; i < fTimeMax; i++) {
1901 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1903 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1905 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1908 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1910 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1916 if ((fd1 == fd2+1) ||
1918 // One of the two fast all the think
1919 if (k2 > (k1+fDifference)) {
1920 //we choose to fill the fd1 with all the values
1922 for (Int_t i = 0; i < fTimeMax; i++) {
1924 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1926 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1930 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1932 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1937 if ((k2+fDifference) < fTimeMax) {
1938 //we choose to fill the fd2 with all the values
1940 for (Int_t i = 0; i < fTimeMax; i++) {
1942 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1944 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1948 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1950 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1956 // Two zones voisines sinon rien!
1957 if (fCalibraMode->GetNfragZ(1) > 1) {
1959 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1960 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1961 // One of the two fast all the think
1962 if (k2 > (k1+fDifference)) {
1963 //we choose to fill the fd1 with all the values
1965 for (Int_t i = 0; i < fTimeMax; i++) {
1967 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1969 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1973 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1975 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1980 if ((k2+fDifference) < fTimeMax) {
1981 //we choose to fill the fd2 with all the values
1983 for (Int_t i = 0; i < fTimeMax; i++) {
1985 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1987 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1991 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1993 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2000 // Two zones voisines sinon rien!
2002 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2003 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2004 // One of the two fast all the think
2005 if (k2 > (k1 + fDifference)) {
2006 //we choose to fill the fd1 with all the values
2008 for (Int_t i = 0; i < fTimeMax; i++) {
2010 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2012 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2016 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2018 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2023 if ((k2+fDifference) < fTimeMax) {
2024 //we choose to fill the fd2 with all the values
2026 for (Int_t i = 0; i < fTimeMax; i++) {
2028 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2030 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2034 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2036 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2048 //////////////////////////////////////////////////////////////////////////////////////////
2049 // DAQ process functions
2050 /////////////////////////////////////////////////////////////////////////////////////////
2051 //_____________________________________________________________________
2052 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2055 // Event Processing loop - AliTRDrawStream
2057 // 0 timebin problem
2060 // Same algorithm as TestBeam but different reader
2063 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2065 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2066 digitsManager->CreateArrays();
2068 Int_t withInput = 1;
2070 Double_t phvalue[16][144][36];
2071 for(Int_t k = 0; k < 36; k++){
2072 for(Int_t j = 0; j < 16; j++){
2073 for(Int_t c = 0; c < 144; c++){
2074 phvalue[j][c][k] = 0.0;
2079 fDetectorPreviousTrack = -1;
2083 Int_t nbtimebin = 0;
2084 Int_t baseline = 10;
2090 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2092 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2093 // printf("there is ADC data on this chamber!\n");
2095 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2096 if (digits->HasData()) { //array
2098 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2099 if (indexes->IsAllocated() == kFALSE) {
2100 AliError("Indexes do not exist!");
2105 indexes->ResetCounters();
2107 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2108 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2109 //while (rawStream->Next()) {
2111 Int_t idetector = det; // current detector
2112 //Int_t imcm = rawStream->GetMCM(); // current MCM
2113 //Int_t irob = rawStream->GetROB(); // current ROB
2116 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2118 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2121 for(Int_t k = 0; k < 36; k++){
2122 for(Int_t j = 0; j < 16; j++){
2123 for(Int_t c = 0; c < 144; c++){
2124 phvalue[j][c][k] = 0.0;
2130 fDetectorPreviousTrack = idetector;
2131 //fMCMPrevious = imcm;
2132 //fROBPrevious = irob;
2134 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2135 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2136 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2137 baseline = digitParam->GetADCbaseline(det); // baseline
2139 if(nbtimebin == 0) return 0;
2140 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2141 fTimeMax = nbtimebin;
2143 fNumberClustersf = fTimeMax;
2144 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2147 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2148 // phvalue[row][col][itime] = signal[itime]-baseline;
2149 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2150 /*if(phvalue[iRow][iCol][itime] >= 20) {
2151 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2155 (Short_t)(digits->GetData(iRow,iCol,itime)),
2162 // fill the last one
2163 if(fDetectorPreviousTrack != -1){
2166 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2167 // printf("\n ---> withinput %d\n\n",withInput);
2169 for(Int_t k = 0; k < 36; k++){
2170 for(Int_t j = 0; j < 16; j++){
2171 for(Int_t c = 0; c < 144; c++){
2172 phvalue[j][c][k] = 0.0;
2180 digitsManager->ClearArrays(det);
2182 delete digitsManager;
2187 //_____________________________________________________________________
2188 //////////////////////////////////////////////////////////////////////////////
2189 // Routine inside the DAQ process
2190 /////////////////////////////////////////////////////////////////////////////
2191 //_______________________________________________________________________
2192 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2195 // Look for the maximum by collapsing over the time
2196 // Sum over four pad col and two pad row
2202 Int_t idect = fDetectorPreviousTrack;
2203 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2205 for(Int_t tb = 0; tb < 36; tb++){
2209 //fGoodTracklet = kTRUE;
2210 //fDetectorPreviousTrack = 0;
2213 ///////////////////////////
2215 /////////////////////////
2219 Double_t integralMax = -1;
2221 for (Int_t ir = 1; ir <= 15; ir++)
2223 for (Int_t ic = 2; ic <= 142; ic++)
2225 Double_t integral = 0;
2226 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2228 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2230 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2231 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2234 for(Int_t tb = 0; tb< fTimeMax; tb++){
2235 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2240 if (integralMax < integral)
2244 integralMax = integral;
2246 } // check max integral
2250 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2251 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2256 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2260 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2261 //if(!fGoodTracklet) used = 1;;
2263 // /////////////////////////////////////////////////////
2264 // sum ober 2 row and 4 pad cols for each time bins
2265 // ////////////////////////////////////////////////////
2269 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2271 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2273 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2274 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2276 for(Int_t it = 0; it < fTimeMax; it++){
2277 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2284 Double_t sumcharge = 0.0;
2285 for(Int_t it = 0; it < fTimeMax; it++){
2286 sumcharge += sum[it];
2287 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2291 /////////////////////////////////////////////////////////
2293 ////////////////////////////////////////////////////////
2294 if(fDebugLevel > 1){
2295 if ( !fDebugStreamer ) {
2297 TDirectory *backup = gDirectory;
2298 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2299 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2302 Double_t amph0 = sum[0];
2303 Double_t amphlast = sum[fTimeMax-1];
2304 Double_t rms = TMath::RMS(fTimeMax,sum);
2305 Int_t goodtracklet = (Int_t) fGoodTracklet;
2306 for(Int_t it = 0; it < fTimeMax; it++){
2307 Double_t clustera = sum[it];
2309 (* fDebugStreamer) << "FillDAQa"<<
2310 "ampTotal="<<sumcharge<<
2313 "detector="<<idect<<
2315 "amphlast="<<amphlast<<
2316 "goodtracklet="<<goodtracklet<<
2317 "clustera="<<clustera<<
2325 ////////////////////////////////////////////////////////
2327 ///////////////////////////////////////////////////////
2328 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2329 if(sum[0] > 100.0) used = 1;
2330 if(nbcl < fNumberClusters) used = 1;
2331 if(nbcl > fNumberClustersf) used = 1;
2333 //if(fDetectorPreviousTrack == 15){
2334 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2336 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2338 for(Int_t it = 0; it < fTimeMax; it++){
2339 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2341 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2343 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2345 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2350 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2352 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2359 //____________Online trackling in AliTRDtrigger________________________________
2360 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2364 // Fill a simple average pulse height
2368 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2374 //____________Write_____________________________________________________
2375 //_____________________________________________________________________
2376 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2379 // Write infos to file
2383 if ( fDebugStreamer ) {
2384 delete fDebugStreamer;
2385 fDebugStreamer = 0x0;
2388 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2393 ,fNumberUsedPh[1]));
2395 TDirectory *backup = gDirectory;
2401 option = "recreate";
2403 TFile f(filename,option.Data());
2405 TStopwatch stopwatch;
2408 f.WriteTObject(fCalibraVector);
2413 f.WriteTObject(fCH2d);
2418 f.WriteTObject(fPH2d);
2423 f.WriteTObject(fPRF2d);
2426 if(fLinearFitterOn){
2427 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2428 f.WriteTObject(fLinearVdriftFit);
2433 if ( backup ) backup->cd();
2435 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2436 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2438 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2440 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2441 //___________________________________________probe the histos__________________________________________________
2442 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2445 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2446 // debug mode with 2 for TH2I and 3 for TProfile2D
2447 // It gives a pointer to a Double_t[7] with the info following...
2448 // [0] : number of calibration groups with entries
2449 // [1] : minimal number of entries found
2450 // [2] : calibration group number of the min
2451 // [3] : maximal number of entries found
2452 // [4] : calibration group number of the max
2453 // [5] : mean number of entries found
2454 // [6] : mean relative error
2457 Double_t *info = new Double_t[7];
2459 // Number of Xbins (detectors or groups of pads)
2460 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2461 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2464 Double_t nbwe = 0; //number of calibration groups with entries
2465 Double_t minentries = 0; //minimal number of entries found
2466 Double_t maxentries = 0; //maximal number of entries found
2467 Double_t placemin = 0; //calibration group number of the min
2468 Double_t placemax = -1; //calibration group number of the max
2469 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2470 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2472 Double_t counter = 0;
2475 TH1F *nbEntries = 0x0;//distribution of the number of entries
2476 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2477 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2479 // Beginning of the loop over the calibration groups
2480 for (Int_t idect = 0; idect < nbins; idect++) {
2482 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2483 projch->SetDirectory(0);
2485 // Number of entries for this calibration group
2486 Double_t nentries = 0.0;
2488 for (Int_t k = 0; k < nxbins; k++) {
2489 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2493 for (Int_t k = 0; k < nxbins; k++) {
2494 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2495 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2496 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2505 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2506 nbEntries->SetDirectory(0);
2507 nbEntries->Fill(nentries);
2508 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2509 nbEntriesPerGroup->SetDirectory(0);
2510 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2511 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));
2512 nbEntriesPerSp->SetDirectory(0);
2513 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2518 if(nentries > maxentries){
2519 maxentries = nentries;
2523 minentries = nentries;
2525 if(nentries < minentries){
2526 minentries = nentries;
2532 meanstats += nentries;
2534 }//calibration groups loop
2536 if(nbwe > 0) meanstats /= nbwe;
2537 if(counter > 0) meanrelativerror /= counter;
2539 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2540 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2541 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2542 AliInfo(Form("The mean number of entries is %f",meanstats));
2543 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2546 info[1] = minentries;
2548 info[3] = maxentries;
2550 info[5] = meanstats;
2551 info[6] = meanrelativerror;
2553 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2554 gStyle->SetPalette(1);
2555 gStyle->SetOptStat(1111);
2556 gStyle->SetPadBorderMode(0);
2557 gStyle->SetCanvasColor(10);
2558 gStyle->SetPadLeftMargin(0.13);
2559 gStyle->SetPadRightMargin(0.01);
2560 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2563 nbEntries->Draw("");
2565 nbEntriesPerSp->SetStats(0);
2566 nbEntriesPerSp->Draw("");
2567 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2569 nbEntriesPerGroup->SetStats(0);
2570 nbEntriesPerGroup->Draw("");
2576 //____________________________________________________________________________
2577 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2580 // Return a Int_t[4] with:
2581 // 0 Mean number of entries
2582 // 1 median of number of entries
2583 // 2 rms of number of entries
2584 // 3 number of group with entries
2587 Double_t *stat = new Double_t[4];
2590 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2592 Double_t *weight = new Double_t[nbofgroups];
2593 Double_t *nonul = new Double_t[nbofgroups];
2595 for(Int_t k = 0; k < nbofgroups; k++){
2596 if(fEntriesCH[k] > 0) {
2598 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2601 else weight[k] = 0.0;
2603 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2604 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2605 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2613 //____________________________________________________________________________
2614 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2617 // Return a Int_t[4] with:
2618 // 0 Mean number of entries
2619 // 1 median of number of entries
2620 // 2 rms of number of entries
2621 // 3 number of group with entries
2624 Double_t *stat = new Double_t[4];
2627 Int_t nbofgroups = 540;
2628 Double_t *weight = new Double_t[nbofgroups];
2629 Int_t *nonul = new Int_t[nbofgroups];
2631 for(Int_t k = 0; k < nbofgroups; k++){
2632 if(fEntriesLinearFitter[k] > 0) {
2634 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2637 else weight[k] = 0.0;
2639 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2640 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2641 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2649 //////////////////////////////////////////////////////////////////////////////////////
2651 //////////////////////////////////////////////////////////////////////////////////////
2652 //_____________________________________________________________________________
2653 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2656 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2657 // If fNgroupprf is zero then no binning in tnp
2661 name += fCalibraMode->GetNz(2);
2663 name += fCalibraMode->GetNrphi(2);
2667 if(fNgroupprf != 0){
2669 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2670 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2671 fPRF2d->SetYTitle("Det/pad groups");
2672 fPRF2d->SetXTitle("Position x/W [pad width units]");
2673 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2674 fPRF2d->SetStats(0);
2677 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2678 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2679 fPRF2d->SetYTitle("Det/pad groups");
2680 fPRF2d->SetXTitle("Position x/W [pad width units]");
2681 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2682 fPRF2d->SetStats(0);
2687 //_____________________________________________________________________________
2688 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2691 // Create the 2D histos
2694 TString name("Ver");
2695 name += fVersionVdriftUsed;
2697 name += fSubVersionVdriftUsed;
2699 name += fFirstRunVdrift;
2701 name += fCalibraMode->GetNz(1);
2703 name += fCalibraMode->GetNrphi(1);
2705 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2706 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2708 fPH2d->SetYTitle("Det/pad groups");
2709 fPH2d->SetXTitle("time [#mus]");
2710 fPH2d->SetZTitle("<PH> [a.u.]");
2714 //_____________________________________________________________________________
2715 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2718 // Create the 2D histos
2721 TString name("Ver");
2722 name += fVersionGainUsed;
2724 name += fSubVersionGainUsed;
2726 name += fFirstRunGain;
2728 name += fCalibraMode->GetNz(0);
2730 name += fCalibraMode->GetNrphi(0);
2732 fCH2d = new TH2I("CH2d",(const Char_t *) name
2733 ,(Int_t)fNumberBinCharge,0,fRangeHistoCharge,nn,0,nn);
2734 fCH2d->SetYTitle("Det/pad groups");
2735 fCH2d->SetXTitle("charge deposit [a.u]");
2736 fCH2d->SetZTitle("counts");
2741 //////////////////////////////////////////////////////////////////////////////////
2742 // Set relative scale
2743 /////////////////////////////////////////////////////////////////////////////////
2744 //_____________________________________________________________________________
2745 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2748 // Set the factor that will divide the deposited charge
2749 // to fit in the histo range [0,fRangeHistoCharge]
2752 if (RelativeScale > 0.0) {
2753 fRelativeScale = RelativeScale;
2756 AliInfo("RelativeScale must be strict positif!");
2760 //////////////////////////////////////////////////////////////////////////////////
2761 // Quick way to fill a histo
2762 //////////////////////////////////////////////////////////////////////////////////
2763 //_____________________________________________________________________
2764 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2767 // FillCH2d: Marian style
2770 //skip simply the value out of range
2771 if((y>=fRangeHistoCharge) || (y<0.0)) return;
2772 if(fRangeHistoCharge < 0.0) return;
2774 //Calcul the y place
2775 Int_t yplace = (Int_t) (fNumberBinCharge*y/fRangeHistoCharge)+1;
2776 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2779 fCH2d->GetArray()[place]++;
2783 //////////////////////////////////////////////////////////////////////////////////
2784 // Geometrical functions
2785 ///////////////////////////////////////////////////////////////////////////////////
2786 //_____________________________________________________________________________
2787 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2790 // Reconstruct the layer number from the detector number
2793 return ((Int_t) (d % 6));
2797 //_____________________________________________________________________________
2798 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2801 // Reconstruct the stack number from the detector number
2803 const Int_t kNlayer = 6;
2805 return ((Int_t) (d % 30) / kNlayer);
2809 //_____________________________________________________________________________
2810 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2813 // Reconstruct the sector number from the detector number
2817 return ((Int_t) (d / fg));
2820 ///////////////////////////////////////////////////////////////////////////////////
2821 // Getter functions for DAQ of the CH2d and the PH2d
2822 //////////////////////////////////////////////////////////////////////////////////
2823 //_____________________________________________________________________
2824 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2827 // return pointer to fPH2d TProfile2D
2828 // create a new TProfile2D if it doesn't exist allready
2834 fTimeMax = nbtimebin;
2835 fSf = samplefrequency;
2841 //_____________________________________________________________________
2842 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2845 // return pointer to fCH2d TH2I
2846 // create a new TH2I if it doesn't exist allready
2855 ////////////////////////////////////////////////////////////////////////////////////////////
2856 // Drift velocity calibration
2857 ///////////////////////////////////////////////////////////////////////////////////////////
2858 //_____________________________________________________________________
2859 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2862 // return pointer to TLinearFitter Calibration
2863 // if force is true create a new TLinearFitter if it doesn't exist allready
2866 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2867 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2870 // if we are forced and TLinearFitter doesn't yet exist create it
2872 // new TLinearFitter
2873 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2874 fLinearFitterArray.AddAt(linearfitter,detector);
2875 return linearfitter;
2878 //____________________________________________________________________________
2879 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2882 // Analyse array of linear fitter because can not be written
2883 // Store two arrays: one with the param the other one with the error param + number of entries
2886 for(Int_t k = 0; k < 540; k++){
2887 TLinearFitter *linearfitter = GetLinearFitter(k);
2888 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2889 TVectorD *par = new TVectorD(2);
2890 TVectorD pare = TVectorD(2);
2891 TVectorD *parE = new TVectorD(3);
2892 if((linearfitter->EvalRobust(0.8)==0)) {
2893 //linearfitter->Eval();
2894 linearfitter->GetParameters(*par);
2895 //linearfitter->GetErrors(pare);
2896 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2897 //(*parE)[0] = pare[0]*ppointError;
2898 //(*parE)[1] = pare[1]*ppointError;
2902 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2903 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2904 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);