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 "AliTRDCalibraFillHisto.h"
57 #include "AliTRDCalibraMode.h"
58 #include "AliTRDCalibraVector.h"
59 #include "AliTRDCalibraVdriftLinearFit.h"
60 #include "AliTRDCalibraExbAltFit.h"
61 #include "AliTRDcalibDB.h"
62 #include "AliTRDCommonParam.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcluster.h"
65 #include "AliTRDtrackV1.h"
66 #include "AliRawReader.h"
67 #include "AliRawReaderDate.h"
68 #include "AliTRDgeometry.h"
69 #include "AliTRDfeeParam.h"
70 #include "./Cal/AliTRDCalROC.h"
71 #include "./Cal/AliTRDCalPad.h"
72 #include "./Cal/AliTRDCalDet.h"
74 #include "AliTRDdigitsManager.h"
75 #include "AliTRDdigitsParam.h"
76 #include "AliTRDSignalIndex.h"
77 #include "AliTRDarrayADC.h"
79 #include "AliTRDrawStream.h"
81 #include "AliCDBEntry.h"
82 #include "AliCDBManager.h"
89 ClassImp(AliTRDCalibraFillHisto)
91 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
92 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
94 //_____________singleton implementation_________________________________________________
95 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
98 // Singleton implementation
101 if (fgTerminated != kFALSE) {
105 if (fgInstance == 0) {
106 fgInstance = new AliTRDCalibraFillHisto();
113 //______________________________________________________________________________________
114 void AliTRDCalibraFillHisto::Terminate()
117 // Singleton implementation
118 // Deletes the instance of this class
121 fgTerminated = kTRUE;
123 if (fgInstance != 0) {
130 //______________________________________________________________________________________
131 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
141 ,fLinearFitterOn(kFALSE)
142 ,fLinearFitterDebugOn(kFALSE)
143 ,fExbAltFitOn(kFALSE)
145 ,fThresholdClusterPRF2(15.0)
146 ,fLimitChargeIntegration(kFALSE)
147 ,fFillWithZero(kFALSE)
148 ,fNormalizeNbOfCluster(kFALSE)
153 ,fSubVersionGainUsed(0)
154 ,fFirstRunGainLocal(0)
155 ,fVersionGainLocalUsed(0)
156 ,fSubVersionGainLocalUsed(0)
158 ,fVersionVdriftUsed(0)
159 ,fSubVersionVdriftUsed(0)
162 ,fSubVersionExBUsed(0)
163 ,fCalibraMode(new AliTRDCalibraMode())
166 ,fDetectorPreviousTrack(-1)
170 ,fNumberClustersf(30)
171 ,fNumberClustersProcent(0.5)
172 ,fThresholdClustersDAQ(120.0)
180 ,fNumberBinCharge(50)
186 ,fGoodTracklet(kTRUE)
187 ,fLinearFitterTracklet(0x0)
189 ,fEntriesLinearFitter(0x0)
194 ,fLinearFitterArray(540)
195 ,fLinearVdriftFit(0x0)
201 // Default constructor
205 // Init some default values
208 fNumberUsedCh[0] = 0;
209 fNumberUsedCh[1] = 0;
210 fNumberUsedPh[0] = 0;
211 fNumberUsedPh[1] = 0;
213 fGeo = new AliTRDgeometry();
214 fCalibDB = AliTRDcalibDB::Instance();
217 //______________________________________________________________________________________
218 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
225 ,fPRF2dOn(c.fPRF2dOn)
226 ,fHisto2d(c.fHisto2d)
227 ,fVector2d(c.fVector2d)
228 ,fLinearFitterOn(c.fLinearFitterOn)
229 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
230 ,fExbAltFitOn(c.fExbAltFitOn)
231 ,fRelativeScale(c.fRelativeScale)
232 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
233 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
234 ,fFillWithZero(c.fFillWithZero)
235 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
236 ,fMaxCluster(c.fMaxCluster)
237 ,fNbMaxCluster(c.fNbMaxCluster)
238 ,fFirstRunGain(c.fFirstRunGain)
239 ,fVersionGainUsed(c.fVersionGainUsed)
240 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
241 ,fFirstRunGainLocal(c.fFirstRunGainLocal)
242 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
243 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
244 ,fFirstRunVdrift(c.fFirstRunVdrift)
245 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
246 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
247 ,fFirstRunExB(c.fFirstRunExB)
248 ,fVersionExBUsed(c.fVersionExBUsed)
249 ,fSubVersionExBUsed(c.fSubVersionExBUsed)
252 ,fDebugLevel(c.fDebugLevel)
253 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
254 ,fMCMPrevious(c.fMCMPrevious)
255 ,fROBPrevious(c.fROBPrevious)
256 ,fNumberClusters(c.fNumberClusters)
257 ,fNumberClustersf(c.fNumberClustersf)
258 ,fNumberClustersProcent(c.fNumberClustersProcent)
259 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
260 ,fNumberRowDAQ(c.fNumberRowDAQ)
261 ,fNumberColDAQ(c.fNumberColDAQ)
262 ,fProcent(c.fProcent)
263 ,fDifference(c.fDifference)
264 ,fNumberTrack(c.fNumberTrack)
265 ,fTimeMax(c.fTimeMax)
267 ,fNumberBinCharge(c.fNumberBinCharge)
268 ,fNumberBinPRF(c.fNumberBinPRF)
269 ,fNgroupprf(c.fNgroupprf)
273 ,fGoodTracklet(c.fGoodTracklet)
274 ,fLinearFitterTracklet(0x0)
276 ,fEntriesLinearFitter(0x0)
281 ,fLinearFitterArray(540)
282 ,fLinearVdriftFit(0x0)
290 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
291 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
293 fPH2d = new TProfile2D(*c.fPH2d);
294 fPH2d->SetDirectory(0);
297 fPRF2d = new TProfile2D(*c.fPRF2d);
298 fPRF2d->SetDirectory(0);
301 fCH2d = new TH2I(*c.fCH2d);
302 fCH2d->SetDirectory(0);
304 if(c.fLinearVdriftFit){
305 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
308 fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
311 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
312 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
317 fGeo = new AliTRDgeometry();
318 fCalibDB = AliTRDcalibDB::Instance();
320 fNumberUsedCh[0] = 0;
321 fNumberUsedCh[1] = 0;
322 fNumberUsedPh[0] = 0;
323 fNumberUsedPh[1] = 0;
327 //____________________________________________________________________________________
328 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
331 // AliTRDCalibraFillHisto destructor
335 if ( fDebugStreamer ) delete fDebugStreamer;
337 if ( fCalDetGain ) delete fCalDetGain;
338 if ( fCalROCGain ) delete fCalROCGain;
340 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
344 delete [] fEntriesCH;
345 delete [] fEntriesLinearFitter;
348 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
349 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
352 if(fLinearVdriftFit) delete fLinearVdriftFit;
353 if(fExbAltFit) delete fExbAltFit;
359 //_____________________________________________________________________________
360 void AliTRDCalibraFillHisto::Destroy()
371 //_____________________________________________________________________________
372 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
375 // Delete DebugStreamer
378 if ( fDebugStreamer ) delete fDebugStreamer;
379 fDebugStreamer = 0x0;
382 //_____________________________________________________________________________
383 void AliTRDCalibraFillHisto::ClearHistos()
403 //////////////////////////////////////////////////////////////////////////////////
404 // calibration with AliTRDtrackV1: Init, Update
405 //////////////////////////////////////////////////////////////////////////////////
406 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
407 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
410 // Init the histograms and stuff to be filled
415 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
417 AliInfo("Could not get calibDB");
420 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
422 AliInfo("Could not get CommonParam");
427 if(nboftimebin > 0) fTimeMax = nboftimebin;
428 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
429 if(fTimeMax <= 0) fTimeMax = 30;
430 printf("////////////////////////////////////////////\n");
431 printf("Number of time bins in calibration component %d\n",fTimeMax);
432 printf("////////////////////////////////////////////\n");
433 fSf = parCom->GetSamplingFrequency();
434 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
435 else fRelativeScale = 1.18;
436 fNumberClustersf = fTimeMax;
437 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
439 // Init linear fitter
440 if(!fLinearFitterTracklet) {
441 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
442 fLinearFitterTracklet->StoreData(kTRUE);
445 // Calcul Xbins Chambd0, Chamb2
446 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
447 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
448 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
450 // If vector method On initialised all the stuff
452 fCalibraVector = new AliTRDCalibraVector();
453 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
454 fCalibraVector->SetTimeMax(fTimeMax);
455 if(fNgroupprf != 0) {
456 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
457 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
460 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
461 fCalibraVector->SetPRFRange(1.5);
463 for(Int_t k = 0; k < 3; k++){
464 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
465 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
467 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
468 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
469 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
470 fCalibraVector->SetNbGroupPRF(fNgroupprf);
473 // Create the 2D histos corresponding to the pad groupCalibration mode
476 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
477 ,fCalibraMode->GetNz(0)
478 ,fCalibraMode->GetNrphi(0)));
480 // Create the 2D histo
485 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
486 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
490 fEntriesCH = new Int_t[ntotal0];
491 for(Int_t k = 0; k < ntotal0; k++){
498 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
499 ,fCalibraMode->GetNz(1)
500 ,fCalibraMode->GetNrphi(1)));
502 // Create the 2D histo
507 fPHPlace = new Short_t[fTimeMax];
508 for (Int_t k = 0; k < fTimeMax; k++) {
511 fPHValue = new Float_t[fTimeMax];
512 for (Int_t k = 0; k < fTimeMax; k++) {
516 if (fLinearFitterOn) {
517 if(fLinearFitterDebugOn) {
518 fLinearFitterArray.SetName("ArrayLinearFitters");
519 fEntriesLinearFitter = new Int_t[540];
520 for(Int_t k = 0; k < 540; k++){
521 fEntriesLinearFitter[k] = 0;
524 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
525 TString nameee("Ver");
526 nameee += fVersionExBUsed;
528 nameee += fSubVersionExBUsed;
529 nameee += "FirstRun";
530 nameee += fFirstRunExB;
532 fLinearVdriftFit->SetNameCalibUsed(nameee);
535 fExbAltFit = new AliTRDCalibraExbAltFit();
540 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
541 ,fCalibraMode->GetNz(2)
542 ,fCalibraMode->GetNrphi(2)));
543 // Create the 2D histo
545 CreatePRF2d(ntotal2);
552 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
553 Bool_t AliTRDCalibraFillHisto::InitCalDet()
556 // Init the Gain Cal Det
561 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
563 AliError("No gain det calibration entry found");
566 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
568 AliError("No calDet gain found");
574 fCalDetGain->~AliTRDCalDet();
575 new(fCalDetGain) AliTRDCalDet(*(calDet));
576 }else fCalDetGain = new AliTRDCalDet(*(calDet));
581 name += fVersionGainUsed;
583 name += fSubVersionGainUsed;
585 name += fFirstRunGain;
587 name += fCalibraMode->GetNz(0);
589 name += fCalibraMode->GetNrphi(0);
591 fCH2d->SetTitle(name);
594 TString namee("Ver");
595 namee += fVersionVdriftUsed;
597 namee += fSubVersionVdriftUsed;
599 namee += fFirstRunVdrift;
601 namee += fCalibraMode->GetNz(1);
603 namee += fCalibraMode->GetNrphi(1);
605 fPH2d->SetTitle(namee);
607 // title AliTRDCalibraVdriftLinearFit
608 TString nameee("Ver");
609 nameee += fVersionExBUsed;
611 nameee += fSubVersionExBUsed;
612 nameee += "FirstRun";
613 nameee += fFirstRunExB;
617 fLinearVdriftFit->SetNameCalibUsed(nameee);
624 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
625 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
628 // Init the Gain Cal Pad
633 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainLocalUsed,fSubVersionGainLocalUsed);
635 AliError("No gain pad calibration entry found");
638 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
640 AliError("No calPad gain found");
643 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
645 AliError("No calRoc gain found");
650 fCalROCGain->~AliTRDCalROC();
651 new(fCalROCGain) AliTRDCalROC(*(calRoc));
652 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
661 //____________Offline tracking in the AliTRDtracker____________________________
662 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t)
665 // Use AliTRDtrackV1 for the calibration
669 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
670 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
671 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
672 Bool_t newtr = kTRUE; // new track
675 // AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
678 AliInfo("Could not get calibDB");
683 AliInfo("Could not get calibDB");
688 ///////////////////////////
689 // loop over the tracklet
690 ///////////////////////////
691 for(Int_t itr = 0; itr < 6; itr++){
693 if(!(tracklet = t->GetTracklet(itr))) continue;
694 if(!tracklet->IsOK()) continue;
696 ResetfVariablestracklet();
698 //////////////////////////////////////////
699 // localisation of the tracklet and dqdl
700 //////////////////////////////////////////
701 Int_t layer = tracklet->GetPlane();
703 while(!(cl = tracklet->GetClusters(ic++))) continue;
704 Int_t detector = cl->GetDetector();
705 if (detector != fDetectorPreviousTrack) {
706 // if not a new track
708 // don't use the rest of this track if in the same plane
709 if (layer == GetLayer(fDetectorPreviousTrack)) {
710 //printf("bad tracklet, same layer for detector %d\n",detector);
714 //Localise the detector bin
715 LocalisationDetectorXbins(detector);
717 if(!fIsHLT) InitCalPad(detector);
720 fDetectorPreviousTrack = detector;
724 ////////////////////////////
725 // loop over the clusters
726 ////////////////////////////
727 Int_t nbclusters = 0;
728 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
729 if(!(cl = tracklet->GetClusters(jc))) continue;
732 // Store the info bis of the tracklet
733 Int_t row = cl->GetPadRow();
734 Int_t col = cl->GetPadCol();
735 CheckGoodTrackletV1(cl);
736 Int_t group[2] = {0,0};
737 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
738 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
739 // Add the charge if shared cluster
740 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
742 StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
745 ////////////////////////////////////////
746 // Fill the stuffs if a good tracklet
747 ////////////////////////////////////////
750 // drift velocity unables to cut bad tracklets
751 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
753 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
757 FillTheInfoOfTheTrackCH(nbclusters);
762 FillTheInfoOfTheTrackPH();
765 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
767 } // if a good tracklet
773 ///////////////////////////////////////////////////////////////////////////////////
774 // Routine inside the update with AliTRDtrack
775 ///////////////////////////////////////////////////////////////////////////////////
776 //____________Offine tracking in the AliTRDtracker_____________________________
777 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
780 // Drift velocity calibration:
781 // Fit the clusters with a straight line
782 // From the slope find the drift velocity
785 ////////////////////////////////////////////////
786 //Number of points: if less than 3 return kFALSE
787 /////////////////////////////////////////////////
788 if(nbclusters <= 2) return kFALSE;
793 // results of the linear fit
794 Double_t dydt = 0.0; // dydt tracklet after straight line fit
795 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
796 Double_t pointError = 0.0; // error after straight line fit
797 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
798 Int_t crossrow = 0; // if it crosses a pad row
799 Int_t rowp = -1; // if it crosses a pad row
800 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
801 fLinearFitterTracklet->ClearPoints();
804 ///////////////////////////////////////////
805 // Take the parameters of the track
806 //////////////////////////////////////////
807 // take now the snp, tnp and tgl from the track
808 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
809 Double_t tnp = 0.0; // dy/dx at the end of the chamber
810 if( TMath::Abs(snp) < 1.){
811 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
813 Double_t tgl = tracklet->GetTgl(); // dz/dl
814 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
816 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
817 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
818 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
819 // at the end with correction due to linear fit
820 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
821 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
824 ////////////////////////////
825 // loop over the clusters
826 ////////////////////////////
828 AliTRDcluster *cl = 0x0;
829 //////////////////////////////
830 // Check no shared clusters
831 //////////////////////////////
832 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
833 cl = tracklet->GetClusters(icc);
836 //////////////////////////////////
838 //////////////////////////////////
840 Float_t sigArr[AliTRDfeeParam::GetNcol()];
841 memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
842 Int_t ncl=0, tbf=0, tbl=0;
844 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
845 if(!(cl = tracklet->GetClusters(ic))) continue;
850 Int_t col = cl->GetPadCol();
851 for(int ip=-1, jp=2; jp<5; ip++, jp++){
853 if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
854 sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
857 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
859 Double_t ycluster = cl->GetY();
860 Int_t time = cl->GetPadTime();
861 Double_t timeis = time/fSf;
862 //See if cross two pad rows
863 Int_t row = cl->GetPadRow();
864 if(rowp==-1) rowp = row;
865 if(row != rowp) crossrow = 1;
867 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
873 ////////////////////////////////////
874 // Do the straight line fit now
875 ///////////////////////////////////
877 fLinearFitterTracklet->ClearPoints();
881 fLinearFitterTracklet->Eval();
882 fLinearFitterTracklet->GetParameters(pars);
883 pointError = TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
884 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
886 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
887 fLinearFitterTracklet->ClearPoints();
889 ////////////////////////////////////
890 // Calc the projection of the clusters on the y direction
891 ///////////////////////////////////
893 Float_t signalSum(0.);
894 Float_t mean = 0.0, rms = 0.0;
895 Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
896 Float_t dz = dzdx*(tbl-tbf)/10;
898 for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
899 signalSum+=sigArr[ip];
902 if(signalSum > 0.0) mean/=signalSum;
904 for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
905 rms+=sigArr[ip]*(ip-mean)*(ip-mean);
907 if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
909 rms -= TMath::Abs(dz*tilt);
913 ////////////////////////////////
915 ///////////////////////////////
918 //if(fDebugLevel > 0){
919 if ( !fDebugStreamer ) {
921 TDirectory *backup = gDirectory;
922 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
923 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
926 float xcoord = tnp-dzdx*tnt;
927 float pt = tracklet->GetPt();
928 Int_t layer = GetLayer(fDetectorPreviousTrack);
930 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
931 //"snpright="<<snpright<<
933 "nbclusters="<<nbclusters<<
934 "detector="<<fDetectorPreviousTrack<<
942 "crossrow="<<crossrow<<
943 "errorpar="<<errorpar<<
944 "pointError="<<pointError<<
956 /////////////////////////
958 ////////////////////////
960 if(nbclusters < fNumberClusters) return kFALSE;
961 if(nbclusters > fNumberClustersf) return kFALSE;
962 if(pointError >= 0.3) return kFALSE;
963 if(crossrow == 1) return kTRUE;
965 ///////////////////////
967 //////////////////////
970 //Add to the linear fitter of the detector
971 if( TMath::Abs(snp) < 1.){
972 Double_t x = tnp-dzdx*tnt;
973 if(fLinearFitterDebugOn) {
974 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
975 fEntriesLinearFitter[fDetectorPreviousTrack]++;
977 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
981 fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
986 //____________Offine tracking in the AliTRDtracker_____________________________
987 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
990 // PRF width calibration
991 // Assume a Gaussian shape: determinate the position of the three pad clusters
992 // Fit with a straight line
993 // Take the fitted values for all the clusters (3 or 2 pad clusters)
994 // Fill the PRF as function of angle of the track
999 ///////////////////////////////////////////
1000 // Take the parameters of the track
1001 //////////////////////////////////////////
1002 // take now the snp, tnp and tgl from the track
1003 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1004 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1005 if( TMath::Abs(snp) < 1.){
1006 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1008 Double_t tgl = tracklet->GetTgl(); // dz/dl
1009 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1011 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1012 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1013 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1014 // at the end with correction due to linear fit
1015 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1016 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1018 ///////////////////////////////
1019 // Calculate tnp group shift
1020 ///////////////////////////////
1021 Bool_t echec = kFALSE;
1022 Double_t shift = 0.0;
1023 //Calculate the shift in x coresponding to this tnp
1024 if(fNgroupprf != 0.0){
1025 shift = -3.0*(fNgroupprf-1)-1.5;
1026 Double_t limithigh = -0.2*(fNgroupprf-1);
1027 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1029 while(tnp > limithigh){
1035 // do nothing if out of tnp range
1036 //printf("echec %d\n",(Int_t)echec);
1037 if(echec) return kFALSE;
1039 ///////////////////////
1041 //////////////////////
1043 Int_t nb3pc = 0; // number of three pads clusters used for fit
1044 // to see the difference between the fit and the 3 pad clusters position
1045 Double_t padPositions[AliTRDseedV1::kNtb];
1046 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1047 fLinearFitterTracklet->ClearPoints();
1049 //printf("loop clusters \n");
1050 ////////////////////////////
1051 // loop over the clusters
1052 ////////////////////////////
1053 AliTRDcluster *cl = 0x0;
1054 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1055 // reject shared clusters on pad row
1056 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1057 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1060 cl = tracklet->GetClusters(ic);
1062 Double_t time = cl->GetPadTime();
1063 if((time<=7) || (time>=21)) continue;
1064 Short_t *signals = cl->GetSignals();
1065 Float_t xcenter = 0.0;
1066 Bool_t echec1 = kTRUE;
1068 /////////////////////////////////////////////////////////////
1069 // Center 3 balanced: position with the center of the pad
1070 /////////////////////////////////////////////////////////////
1071 if ((((Float_t) signals[3]) > 0.0) &&
1072 (((Float_t) signals[2]) > 0.0) &&
1073 (((Float_t) signals[4]) > 0.0)) {
1075 // Security if the denomiateur is 0
1076 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1077 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1078 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1079 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1080 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1086 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1087 if(echec1) continue;
1089 ////////////////////////////////////////////////////////
1090 //if no echec1: calculate with the position of the pad
1091 // Position of the cluster
1092 // fill the linear fitter
1093 ///////////////////////////////////////////////////////
1094 Double_t padPosition = xcenter + cl->GetPadCol();
1095 padPositions[ic] = padPosition;
1097 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1102 //printf("Fin loop clusters \n");
1103 //////////////////////////////
1104 // fit with a straight line
1105 /////////////////////////////
1107 fLinearFitterTracklet->ClearPoints();
1110 fLinearFitterTracklet->Eval();
1112 fLinearFitterTracklet->GetParameters(line);
1113 Float_t pointError = -1.0;
1114 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1115 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1117 fLinearFitterTracklet->ClearPoints();
1119 //printf("PRF second loop \n");
1120 ////////////////////////////////////////////////
1121 // Fill the PRF: Second loop over clusters
1122 //////////////////////////////////////////////
1123 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1124 // reject shared clusters on pad row
1125 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1126 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1128 cl = tracklet->GetClusters(ic);
1131 Short_t *signals = cl->GetSignals(); // signal
1132 Double_t time = cl->GetPadTime(); // time bin
1133 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1134 Float_t padPos = cl->GetPadCol(); // middle pad
1135 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1136 Float_t ycenter = 0.0; // relative center charge
1137 Float_t ymin = 0.0; // relative left charge
1138 Float_t ymax = 0.0; // relative right charge
1140 ////////////////////////////////////////////////////////////////
1141 // Calculate ycenter, ymin and ymax even for two pad clusters
1142 ////////////////////////////////////////////////////////////////
1143 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1144 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1145 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1146 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1147 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1148 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1151 /////////////////////////
1152 // Calibration group
1153 ////////////////////////
1154 Int_t rowcl = cl->GetPadRow(); // row of cluster
1155 Int_t colcl = cl->GetPadCol(); // col of cluster
1156 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1157 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1158 Float_t xcl = cl->GetY(); // y cluster
1159 Float_t qcl = cl->GetQ(); // charge cluster
1160 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1161 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1162 Double_t xdiff = dpad; // reconstructed position constant
1163 Double_t x = dpad; // reconstructed position moved
1164 Float_t ep = pointError; // error of fit
1165 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1166 Float_t signal3 = (Float_t)signals[3]; // signal
1167 Float_t signal2 = (Float_t)signals[2]; // signal
1168 Float_t signal4 = (Float_t)signals[4]; // signal
1169 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1173 /////////////////////
1175 ////////////////////
1177 if(fDebugLevel > 0){
1178 if ( !fDebugStreamer ) {
1180 TDirectory *backup = gDirectory;
1181 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1182 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1187 Float_t y = ycenter;
1188 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1189 "caligroup="<<caligroup<<
1190 "detector="<<fDetectorPreviousTrack<<
1193 "npoints="<<nbclusters<<
1202 "padPosition="<<padPositions[ic]<<
1203 "padPosTracklet="<<padPosTracklet<<
1208 "signal1="<<signal1<<
1209 "signal2="<<signal2<<
1210 "signal3="<<signal3<<
1211 "signal4="<<signal4<<
1212 "signal5="<<signal5<<
1218 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1219 "caligroup="<<caligroup<<
1220 "detector="<<fDetectorPreviousTrack<<
1223 "npoints="<<nbclusters<<
1232 "padPosition="<<padPositions[ic]<<
1233 "padPosTracklet="<<padPosTracklet<<
1238 "signal1="<<signal1<<
1239 "signal2="<<signal2<<
1240 "signal3="<<signal3<<
1241 "signal4="<<signal4<<
1242 "signal5="<<signal5<<
1248 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1249 "caligroup="<<caligroup<<
1250 "detector="<<fDetectorPreviousTrack<<
1253 "npoints="<<nbclusters<<
1262 "padPosition="<<padPositions[ic]<<
1263 "padPosTracklet="<<padPosTracklet<<
1268 "signal1="<<signal1<<
1269 "signal2="<<signal2<<
1270 "signal3="<<signal3<<
1271 "signal4="<<signal4<<
1272 "signal5="<<signal5<<
1278 /////////////////////
1280 /////////////////////
1281 if(nbclusters < fNumberClusters) continue;
1282 if(nbclusters > fNumberClustersf) continue;
1283 if(nb3pc <= 5) continue;
1284 if((time >= 21) || (time < 7)) continue;
1285 if(TMath::Abs(qcl) < 80) continue;
1286 if( TMath::Abs(snp) > 1.) continue;
1289 ////////////////////////
1291 ///////////////////////
1293 if(TMath::Abs(dpad) < 1.5) {
1294 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1295 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1296 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1298 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1299 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1300 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1302 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1303 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1304 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1309 if(TMath::Abs(dpad) < 1.5) {
1310 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1311 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1313 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1314 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1315 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1317 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1318 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1319 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1322 } // second loop over clusters
1327 ///////////////////////////////////////////////////////////////////////////////////////
1328 // Pad row col stuff: see if masked or not
1329 ///////////////////////////////////////////////////////////////////////////////////////
1330 //_____________________________________________________________________________
1331 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1334 // See if we are not near a masked pad
1337 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1341 //_____________________________________________________________________________
1342 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1345 // See if we are not near a masked pad
1348 if (!IsPadOn(detector, col, row)) {
1349 fGoodTracklet = kFALSE;
1353 if (!IsPadOn(detector, col-1, row)) {
1354 fGoodTracklet = kFALSE;
1359 if (!IsPadOn(detector, col+1, row)) {
1360 fGoodTracklet = kFALSE;
1365 //_____________________________________________________________________________
1366 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1369 // Look in the choosen database if the pad is On.
1370 // If no the track will be "not good"
1373 // Get the parameter object
1374 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1376 AliInfo("Could not get calibDB");
1380 if (!cal->IsChamberGood(detector) ||
1381 cal->IsChamberNoData(detector) ||
1382 cal->IsPadMasked(detector,col,row)) {
1390 ///////////////////////////////////////////////////////////////////////////////////////
1391 // Calibration groups: calculate the number of groups, localise...
1392 ////////////////////////////////////////////////////////////////////////////////////////
1393 //_____________________________________________________________________________
1394 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1397 // Calculate the calibration group number for i
1400 // Row of the cluster and position in the pad groups
1402 if (fCalibraMode->GetNnZ(i) != 0) {
1403 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1407 // Col of the cluster and position in the pad groups
1409 if (fCalibraMode->GetNnRphi(i) != 0) {
1410 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1413 return posc*fCalibraMode->GetNfragZ(i)+posr;
1416 //____________________________________________________________________________________
1417 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1420 // Calculate the total number of calibration groups
1426 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1428 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1433 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1435 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1440 fCalibraMode->ModePadCalibration(2,i);
1441 fCalibraMode->ModePadFragmentation(0,2,0,i);
1442 fCalibraMode->SetDetChamb2(i);
1443 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1444 fCalibraMode->ModePadCalibration(0,i);
1445 fCalibraMode->ModePadFragmentation(0,0,0,i);
1446 fCalibraMode->SetDetChamb0(i);
1447 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1448 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1452 //_____________________________________________________________________________
1453 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1456 // Set the mode of calibration group in the z direction for the parameter i
1461 fCalibraMode->SetNz(i, Nz);
1464 AliInfo("You have to choose between 0 and 4");
1468 //_____________________________________________________________________________
1469 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1472 // Set the mode of calibration group in the rphi direction for the parameter i
1477 fCalibraMode->SetNrphi(i ,Nrphi);
1480 AliInfo("You have to choose between 0 and 6");
1485 //_____________________________________________________________________________
1486 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1489 // Set the mode of calibration group all together
1491 if(fVector2d == kTRUE) {
1492 AliInfo("Can not work with the vector method");
1495 fCalibraMode->SetAllTogether(i);
1499 //_____________________________________________________________________________
1500 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1503 // Set the mode of calibration group per supermodule
1505 if(fVector2d == kTRUE) {
1506 AliInfo("Can not work with the vector method");
1509 fCalibraMode->SetPerSuperModule(i);
1513 //____________Set the pad calibration variables for the detector_______________
1514 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1517 // For the detector calcul the first Xbins and set the number of row
1518 // and col pads per calibration groups, the number of calibration
1519 // groups in the detector.
1522 // first Xbins of the detector
1524 fCalibraMode->CalculXBins(detector,0);
1527 fCalibraMode->CalculXBins(detector,1);
1530 fCalibraMode->CalculXBins(detector,2);
1533 // fragmentation of idect
1534 for (Int_t i = 0; i < 3; i++) {
1535 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1536 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1537 , (Int_t) GetStack(detector)
1538 , (Int_t) GetSector(detector),i);
1544 //_____________________________________________________________________________
1545 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1548 // Should be between 0 and 6
1551 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1552 AliInfo("The number of groups must be between 0 and 6!");
1555 fNgroupprf = numberGroupsPRF;
1559 ///////////////////////////////////////////////////////////////////////////////////////////
1560 // Per tracklet: store or reset the info, fill the histos with the info
1561 //////////////////////////////////////////////////////////////////////////////////////////
1562 //_____________________________________________________________________________
1563 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(const AliTRDcluster *cl,const Double_t dqdl,const Int_t *group,const Int_t row,const Int_t col,const AliTRDcluster *cls)
1566 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1567 // Correct from the gain correction before
1568 // cls is shared cluster if any
1571 //printf("StoreInfoCHPHtrack\n");
1573 // time bin of the cluster not corrected
1574 Int_t time = cl->GetPadTime();
1575 Float_t charge = TMath::Abs(cl->GetQ());
1577 charge += TMath::Abs(cls->GetQ());
1578 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1581 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1583 //Correct for the gain coefficient used in the database for reconstruction
1584 Float_t correctthegain = 1.0;
1585 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1586 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1587 Float_t correction = 1.0;
1588 Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1589 // we divide with gain in AliTRDclusterizer::Transform...
1590 if( correctthegain > 0 ) normalisation /= correctthegain;
1593 // take dd/dl corrected from the angle of the track
1594 correction = dqdl / (normalisation);
1597 // Fill the fAmpTotal with the charge
1599 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1600 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1601 fAmpTotal[(Int_t) group[0]] += correction;
1605 // Fill the fPHPlace and value
1607 fPHPlace[time] = group[1];
1608 fPHValue[time] = charge;
1612 //____________Offine tracking in the AliTRDtracker_____________________________
1613 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1616 // Reset values per tracklet
1619 //Reset good tracklet
1620 fGoodTracklet = kTRUE;
1622 // Reset the fPHValue
1624 //Reset the fPHValue and fPHPlace
1625 for (Int_t k = 0; k < fTimeMax; k++) {
1631 // Reset the fAmpTotal where we put value
1633 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1638 //____________Offine tracking in the AliTRDtracker_____________________________
1639 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1642 // For the offline tracking or mcm tracklets
1643 // This function will be called in the functions UpdateHistogram...
1644 // to fill the info of a track for the relativ gain calibration
1647 Int_t nb = 0; // Nombre de zones traversees
1648 Int_t fd = -1; // Premiere zone non nulle
1649 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1651 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1653 if(nbclusters < fNumberClusters) return;
1654 if(nbclusters > fNumberClustersf) return;
1657 // Normalize with the number of clusters
1658 Double_t normalizeCst = fRelativeScale;
1659 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1661 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1663 // See if the track goes through different zones
1664 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1665 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1666 if (fAmpTotal[k] > 0.0) {
1667 totalcharge += fAmpTotal[k];
1675 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1681 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1683 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1684 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1687 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1691 if ((fAmpTotal[fd] > 0.0) &&
1692 (fAmpTotal[fd+1] > 0.0)) {
1693 // One of the two very big
1694 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1696 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1697 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1700 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1703 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1705 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1707 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1708 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1711 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1714 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1717 if (fCalibraMode->GetNfragZ(0) > 1) {
1718 if (fAmpTotal[fd] > 0.0) {
1719 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1720 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1721 // One of the two very big
1722 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1724 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1725 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1728 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1731 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1733 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1735 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1736 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1739 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1741 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1752 //____________Offine tracking in the AliTRDtracker_____________________________
1753 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1756 // For the offline tracking or mcm tracklets
1757 // This function will be called in the functions UpdateHistogram...
1758 // to fill the info of a track for the drift velocity calibration
1761 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1762 Int_t fd1 = -1; // Premiere zone non nulle
1763 Int_t fd2 = -1; // Deuxieme zone non nulle
1764 Int_t k1 = -1; // Debut de la premiere zone
1765 Int_t k2 = -1; // Debut de la seconde zone
1766 Int_t nbclusters = 0; // number of clusters
1770 // See if the track goes through different zones
1771 for (Int_t k = 0; k < fTimeMax; k++) {
1772 if (fPHValue[k] > 0.0) {
1778 if (fPHPlace[k] != fd1) {
1784 if (fPHPlace[k] != fd2) {
1791 // See if noise before and after
1792 if(fMaxCluster > 0) {
1793 if(fPHValue[0] > fMaxCluster) return;
1794 if(fTimeMax > fNbMaxCluster) {
1795 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1796 if(fPHValue[k] > fMaxCluster) return;
1801 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1803 if(nbclusters < fNumberClusters) return;
1804 if(nbclusters > fNumberClustersf) return;
1810 for (Int_t i = 0; i < fTimeMax; i++) {
1812 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1814 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1816 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1819 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1821 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1827 if ((fd1 == fd2+1) ||
1829 // One of the two fast all the think
1830 if (k2 > (k1+fDifference)) {
1831 //we choose to fill the fd1 with all the values
1833 for (Int_t i = 0; i < fTimeMax; i++) {
1835 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1837 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1841 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1843 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1848 if ((k2+fDifference) < fTimeMax) {
1849 //we choose to fill the fd2 with all the values
1851 for (Int_t i = 0; i < fTimeMax; i++) {
1853 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1855 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1859 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1861 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1867 // Two zones voisines sinon rien!
1868 if (fCalibraMode->GetNfragZ(1) > 1) {
1870 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1871 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1872 // One of the two fast all the think
1873 if (k2 > (k1+fDifference)) {
1874 //we choose to fill the fd1 with all the values
1876 for (Int_t i = 0; i < fTimeMax; i++) {
1878 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1880 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1884 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1886 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1891 if ((k2+fDifference) < fTimeMax) {
1892 //we choose to fill the fd2 with all the values
1894 for (Int_t i = 0; i < fTimeMax; i++) {
1896 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1898 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1902 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1904 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1911 // Two zones voisines sinon rien!
1913 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1914 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1915 // One of the two fast all the think
1916 if (k2 > (k1 + fDifference)) {
1917 //we choose to fill the fd1 with all the values
1919 for (Int_t i = 0; i < fTimeMax; i++) {
1921 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1923 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1927 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1929 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1934 if ((k2+fDifference) < fTimeMax) {
1935 //we choose to fill the fd2 with all the values
1937 for (Int_t i = 0; i < fTimeMax; i++) {
1939 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1941 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1945 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1947 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1959 //////////////////////////////////////////////////////////////////////////////////////////
1960 // DAQ process functions
1961 /////////////////////////////////////////////////////////////////////////////////////////
1962 //_____________________________________________________________________
1963 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
1966 // Event Processing loop - AliTRDrawStream
1968 // 0 timebin problem
1971 // Same algorithm as TestBeam but different reader
1974 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
1976 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
1977 digitsManager->CreateArrays();
1979 Int_t withInput = 1;
1981 Double_t phvalue[16][144][36];
1982 for(Int_t k = 0; k < 36; k++){
1983 for(Int_t j = 0; j < 16; j++){
1984 for(Int_t c = 0; c < 144; c++){
1985 phvalue[j][c][k] = 0.0;
1990 fDetectorPreviousTrack = -1;
1994 Int_t nbtimebin = 0;
1995 Int_t baseline = 10;
2001 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2003 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2004 // printf("there is ADC data on this chamber!\n");
2006 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2007 if (digits->HasData()) { //array
2009 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2010 if (indexes->IsAllocated() == kFALSE) {
2011 AliError("Indexes do not exist!");
2016 indexes->ResetCounters();
2018 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2019 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2020 //while (rawStream->Next()) {
2022 Int_t idetector = det; // current detector
2023 //Int_t imcm = rawStream->GetMCM(); // current MCM
2024 //Int_t irob = rawStream->GetROB(); // current ROB
2027 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2029 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2032 for(Int_t k = 0; k < 36; k++){
2033 for(Int_t j = 0; j < 16; j++){
2034 for(Int_t c = 0; c < 144; c++){
2035 phvalue[j][c][k] = 0.0;
2041 fDetectorPreviousTrack = idetector;
2042 //fMCMPrevious = imcm;
2043 //fROBPrevious = irob;
2045 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2046 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2047 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2048 baseline = digitParam->GetADCbaseline(det); // baseline
2050 if(nbtimebin == 0) return 0;
2051 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2052 fTimeMax = nbtimebin;
2054 fNumberClustersf = fTimeMax;
2055 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2058 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2059 // phvalue[row][col][itime] = signal[itime]-baseline;
2060 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2061 /*if(phvalue[iRow][iCol][itime] >= 20) {
2062 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2066 (Short_t)(digits->GetData(iRow,iCol,itime)),
2073 // fill the last one
2074 if(fDetectorPreviousTrack != -1){
2077 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2078 // printf("\n ---> withinput %d\n\n",withInput);
2080 for(Int_t k = 0; k < 36; k++){
2081 for(Int_t j = 0; j < 16; j++){
2082 for(Int_t c = 0; c < 144; c++){
2083 phvalue[j][c][k] = 0.0;
2091 digitsManager->ClearArrays(det);
2093 delete digitsManager;
2098 //_____________________________________________________________________
2099 //////////////////////////////////////////////////////////////////////////////
2100 // Routine inside the DAQ process
2101 /////////////////////////////////////////////////////////////////////////////
2102 //_______________________________________________________________________
2103 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2106 // Look for the maximum by collapsing over the time
2107 // Sum over four pad col and two pad row
2113 Int_t idect = fDetectorPreviousTrack;
2114 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2116 for(Int_t tb = 0; tb < 36; tb++){
2120 //fGoodTracklet = kTRUE;
2121 //fDetectorPreviousTrack = 0;
2124 ///////////////////////////
2126 /////////////////////////
2130 Double_t integralMax = -1;
2132 for (Int_t ir = 1; ir <= 15; ir++)
2134 for (Int_t ic = 2; ic <= 142; ic++)
2136 Double_t integral = 0;
2137 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2139 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2141 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2142 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2145 for(Int_t tb = 0; tb< fTimeMax; tb++){
2146 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2151 if (integralMax < integral)
2155 integralMax = integral;
2157 } // check max integral
2161 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2162 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2167 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2171 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2172 //if(!fGoodTracklet) used = 1;;
2174 // /////////////////////////////////////////////////////
2175 // sum ober 2 row and 4 pad cols for each time bins
2176 // ////////////////////////////////////////////////////
2180 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2182 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2184 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2185 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2187 for(Int_t it = 0; it < fTimeMax; it++){
2188 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2195 Double_t sumcharge = 0.0;
2196 for(Int_t it = 0; it < fTimeMax; it++){
2197 sumcharge += sum[it];
2198 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2202 /////////////////////////////////////////////////////////
2204 ////////////////////////////////////////////////////////
2205 if(fDebugLevel > 0){
2206 if ( !fDebugStreamer ) {
2208 TDirectory *backup = gDirectory;
2209 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2210 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2213 Double_t amph0 = sum[0];
2214 Double_t amphlast = sum[fTimeMax-1];
2215 Double_t rms = TMath::RMS(fTimeMax,sum);
2216 Int_t goodtracklet = (Int_t) fGoodTracklet;
2217 for(Int_t it = 0; it < fTimeMax; it++){
2218 Double_t clustera = sum[it];
2220 (* fDebugStreamer) << "FillDAQa"<<
2221 "ampTotal="<<sumcharge<<
2224 "detector="<<idect<<
2226 "amphlast="<<amphlast<<
2227 "goodtracklet="<<goodtracklet<<
2228 "clustera="<<clustera<<
2236 ////////////////////////////////////////////////////////
2238 ///////////////////////////////////////////////////////
2239 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2240 if(sum[0] > 100.0) used = 1;
2241 if(nbcl < fNumberClusters) used = 1;
2242 if(nbcl > fNumberClustersf) used = 1;
2244 //if(fDetectorPreviousTrack == 15){
2245 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2247 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2249 for(Int_t it = 0; it < fTimeMax; it++){
2250 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2252 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2254 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2256 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2261 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2263 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2270 //____________Online trackling in AliTRDtrigger________________________________
2271 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2275 // Fill a simple average pulse height
2279 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2285 //____________Write_____________________________________________________
2286 //_____________________________________________________________________
2287 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2290 // Write infos to file
2294 if ( fDebugStreamer ) {
2295 delete fDebugStreamer;
2296 fDebugStreamer = 0x0;
2299 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2304 ,fNumberUsedPh[1]));
2306 TDirectory *backup = gDirectory;
2312 option = "recreate";
2314 TFile f(filename,option.Data());
2316 TStopwatch stopwatch;
2319 f.WriteTObject(fCalibraVector);
2324 f.WriteTObject(fCH2d);
2329 f.WriteTObject(fPH2d);
2334 f.WriteTObject(fPRF2d);
2337 if(fLinearFitterOn){
2338 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2339 f.WriteTObject(fLinearVdriftFit);
2344 if ( backup ) backup->cd();
2346 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2347 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2349 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2351 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2352 //___________________________________________probe the histos__________________________________________________
2353 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2356 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2357 // debug mode with 2 for TH2I and 3 for TProfile2D
2358 // It gives a pointer to a Double_t[7] with the info following...
2359 // [0] : number of calibration groups with entries
2360 // [1] : minimal number of entries found
2361 // [2] : calibration group number of the min
2362 // [3] : maximal number of entries found
2363 // [4] : calibration group number of the max
2364 // [5] : mean number of entries found
2365 // [6] : mean relative error
2368 Double_t *info = new Double_t[7];
2370 // Number of Xbins (detectors or groups of pads)
2371 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2372 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2375 Double_t nbwe = 0; //number of calibration groups with entries
2376 Double_t minentries = 0; //minimal number of entries found
2377 Double_t maxentries = 0; //maximal number of entries found
2378 Double_t placemin = 0; //calibration group number of the min
2379 Double_t placemax = -1; //calibration group number of the max
2380 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2381 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2383 Double_t counter = 0;
2386 TH1F *nbEntries = 0x0;//distribution of the number of entries
2387 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2388 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2390 // Beginning of the loop over the calibration groups
2391 for (Int_t idect = 0; idect < nbins; idect++) {
2393 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2394 projch->SetDirectory(0);
2396 // Number of entries for this calibration group
2397 Double_t nentries = 0.0;
2399 for (Int_t k = 0; k < nxbins; k++) {
2400 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2404 for (Int_t k = 0; k < nxbins; k++) {
2405 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2406 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2407 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2416 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2417 nbEntries->SetDirectory(0);
2418 nbEntries->Fill(nentries);
2419 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2420 nbEntriesPerGroup->SetDirectory(0);
2421 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2422 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));
2423 nbEntriesPerSp->SetDirectory(0);
2424 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2429 if(nentries > maxentries){
2430 maxentries = nentries;
2434 minentries = nentries;
2436 if(nentries < minentries){
2437 minentries = nentries;
2443 meanstats += nentries;
2445 }//calibration groups loop
2447 if(nbwe > 0) meanstats /= nbwe;
2448 if(counter > 0) meanrelativerror /= counter;
2450 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2451 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2452 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2453 AliInfo(Form("The mean number of entries is %f",meanstats));
2454 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2457 info[1] = minentries;
2459 info[3] = maxentries;
2461 info[5] = meanstats;
2462 info[6] = meanrelativerror;
2464 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2465 gStyle->SetPalette(1);
2466 gStyle->SetOptStat(1111);
2467 gStyle->SetPadBorderMode(0);
2468 gStyle->SetCanvasColor(10);
2469 gStyle->SetPadLeftMargin(0.13);
2470 gStyle->SetPadRightMargin(0.01);
2471 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2474 nbEntries->Draw("");
2476 nbEntriesPerSp->SetStats(0);
2477 nbEntriesPerSp->Draw("");
2478 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2480 nbEntriesPerGroup->SetStats(0);
2481 nbEntriesPerGroup->Draw("");
2487 //____________________________________________________________________________
2488 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2491 // Return a Int_t[4] with:
2492 // 0 Mean number of entries
2493 // 1 median of number of entries
2494 // 2 rms of number of entries
2495 // 3 number of group with entries
2498 Double_t *stat = new Double_t[4];
2501 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2503 Double_t *weight = new Double_t[nbofgroups];
2504 Double_t *nonul = new Double_t[nbofgroups];
2506 for(Int_t k = 0; k < nbofgroups; k++){
2507 if(fEntriesCH[k] > 0) {
2509 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2512 else weight[k] = 0.0;
2514 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2515 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2516 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2524 //____________________________________________________________________________
2525 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2528 // Return a Int_t[4] with:
2529 // 0 Mean number of entries
2530 // 1 median of number of entries
2531 // 2 rms of number of entries
2532 // 3 number of group with entries
2535 Double_t *stat = new Double_t[4];
2538 Int_t nbofgroups = 540;
2539 Double_t *weight = new Double_t[nbofgroups];
2540 Int_t *nonul = new Int_t[nbofgroups];
2542 for(Int_t k = 0; k < nbofgroups; k++){
2543 if(fEntriesLinearFitter[k] > 0) {
2545 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2548 else weight[k] = 0.0;
2550 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2551 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2552 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2560 //////////////////////////////////////////////////////////////////////////////////////
2562 //////////////////////////////////////////////////////////////////////////////////////
2563 //_____________________________________________________________________________
2564 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2567 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2568 // If fNgroupprf is zero then no binning in tnp
2572 name += fCalibraMode->GetNz(2);
2574 name += fCalibraMode->GetNrphi(2);
2578 if(fNgroupprf != 0){
2580 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2581 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2582 fPRF2d->SetYTitle("Det/pad groups");
2583 fPRF2d->SetXTitle("Position x/W [pad width units]");
2584 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2585 fPRF2d->SetStats(0);
2588 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2589 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2590 fPRF2d->SetYTitle("Det/pad groups");
2591 fPRF2d->SetXTitle("Position x/W [pad width units]");
2592 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2593 fPRF2d->SetStats(0);
2598 //_____________________________________________________________________________
2599 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2602 // Create the 2D histos
2605 TString name("Ver");
2606 name += fVersionVdriftUsed;
2608 name += fSubVersionVdriftUsed;
2610 name += fFirstRunVdrift;
2612 name += fCalibraMode->GetNz(1);
2614 name += fCalibraMode->GetNrphi(1);
2616 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2617 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2619 fPH2d->SetYTitle("Det/pad groups");
2620 fPH2d->SetXTitle("time [#mus]");
2621 fPH2d->SetZTitle("<PH> [a.u.]");
2625 //_____________________________________________________________________________
2626 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2629 // Create the 2D histos
2632 TString name("Ver");
2633 name += fVersionGainUsed;
2635 name += fSubVersionGainUsed;
2637 name += fFirstRunGain;
2639 name += fCalibraMode->GetNz(0);
2641 name += fCalibraMode->GetNrphi(0);
2643 fCH2d = new TH2I("CH2d",(const Char_t *) name
2644 ,fNumberBinCharge,0,300,nn,0,nn);
2645 fCH2d->SetYTitle("Det/pad groups");
2646 fCH2d->SetXTitle("charge deposit [a.u]");
2647 fCH2d->SetZTitle("counts");
2652 //////////////////////////////////////////////////////////////////////////////////
2653 // Set relative scale
2654 /////////////////////////////////////////////////////////////////////////////////
2655 //_____________________________________________________________________________
2656 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2659 // Set the factor that will divide the deposited charge
2660 // to fit in the histo range [0,300]
2663 if (RelativeScale > 0.0) {
2664 fRelativeScale = RelativeScale;
2667 AliInfo("RelativeScale must be strict positif!");
2671 //////////////////////////////////////////////////////////////////////////////////
2672 // Quick way to fill a histo
2673 //////////////////////////////////////////////////////////////////////////////////
2674 //_____________________________________________________________________
2675 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2678 // FillCH2d: Marian style
2681 //skip simply the value out of range
2682 if((y>=300.0) || (y<0.0)) return;
2684 //Calcul the y place
2685 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2686 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2689 fCH2d->GetArray()[place]++;
2693 //////////////////////////////////////////////////////////////////////////////////
2694 // Geometrical functions
2695 ///////////////////////////////////////////////////////////////////////////////////
2696 //_____________________________________________________________________________
2697 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2700 // Reconstruct the layer number from the detector number
2703 return ((Int_t) (d % 6));
2707 //_____________________________________________________________________________
2708 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2711 // Reconstruct the stack number from the detector number
2713 const Int_t kNlayer = 6;
2715 return ((Int_t) (d % 30) / kNlayer);
2719 //_____________________________________________________________________________
2720 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2723 // Reconstruct the sector number from the detector number
2727 return ((Int_t) (d / fg));
2730 ///////////////////////////////////////////////////////////////////////////////////
2731 // Getter functions for DAQ of the CH2d and the PH2d
2732 //////////////////////////////////////////////////////////////////////////////////
2733 //_____________________________________________________________________
2734 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2737 // return pointer to fPH2d TProfile2D
2738 // create a new TProfile2D if it doesn't exist allready
2744 fTimeMax = nbtimebin;
2745 fSf = samplefrequency;
2751 //_____________________________________________________________________
2752 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2755 // return pointer to fCH2d TH2I
2756 // create a new TH2I if it doesn't exist allready
2765 ////////////////////////////////////////////////////////////////////////////////////////////
2766 // Drift velocity calibration
2767 ///////////////////////////////////////////////////////////////////////////////////////////
2768 //_____________________________________________________________________
2769 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2772 // return pointer to TLinearFitter Calibration
2773 // if force is true create a new TLinearFitter if it doesn't exist allready
2776 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2777 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2780 // if we are forced and TLinearFitter doesn't yet exist create it
2782 // new TLinearFitter
2783 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2784 fLinearFitterArray.AddAt(linearfitter,detector);
2785 return linearfitter;
2788 //____________________________________________________________________________
2789 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2792 // Analyse array of linear fitter because can not be written
2793 // Store two arrays: one with the param the other one with the error param + number of entries
2796 for(Int_t k = 0; k < 540; k++){
2797 TLinearFitter *linearfitter = GetLinearFitter(k);
2798 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2799 TVectorD *par = new TVectorD(2);
2800 TVectorD pare = TVectorD(2);
2801 TVectorD *parE = new TVectorD(3);
2802 if((linearfitter->EvalRobust(0.8)==0)) {
2803 //linearfitter->Eval();
2804 linearfitter->GetParameters(*par);
2805 //linearfitter->GetErrors(pare);
2806 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2807 //(*parE)[0] = pare[0]*ppointError;
2808 //(*parE)[1] = pare[1]*ppointError;
2812 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2813 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2814 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);