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 ,fNumberBinCharge(50)
191 ,fGoodTracklet(kTRUE)
192 ,fLinearFitterTracklet(0x0)
194 ,fEntriesLinearFitter(0x0)
199 ,fLinearFitterArray(540)
200 ,fLinearVdriftFit(0x0)
206 // Default constructor
210 // Init some default values
213 fNumberUsedCh[0] = 0;
214 fNumberUsedCh[1] = 0;
215 fNumberUsedPh[0] = 0;
216 fNumberUsedPh[1] = 0;
218 fGeo = new AliTRDgeometry();
219 fCalibDB = AliTRDcalibDB::Instance();
222 //______________________________________________________________________________________
223 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
230 ,fPRF2dOn(c.fPRF2dOn)
231 ,fHisto2d(c.fHisto2d)
232 ,fVector2d(c.fVector2d)
233 ,fLinearFitterOn(c.fLinearFitterOn)
234 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
235 ,fExbAltFitOn(c.fExbAltFitOn)
236 ,fScaleWithTPCSignal(c.fScaleWithTPCSignal)
237 ,fRelativeScale(c.fRelativeScale)
238 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
239 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
240 ,fFillWithZero(c.fFillWithZero)
241 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
242 ,fMaxCluster(c.fMaxCluster)
243 ,fNbMaxCluster(c.fNbMaxCluster)
244 ,fCutWithVdriftCalib(c.fCutWithVdriftCalib)
245 ,fMinNbTRDtracklets(c.fMinNbTRDtracklets)
246 ,fMinTRDMomentum(c.fMinTRDMomentum)
247 ,fFirstRunGain(c.fFirstRunGain)
248 ,fVersionGainUsed(c.fVersionGainUsed)
249 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
250 ,fFirstRunGainLocal(c.fFirstRunGainLocal)
251 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
252 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
253 ,fFirstRunVdrift(c.fFirstRunVdrift)
254 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
255 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
256 ,fFirstRunExB(c.fFirstRunExB)
257 ,fVersionExBUsed(c.fVersionExBUsed)
258 ,fSubVersionExBUsed(c.fSubVersionExBUsed)
261 ,fDebugLevel(c.fDebugLevel)
262 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
263 ,fMCMPrevious(c.fMCMPrevious)
264 ,fROBPrevious(c.fROBPrevious)
265 ,fNumberClusters(c.fNumberClusters)
266 ,fNumberClustersf(c.fNumberClustersf)
267 ,fNumberClustersProcent(c.fNumberClustersProcent)
268 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
269 ,fNumberRowDAQ(c.fNumberRowDAQ)
270 ,fNumberColDAQ(c.fNumberColDAQ)
271 ,fProcent(c.fProcent)
272 ,fDifference(c.fDifference)
273 ,fNumberTrack(c.fNumberTrack)
274 ,fTimeMax(c.fTimeMax)
276 ,fNumberBinCharge(c.fNumberBinCharge)
277 ,fNumberBinPRF(c.fNumberBinPRF)
278 ,fNgroupprf(c.fNgroupprf)
282 ,fGoodTracklet(c.fGoodTracklet)
283 ,fLinearFitterTracklet(0x0)
285 ,fEntriesLinearFitter(0x0)
290 ,fLinearFitterArray(540)
291 ,fLinearVdriftFit(0x0)
299 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
300 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
302 fPH2d = new TProfile2D(*c.fPH2d);
303 fPH2d->SetDirectory(0);
306 fPRF2d = new TProfile2D(*c.fPRF2d);
307 fPRF2d->SetDirectory(0);
310 fCH2d = new TH2I(*c.fCH2d);
311 fCH2d->SetDirectory(0);
313 if(c.fLinearVdriftFit){
314 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
317 fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
320 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
321 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
326 fGeo = new AliTRDgeometry();
327 fCalibDB = AliTRDcalibDB::Instance();
329 fNumberUsedCh[0] = 0;
330 fNumberUsedCh[1] = 0;
331 fNumberUsedPh[0] = 0;
332 fNumberUsedPh[1] = 0;
336 //____________________________________________________________________________________
337 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
340 // AliTRDCalibraFillHisto destructor
344 if ( fDebugStreamer ) delete fDebugStreamer;
346 if ( fCalDetGain ) delete fCalDetGain;
347 if ( fCalROCGain ) delete fCalROCGain;
349 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
353 delete [] fEntriesCH;
354 delete [] fEntriesLinearFitter;
357 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
358 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
361 if(fLinearVdriftFit) delete fLinearVdriftFit;
362 if(fExbAltFit) delete fExbAltFit;
368 //_____________________________________________________________________________
369 void AliTRDCalibraFillHisto::Destroy()
380 //_____________________________________________________________________________
381 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
384 // Delete DebugStreamer
387 if ( fDebugStreamer ) delete fDebugStreamer;
388 fDebugStreamer = 0x0;
391 //_____________________________________________________________________________
392 void AliTRDCalibraFillHisto::ClearHistos()
412 //////////////////////////////////////////////////////////////////////////////////
413 // calibration with AliTRDtrackV1: Init, Update
414 //////////////////////////////////////////////////////////////////////////////////
415 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
416 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
419 // Init the histograms and stuff to be filled
424 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
426 AliInfo("Could not get calibDB");
429 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
431 AliInfo("Could not get CommonParam");
436 if(nboftimebin > 0) fTimeMax = nboftimebin;
437 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
438 if(fTimeMax <= 0) fTimeMax = 30;
439 printf("////////////////////////////////////////////\n");
440 printf("Number of time bins in calibration component %d\n",fTimeMax);
441 printf("////////////////////////////////////////////\n");
442 fSf = parCom->GetSamplingFrequency();
443 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
444 else fRelativeScale = 1.18;
445 fNumberClustersf = fTimeMax;
446 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
448 // Init linear fitter
449 if(!fLinearFitterTracklet) {
450 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
451 fLinearFitterTracklet->StoreData(kTRUE);
454 // Calcul Xbins Chambd0, Chamb2
455 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
456 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
457 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
459 // If vector method On initialised all the stuff
461 fCalibraVector = new AliTRDCalibraVector();
462 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
463 fCalibraVector->SetTimeMax(fTimeMax);
464 if(fNgroupprf != 0) {
465 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
466 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
469 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
470 fCalibraVector->SetPRFRange(1.5);
472 for(Int_t k = 0; k < 3; k++){
473 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
474 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
476 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
477 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
478 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
479 fCalibraVector->SetNbGroupPRF(fNgroupprf);
482 // Create the 2D histos corresponding to the pad groupCalibration mode
485 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
486 ,fCalibraMode->GetNz(0)
487 ,fCalibraMode->GetNrphi(0)));
489 // Create the 2D histo
494 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
495 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
499 fEntriesCH = new Int_t[ntotal0];
500 for(Int_t k = 0; k < ntotal0; k++){
507 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
508 ,fCalibraMode->GetNz(1)
509 ,fCalibraMode->GetNrphi(1)));
511 // Create the 2D histo
516 fPHPlace = new Short_t[fTimeMax];
517 for (Int_t k = 0; k < fTimeMax; k++) {
520 fPHValue = new Float_t[fTimeMax];
521 for (Int_t k = 0; k < fTimeMax; k++) {
525 if (fLinearFitterOn) {
526 if(fLinearFitterDebugOn) {
527 fLinearFitterArray.SetName("ArrayLinearFitters");
528 fEntriesLinearFitter = new Int_t[540];
529 for(Int_t k = 0; k < 540; k++){
530 fEntriesLinearFitter[k] = 0;
533 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
534 TString nameee("Ver");
535 nameee += fVersionExBUsed;
537 nameee += fSubVersionExBUsed;
538 nameee += "FirstRun";
539 nameee += fFirstRunExB;
541 fLinearVdriftFit->SetNameCalibUsed(nameee);
544 fExbAltFit = new AliTRDCalibraExbAltFit();
549 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
550 ,fCalibraMode->GetNz(2)
551 ,fCalibraMode->GetNrphi(2)));
552 // Create the 2D histo
554 CreatePRF2d(ntotal2);
561 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
562 Bool_t AliTRDCalibraFillHisto::InitCalDet()
565 // Init the Gain Cal Det
570 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
572 AliError("No gain det calibration entry found");
575 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
577 AliError("No calDet gain found");
583 fCalDetGain->~AliTRDCalDet();
584 new(fCalDetGain) AliTRDCalDet(*(calDet));
585 }else fCalDetGain = new AliTRDCalDet(*(calDet));
590 name += fVersionGainUsed;
592 name += fSubVersionGainUsed;
594 name += fFirstRunGain;
596 name += fCalibraMode->GetNz(0);
598 name += fCalibraMode->GetNrphi(0);
600 fCH2d->SetTitle(name);
603 TString namee("Ver");
604 namee += fVersionVdriftUsed;
606 namee += fSubVersionVdriftUsed;
608 namee += fFirstRunVdrift;
610 namee += fCalibraMode->GetNz(1);
612 namee += fCalibraMode->GetNrphi(1);
614 fPH2d->SetTitle(namee);
616 // title AliTRDCalibraVdriftLinearFit
617 TString nameee("Ver");
618 nameee += fVersionExBUsed;
620 nameee += fSubVersionExBUsed;
621 nameee += "FirstRun";
622 nameee += fFirstRunExB;
626 fLinearVdriftFit->SetNameCalibUsed(nameee);
633 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
634 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
637 // Init the Gain Cal Pad
642 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainLocalUsed,fSubVersionGainLocalUsed);
644 AliError("No gain pad calibration entry found");
647 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
649 AliError("No calPad gain found");
652 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
654 AliError("No calRoc gain found");
659 fCalROCGain->~AliTRDCalROC();
660 new(fCalROCGain) AliTRDCalROC(*(calRoc));
661 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
670 //____________Offline tracking in the AliTRDtracker____________________________
671 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t,const AliESDtrack *esdtrack)
674 // Use AliTRDtrackV1 for the calibration
678 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
679 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
680 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
681 Bool_t newtr = kTRUE; // new track
685 // Cut on the number of TRD tracklets
687 Int_t numberoftrdtracklets = t->GetNumberOfTracklets();
688 if(numberoftrdtracklets < fMinNbTRDtracklets) return kFALSE;
690 Double_t tpcsignal = 1.0;
691 if(esdtrack) tpcsignal = esdtrack->GetTPCsignal()/50.0;
692 if(fScaleWithTPCSignal && tpcsignal <0.00001) return kFALSE;
696 AliInfo("Could not get calibDB");
701 ///////////////////////////
702 // loop over the tracklet
703 ///////////////////////////
704 for(Int_t itr = 0; itr < 6; itr++){
706 if(!(tracklet = t->GetTracklet(itr))) continue;
707 if(!tracklet->IsOK()) continue;
709 ResetfVariablestracklet();
710 Float_t momentum = t->GetMomentum(itr);
711 if(TMath::Abs(momentum) < fMinTRDMomentum) continue;
714 //////////////////////////////////////////
715 // localisation of the tracklet and dqdl
716 //////////////////////////////////////////
717 Int_t layer = tracklet->GetPlane();
719 while(!(cl = tracklet->GetClusters(ic++))) continue;
720 Int_t detector = cl->GetDetector();
721 if (detector != fDetectorPreviousTrack) {
722 // if not a new track
724 // don't use the rest of this track if in the same plane
725 if (layer == GetLayer(fDetectorPreviousTrack)) {
726 //printf("bad tracklet, same layer for detector %d\n",detector);
730 //Localise the detector bin
731 LocalisationDetectorXbins(detector);
733 if(!fIsHLT) InitCalPad(detector);
736 fDetectorPreviousTrack = detector;
740 ////////////////////////////
741 // loop over the clusters
742 ////////////////////////////
743 Double_t chargeQ = 0.0;
744 Int_t nbclusters = 0;
745 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
746 if(!(cl = tracklet->GetClusters(jc))) continue;
749 // Store the info bis of the tracklet
750 Int_t row = cl->GetPadRow();
751 Int_t col = cl->GetPadCol();
752 CheckGoodTrackletV1(cl);
753 Int_t group[2] = {0,0};
754 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
755 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
756 // Add the charge if shared cluster
757 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
759 //Scale with TPC signal or not
760 if(!fScaleWithTPCSignal) {
761 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
762 //printf("Do not scale now\n");
765 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc)/tpcsignal,group,row,col,cls); //tracklet->GetdQdl(jc)
770 ////////////////////////////////////////
771 // Fill the stuffs if a good tracklet
772 ////////////////////////////////////////
775 // drift velocity unables to cut bad tracklets
776 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
778 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
782 if(fCutWithVdriftCalib) {
783 if(pass) FillTheInfoOfTheTrackCH(nbclusters);
785 FillTheInfoOfTheTrackCH(nbclusters);
791 if(fCutWithVdriftCalib) {
792 if(pass) FillTheInfoOfTheTrackPH();
795 FillTheInfoOfTheTrackPH();
799 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
802 /////////////////////////////////////////////////////////
804 ////////////////////////////////////////////////////////
807 if ( !fDebugStreamer ) {
809 TDirectory *backup = gDirectory;
810 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
811 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
814 Int_t stacke = AliTRDgeometry::GetStack(detector);
815 Int_t sme = AliTRDgeometry::GetSector(detector);
816 Int_t layere = AliTRDgeometry::GetLayer(detector);
818 Float_t b[2] = {0.0,0.0};
819 Float_t bCov[3] = {0.0,0.0,0.0};
820 if(esdtrack) esdtrack->GetImpactParameters(b,bCov);
821 if (bCov[0]<=0 || bCov[2]<=0) {
822 bCov[0]=0; bCov[2]=0;
824 Float_t dcaxy = b[0];
826 Int_t tpcnbclusters = 0;
827 if(esdtrack) tpcnbclusters = esdtrack->GetTPCclusters(0);
828 Double_t ttpcsignal = 0.0;
829 if(esdtrack) ttpcsignal = esdtrack->GetTPCsignal();
830 Int_t cutvdriftlinear = 0;
831 if(!pass) cutvdriftlinear = 1;
833 (* fDebugStreamer) << "FillCharge"<<
834 "detector="<<detector<<
840 "nbtpccls="<<tpcnbclusters<<
841 "tpcsignal="<<ttpcsignal<<
842 "cutvdriftlinear="<<cutvdriftlinear<<
844 "nbtrdtracklet="<<numberoftrdtracklets<<
850 } // if a good tracklet
856 ///////////////////////////////////////////////////////////////////////////////////
857 // Routine inside the update with AliTRDtrack
858 ///////////////////////////////////////////////////////////////////////////////////
859 //____________Offine tracking in the AliTRDtracker_____________________________
860 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
863 // Drift velocity calibration:
864 // Fit the clusters with a straight line
865 // From the slope find the drift velocity
868 ////////////////////////////////////////////////
869 //Number of points: if less than 3 return kFALSE
870 /////////////////////////////////////////////////
871 if(nbclusters <= 2) return kFALSE;
876 // results of the linear fit
877 Double_t dydt = 0.0; // dydt tracklet after straight line fit
878 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
879 Double_t pointError = 0.0; // error after straight line fit
880 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
881 Int_t crossrow = 0; // if it crosses a pad row
882 Int_t rowp = -1; // if it crosses a pad row
883 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
884 fLinearFitterTracklet->ClearPoints();
887 ///////////////////////////////////////////
888 // Take the parameters of the track
889 //////////////////////////////////////////
890 // take now the snp, tnp and tgl from the track
891 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
892 Double_t tnp = 0.0; // dy/dx at the end of the chamber
893 if( TMath::Abs(snp) < 1.){
894 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
896 Double_t tgl = tracklet->GetTgl(); // dz/dl
897 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
899 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
900 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
901 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
902 // at the end with correction due to linear fit
903 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
904 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
907 ////////////////////////////
908 // loop over the clusters
909 ////////////////////////////
911 AliTRDcluster *cl = 0x0;
912 //////////////////////////////
913 // Check no shared clusters
914 //////////////////////////////
915 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
916 cl = tracklet->GetClusters(icc);
919 //////////////////////////////////
921 //////////////////////////////////
923 Float_t sigArr[AliTRDfeeParam::GetNcol()];
924 memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
925 Int_t ncl=0, tbf=0, tbl=0;
927 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
928 if(!(cl = tracklet->GetClusters(ic))) continue;
933 Int_t col = cl->GetPadCol();
934 for(int ip=-1, jp=2; jp<5; ip++, jp++){
936 if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
937 sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
940 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
942 Double_t ycluster = cl->GetY();
943 Int_t time = cl->GetPadTime();
944 Double_t timeis = time/fSf;
945 //See if cross two pad rows
946 Int_t row = cl->GetPadRow();
947 if(rowp==-1) rowp = row;
948 if(row != rowp) crossrow = 1;
950 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
956 ////////////////////////////////////
957 // Do the straight line fit now
958 ///////////////////////////////////
960 fLinearFitterTracklet->ClearPoints();
964 fLinearFitterTracklet->Eval();
965 fLinearFitterTracklet->GetParameters(pars);
966 pointError = TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
967 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
969 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
970 fLinearFitterTracklet->ClearPoints();
972 ////////////////////////////////////
973 // Calc the projection of the clusters on the y direction
974 ///////////////////////////////////
976 Float_t signalSum(0.);
977 Float_t mean = 0.0, rms = 0.0;
978 Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
979 Float_t dz = dzdx*(tbl-tbf)/10;
981 for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
982 signalSum+=sigArr[ip];
985 if(signalSum > 0.0) mean/=signalSum;
987 for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
988 rms+=sigArr[ip]*(ip-mean)*(ip-mean);
990 if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
992 rms -= TMath::Abs(dz*tilt);
996 ////////////////////////////////
998 ///////////////////////////////
1001 if(fDebugLevel > 1){
1002 if ( !fDebugStreamer ) {
1004 TDirectory *backup = gDirectory;
1005 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1006 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1009 float xcoord = tnp-dzdx*tnt;
1010 float pt = tracklet->GetPt();
1011 Int_t layer = GetLayer(fDetectorPreviousTrack);
1013 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1014 //"snpright="<<snpright<<
1016 "nbclusters="<<nbclusters<<
1017 "detector="<<fDetectorPreviousTrack<<
1025 "crossrow="<<crossrow<<
1026 "errorpar="<<errorpar<<
1027 "pointError="<<pointError<<
1039 /////////////////////////
1041 ////////////////////////
1043 if(nbclusters < fNumberClusters) return kFALSE;
1044 if(nbclusters > fNumberClustersf) return kFALSE;
1045 if(pointError >= 0.3) return kFALSE;
1046 if(crossrow == 1) return kTRUE;
1048 ///////////////////////
1050 //////////////////////
1052 if(fLinearFitterOn){
1053 //Add to the linear fitter of the detector
1054 if( TMath::Abs(snp) < 1.){
1055 Double_t x = tnp-dzdx*tnt;
1056 if(fLinearFitterDebugOn) {
1057 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1058 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1060 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1064 fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
1069 //____________Offine tracking in the AliTRDtracker_____________________________
1070 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1073 // PRF width calibration
1074 // Assume a Gaussian shape: determinate the position of the three pad clusters
1075 // Fit with a straight line
1076 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1077 // Fill the PRF as function of angle of the track
1081 //printf("begin\n");
1082 ///////////////////////////////////////////
1083 // Take the parameters of the track
1084 //////////////////////////////////////////
1085 // take now the snp, tnp and tgl from the track
1086 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1087 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1088 if( TMath::Abs(snp) < 1.){
1089 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1091 Double_t tgl = tracklet->GetTgl(); // dz/dl
1092 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1094 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1095 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1096 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1097 // at the end with correction due to linear fit
1098 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1099 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1101 ///////////////////////////////
1102 // Calculate tnp group shift
1103 ///////////////////////////////
1104 Bool_t echec = kFALSE;
1105 Double_t shift = 0.0;
1106 //Calculate the shift in x coresponding to this tnp
1107 if(fNgroupprf != 0.0){
1108 shift = -3.0*(fNgroupprf-1)-1.5;
1109 Double_t limithigh = -0.2*(fNgroupprf-1);
1110 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1112 while(tnp > limithigh){
1118 // do nothing if out of tnp range
1119 //printf("echec %d\n",(Int_t)echec);
1120 if(echec) return kFALSE;
1122 ///////////////////////
1124 //////////////////////
1126 Int_t nb3pc = 0; // number of three pads clusters used for fit
1127 // to see the difference between the fit and the 3 pad clusters position
1128 Double_t padPositions[AliTRDseedV1::kNtb];
1129 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1130 fLinearFitterTracklet->ClearPoints();
1132 //printf("loop clusters \n");
1133 ////////////////////////////
1134 // loop over the clusters
1135 ////////////////////////////
1136 AliTRDcluster *cl = 0x0;
1137 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1138 // reject shared clusters on pad row
1139 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1140 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1143 cl = tracklet->GetClusters(ic);
1145 Double_t time = cl->GetPadTime();
1146 if((time<=7) || (time>=21)) continue;
1147 Short_t *signals = cl->GetSignals();
1148 Float_t xcenter = 0.0;
1149 Bool_t echec1 = kTRUE;
1151 /////////////////////////////////////////////////////////////
1152 // Center 3 balanced: position with the center of the pad
1153 /////////////////////////////////////////////////////////////
1154 if ((((Float_t) signals[3]) > 0.0) &&
1155 (((Float_t) signals[2]) > 0.0) &&
1156 (((Float_t) signals[4]) > 0.0)) {
1158 // Security if the denomiateur is 0
1159 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1160 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1161 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1162 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1163 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1169 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1170 if(echec1) continue;
1172 ////////////////////////////////////////////////////////
1173 //if no echec1: calculate with the position of the pad
1174 // Position of the cluster
1175 // fill the linear fitter
1176 ///////////////////////////////////////////////////////
1177 Double_t padPosition = xcenter + cl->GetPadCol();
1178 padPositions[ic] = padPosition;
1180 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1185 //printf("Fin loop clusters \n");
1186 //////////////////////////////
1187 // fit with a straight line
1188 /////////////////////////////
1190 fLinearFitterTracklet->ClearPoints();
1193 fLinearFitterTracklet->Eval();
1195 fLinearFitterTracklet->GetParameters(line);
1196 Float_t pointError = -1.0;
1197 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1198 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1200 fLinearFitterTracklet->ClearPoints();
1202 //printf("PRF second loop \n");
1203 ////////////////////////////////////////////////
1204 // Fill the PRF: Second loop over clusters
1205 //////////////////////////////////////////////
1206 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1207 // reject shared clusters on pad row
1208 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1209 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1211 cl = tracklet->GetClusters(ic);
1214 Short_t *signals = cl->GetSignals(); // signal
1215 Double_t time = cl->GetPadTime(); // time bin
1216 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1217 Float_t padPos = cl->GetPadCol(); // middle pad
1218 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1219 Float_t ycenter = 0.0; // relative center charge
1220 Float_t ymin = 0.0; // relative left charge
1221 Float_t ymax = 0.0; // relative right charge
1223 ////////////////////////////////////////////////////////////////
1224 // Calculate ycenter, ymin and ymax even for two pad clusters
1225 ////////////////////////////////////////////////////////////////
1226 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1227 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1228 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1229 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1230 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1231 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1234 /////////////////////////
1235 // Calibration group
1236 ////////////////////////
1237 Int_t rowcl = cl->GetPadRow(); // row of cluster
1238 Int_t colcl = cl->GetPadCol(); // col of cluster
1239 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1240 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1241 Float_t xcl = cl->GetY(); // y cluster
1242 Float_t qcl = cl->GetQ(); // charge cluster
1243 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1244 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1245 Double_t xdiff = dpad; // reconstructed position constant
1246 Double_t x = dpad; // reconstructed position moved
1247 Float_t ep = pointError; // error of fit
1248 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1249 Float_t signal3 = (Float_t)signals[3]; // signal
1250 Float_t signal2 = (Float_t)signals[2]; // signal
1251 Float_t signal4 = (Float_t)signals[4]; // signal
1252 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1256 /////////////////////
1258 ////////////////////
1260 if(fDebugLevel > 1){
1261 if ( !fDebugStreamer ) {
1263 TDirectory *backup = gDirectory;
1264 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1265 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1270 Float_t y = ycenter;
1271 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1272 "caligroup="<<caligroup<<
1273 "detector="<<fDetectorPreviousTrack<<
1276 "npoints="<<nbclusters<<
1285 "padPosition="<<padPositions[ic]<<
1286 "padPosTracklet="<<padPosTracklet<<
1291 "signal1="<<signal1<<
1292 "signal2="<<signal2<<
1293 "signal3="<<signal3<<
1294 "signal4="<<signal4<<
1295 "signal5="<<signal5<<
1301 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1302 "caligroup="<<caligroup<<
1303 "detector="<<fDetectorPreviousTrack<<
1306 "npoints="<<nbclusters<<
1315 "padPosition="<<padPositions[ic]<<
1316 "padPosTracklet="<<padPosTracklet<<
1321 "signal1="<<signal1<<
1322 "signal2="<<signal2<<
1323 "signal3="<<signal3<<
1324 "signal4="<<signal4<<
1325 "signal5="<<signal5<<
1331 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1332 "caligroup="<<caligroup<<
1333 "detector="<<fDetectorPreviousTrack<<
1336 "npoints="<<nbclusters<<
1345 "padPosition="<<padPositions[ic]<<
1346 "padPosTracklet="<<padPosTracklet<<
1351 "signal1="<<signal1<<
1352 "signal2="<<signal2<<
1353 "signal3="<<signal3<<
1354 "signal4="<<signal4<<
1355 "signal5="<<signal5<<
1361 /////////////////////
1363 /////////////////////
1364 if(nbclusters < fNumberClusters) continue;
1365 if(nbclusters > fNumberClustersf) continue;
1366 if(nb3pc <= 5) continue;
1367 if((time >= 21) || (time < 7)) continue;
1368 if(TMath::Abs(qcl) < 80) continue;
1369 if( TMath::Abs(snp) > 1.) continue;
1372 ////////////////////////
1374 ///////////////////////
1376 if(TMath::Abs(dpad) < 1.5) {
1377 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1378 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1379 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1381 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1382 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1383 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1385 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1386 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1387 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1392 if(TMath::Abs(dpad) < 1.5) {
1393 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1394 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1396 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1397 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1398 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1400 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1401 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1402 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1405 } // second loop over clusters
1410 ///////////////////////////////////////////////////////////////////////////////////////
1411 // Pad row col stuff: see if masked or not
1412 ///////////////////////////////////////////////////////////////////////////////////////
1413 //_____________________________________________________________________________
1414 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1417 // See if we are not near a masked pad
1420 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1424 //_____________________________________________________________________________
1425 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1428 // See if we are not near a masked pad
1431 if (!IsPadOn(detector, col, row)) {
1432 fGoodTracklet = kFALSE;
1436 if (!IsPadOn(detector, col-1, row)) {
1437 fGoodTracklet = kFALSE;
1442 if (!IsPadOn(detector, col+1, row)) {
1443 fGoodTracklet = kFALSE;
1448 //_____________________________________________________________________________
1449 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1452 // Look in the choosen database if the pad is On.
1453 // If no the track will be "not good"
1456 // Get the parameter object
1457 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1459 AliInfo("Could not get calibDB");
1463 if (!cal->IsChamberGood(detector) ||
1464 cal->IsChamberNoData(detector) ||
1465 cal->IsPadMasked(detector,col,row)) {
1473 ///////////////////////////////////////////////////////////////////////////////////////
1474 // Calibration groups: calculate the number of groups, localise...
1475 ////////////////////////////////////////////////////////////////////////////////////////
1476 //_____________________________________________________________________________
1477 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1480 // Calculate the calibration group number for i
1483 // Row of the cluster and position in the pad groups
1485 if (fCalibraMode->GetNnZ(i) != 0) {
1486 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1490 // Col of the cluster and position in the pad groups
1492 if (fCalibraMode->GetNnRphi(i) != 0) {
1493 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1496 return posc*fCalibraMode->GetNfragZ(i)+posr;
1499 //____________________________________________________________________________________
1500 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1503 // Calculate the total number of calibration groups
1509 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1511 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1516 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1518 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1523 fCalibraMode->ModePadCalibration(2,i);
1524 fCalibraMode->ModePadFragmentation(0,2,0,i);
1525 fCalibraMode->SetDetChamb2(i);
1526 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1527 fCalibraMode->ModePadCalibration(0,i);
1528 fCalibraMode->ModePadFragmentation(0,0,0,i);
1529 fCalibraMode->SetDetChamb0(i);
1530 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1531 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1535 //_____________________________________________________________________________
1536 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1539 // Set the mode of calibration group in the z direction for the parameter i
1544 fCalibraMode->SetNz(i, Nz);
1547 AliInfo("You have to choose between 0 and 4");
1551 //_____________________________________________________________________________
1552 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1555 // Set the mode of calibration group in the rphi direction for the parameter i
1560 fCalibraMode->SetNrphi(i ,Nrphi);
1563 AliInfo("You have to choose between 0 and 6");
1568 //_____________________________________________________________________________
1569 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1572 // Set the mode of calibration group all together
1574 if(fVector2d == kTRUE) {
1575 AliInfo("Can not work with the vector method");
1578 fCalibraMode->SetAllTogether(i);
1582 //_____________________________________________________________________________
1583 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1586 // Set the mode of calibration group per supermodule
1588 if(fVector2d == kTRUE) {
1589 AliInfo("Can not work with the vector method");
1592 fCalibraMode->SetPerSuperModule(i);
1596 //____________Set the pad calibration variables for the detector_______________
1597 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1600 // For the detector calcul the first Xbins and set the number of row
1601 // and col pads per calibration groups, the number of calibration
1602 // groups in the detector.
1605 // first Xbins of the detector
1607 fCalibraMode->CalculXBins(detector,0);
1610 fCalibraMode->CalculXBins(detector,1);
1613 fCalibraMode->CalculXBins(detector,2);
1616 // fragmentation of idect
1617 for (Int_t i = 0; i < 3; i++) {
1618 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1619 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1620 , (Int_t) GetStack(detector)
1621 , (Int_t) GetSector(detector),i);
1627 //_____________________________________________________________________________
1628 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1631 // Should be between 0 and 6
1634 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1635 AliInfo("The number of groups must be between 0 and 6!");
1638 fNgroupprf = numberGroupsPRF;
1642 ///////////////////////////////////////////////////////////////////////////////////////////
1643 // Per tracklet: store or reset the info, fill the histos with the info
1644 //////////////////////////////////////////////////////////////////////////////////////////
1645 //_____________________________________________________________________________
1646 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)
1649 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1650 // Correct from the gain correction before
1651 // cls is shared cluster if any
1652 // Return the charge
1656 //printf("StoreInfoCHPHtrack\n");
1658 // time bin of the cluster not corrected
1659 Int_t time = cl->GetPadTime();
1660 Float_t charge = TMath::Abs(cl->GetQ());
1662 charge += TMath::Abs(cls->GetQ());
1663 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1666 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1668 //Correct for the gain coefficient used in the database for reconstruction
1669 Float_t correctthegain = 1.0;
1670 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1671 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1672 Float_t correction = 1.0;
1673 Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1674 // we divide with gain in AliTRDclusterizer::Transform...
1675 if( correctthegain > 0 ) normalisation /= correctthegain;
1678 // take dd/dl corrected from the angle of the track
1679 correction = dqdl / (normalisation);
1682 // Fill the fAmpTotal with the charge
1684 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1685 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1686 fAmpTotal[(Int_t) group[0]] += correction;
1690 // Fill the fPHPlace and value
1692 fPHPlace[time] = group[1];
1693 fPHValue[time] = charge;
1699 //____________Offine tracking in the AliTRDtracker_____________________________
1700 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1703 // Reset values per tracklet
1706 //Reset good tracklet
1707 fGoodTracklet = kTRUE;
1709 // Reset the fPHValue
1711 //Reset the fPHValue and fPHPlace
1712 for (Int_t k = 0; k < fTimeMax; k++) {
1718 // Reset the fAmpTotal where we put value
1720 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1725 //____________Offine tracking in the AliTRDtracker_____________________________
1726 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1729 // For the offline tracking or mcm tracklets
1730 // This function will be called in the functions UpdateHistogram...
1731 // to fill the info of a track for the relativ gain calibration
1734 Int_t nb = 0; // Nombre de zones traversees
1735 Int_t fd = -1; // Premiere zone non nulle
1736 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1738 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1740 if(nbclusters < fNumberClusters) return;
1741 if(nbclusters > fNumberClustersf) return;
1744 // Normalize with the number of clusters
1745 Double_t normalizeCst = fRelativeScale;
1746 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1748 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1750 // See if the track goes through different zones
1751 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1752 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1753 if (fAmpTotal[k] > 0.0) {
1754 totalcharge += fAmpTotal[k];
1762 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1768 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1770 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1771 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1774 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1778 if ((fAmpTotal[fd] > 0.0) &&
1779 (fAmpTotal[fd+1] > 0.0)) {
1780 // One of the two very big
1781 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1783 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1784 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1787 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1790 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1792 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1794 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1795 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1798 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1801 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1804 if (fCalibraMode->GetNfragZ(0) > 1) {
1805 if (fAmpTotal[fd] > 0.0) {
1806 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1807 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1808 // One of the two very big
1809 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1811 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1812 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1815 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1818 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1820 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1822 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1823 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1826 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1828 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1839 //____________Offine tracking in the AliTRDtracker_____________________________
1840 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1843 // For the offline tracking or mcm tracklets
1844 // This function will be called in the functions UpdateHistogram...
1845 // to fill the info of a track for the drift velocity calibration
1848 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1849 Int_t fd1 = -1; // Premiere zone non nulle
1850 Int_t fd2 = -1; // Deuxieme zone non nulle
1851 Int_t k1 = -1; // Debut de la premiere zone
1852 Int_t k2 = -1; // Debut de la seconde zone
1853 Int_t nbclusters = 0; // number of clusters
1857 // See if the track goes through different zones
1858 for (Int_t k = 0; k < fTimeMax; k++) {
1859 if (fPHValue[k] > 0.0) {
1865 if (fPHPlace[k] != fd1) {
1871 if (fPHPlace[k] != fd2) {
1878 // See if noise before and after
1879 if(fMaxCluster > 0) {
1880 if(fPHValue[0] > fMaxCluster) return;
1881 if(fTimeMax > fNbMaxCluster) {
1882 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1883 if(fPHValue[k] > fMaxCluster) return;
1888 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1890 if(nbclusters < fNumberClusters) return;
1891 if(nbclusters > fNumberClustersf) return;
1897 for (Int_t i = 0; i < fTimeMax; i++) {
1899 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1901 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1903 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1906 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1908 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1914 if ((fd1 == fd2+1) ||
1916 // One of the two fast all the think
1917 if (k2 > (k1+fDifference)) {
1918 //we choose to fill the fd1 with all the values
1920 for (Int_t i = 0; i < fTimeMax; i++) {
1922 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1924 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1928 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1930 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1935 if ((k2+fDifference) < fTimeMax) {
1936 //we choose to fill the fd2 with all the values
1938 for (Int_t i = 0; i < fTimeMax; i++) {
1940 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1942 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1946 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1948 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1954 // Two zones voisines sinon rien!
1955 if (fCalibraMode->GetNfragZ(1) > 1) {
1957 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1958 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1959 // One of the two fast all the think
1960 if (k2 > (k1+fDifference)) {
1961 //we choose to fill the fd1 with all the values
1963 for (Int_t i = 0; i < fTimeMax; i++) {
1965 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1967 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1971 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1973 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1978 if ((k2+fDifference) < fTimeMax) {
1979 //we choose to fill the fd2 with all the values
1981 for (Int_t i = 0; i < fTimeMax; i++) {
1983 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1985 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1989 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1991 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1998 // Two zones voisines sinon rien!
2000 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2001 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2002 // One of the two fast all the think
2003 if (k2 > (k1 + fDifference)) {
2004 //we choose to fill the fd1 with all the values
2006 for (Int_t i = 0; i < fTimeMax; i++) {
2008 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2010 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2014 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2016 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2021 if ((k2+fDifference) < fTimeMax) {
2022 //we choose to fill the fd2 with all the values
2024 for (Int_t i = 0; i < fTimeMax; i++) {
2026 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2028 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2032 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2034 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2046 //////////////////////////////////////////////////////////////////////////////////////////
2047 // DAQ process functions
2048 /////////////////////////////////////////////////////////////////////////////////////////
2049 //_____________________________________________________________________
2050 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2053 // Event Processing loop - AliTRDrawStream
2055 // 0 timebin problem
2058 // Same algorithm as TestBeam but different reader
2061 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2063 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2064 digitsManager->CreateArrays();
2066 Int_t withInput = 1;
2068 Double_t phvalue[16][144][36];
2069 for(Int_t k = 0; k < 36; k++){
2070 for(Int_t j = 0; j < 16; j++){
2071 for(Int_t c = 0; c < 144; c++){
2072 phvalue[j][c][k] = 0.0;
2077 fDetectorPreviousTrack = -1;
2081 Int_t nbtimebin = 0;
2082 Int_t baseline = 10;
2088 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2090 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2091 // printf("there is ADC data on this chamber!\n");
2093 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2094 if (digits->HasData()) { //array
2096 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2097 if (indexes->IsAllocated() == kFALSE) {
2098 AliError("Indexes do not exist!");
2103 indexes->ResetCounters();
2105 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2106 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2107 //while (rawStream->Next()) {
2109 Int_t idetector = det; // current detector
2110 //Int_t imcm = rawStream->GetMCM(); // current MCM
2111 //Int_t irob = rawStream->GetROB(); // current ROB
2114 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2116 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2119 for(Int_t k = 0; k < 36; k++){
2120 for(Int_t j = 0; j < 16; j++){
2121 for(Int_t c = 0; c < 144; c++){
2122 phvalue[j][c][k] = 0.0;
2128 fDetectorPreviousTrack = idetector;
2129 //fMCMPrevious = imcm;
2130 //fROBPrevious = irob;
2132 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2133 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2134 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2135 baseline = digitParam->GetADCbaseline(det); // baseline
2137 if(nbtimebin == 0) return 0;
2138 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2139 fTimeMax = nbtimebin;
2141 fNumberClustersf = fTimeMax;
2142 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2145 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2146 // phvalue[row][col][itime] = signal[itime]-baseline;
2147 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2148 /*if(phvalue[iRow][iCol][itime] >= 20) {
2149 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2153 (Short_t)(digits->GetData(iRow,iCol,itime)),
2160 // fill the last one
2161 if(fDetectorPreviousTrack != -1){
2164 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2165 // printf("\n ---> withinput %d\n\n",withInput);
2167 for(Int_t k = 0; k < 36; k++){
2168 for(Int_t j = 0; j < 16; j++){
2169 for(Int_t c = 0; c < 144; c++){
2170 phvalue[j][c][k] = 0.0;
2178 digitsManager->ClearArrays(det);
2180 delete digitsManager;
2185 //_____________________________________________________________________
2186 //////////////////////////////////////////////////////////////////////////////
2187 // Routine inside the DAQ process
2188 /////////////////////////////////////////////////////////////////////////////
2189 //_______________________________________________________________________
2190 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2193 // Look for the maximum by collapsing over the time
2194 // Sum over four pad col and two pad row
2200 Int_t idect = fDetectorPreviousTrack;
2201 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2203 for(Int_t tb = 0; tb < 36; tb++){
2207 //fGoodTracklet = kTRUE;
2208 //fDetectorPreviousTrack = 0;
2211 ///////////////////////////
2213 /////////////////////////
2217 Double_t integralMax = -1;
2219 for (Int_t ir = 1; ir <= 15; ir++)
2221 for (Int_t ic = 2; ic <= 142; ic++)
2223 Double_t integral = 0;
2224 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2226 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2228 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2229 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2232 for(Int_t tb = 0; tb< fTimeMax; tb++){
2233 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2238 if (integralMax < integral)
2242 integralMax = integral;
2244 } // check max integral
2248 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2249 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2254 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2258 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2259 //if(!fGoodTracklet) used = 1;;
2261 // /////////////////////////////////////////////////////
2262 // sum ober 2 row and 4 pad cols for each time bins
2263 // ////////////////////////////////////////////////////
2267 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2269 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2271 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2272 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2274 for(Int_t it = 0; it < fTimeMax; it++){
2275 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2282 Double_t sumcharge = 0.0;
2283 for(Int_t it = 0; it < fTimeMax; it++){
2284 sumcharge += sum[it];
2285 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2289 /////////////////////////////////////////////////////////
2291 ////////////////////////////////////////////////////////
2292 if(fDebugLevel > 1){
2293 if ( !fDebugStreamer ) {
2295 TDirectory *backup = gDirectory;
2296 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2297 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2300 Double_t amph0 = sum[0];
2301 Double_t amphlast = sum[fTimeMax-1];
2302 Double_t rms = TMath::RMS(fTimeMax,sum);
2303 Int_t goodtracklet = (Int_t) fGoodTracklet;
2304 for(Int_t it = 0; it < fTimeMax; it++){
2305 Double_t clustera = sum[it];
2307 (* fDebugStreamer) << "FillDAQa"<<
2308 "ampTotal="<<sumcharge<<
2311 "detector="<<idect<<
2313 "amphlast="<<amphlast<<
2314 "goodtracklet="<<goodtracklet<<
2315 "clustera="<<clustera<<
2323 ////////////////////////////////////////////////////////
2325 ///////////////////////////////////////////////////////
2326 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2327 if(sum[0] > 100.0) used = 1;
2328 if(nbcl < fNumberClusters) used = 1;
2329 if(nbcl > fNumberClustersf) used = 1;
2331 //if(fDetectorPreviousTrack == 15){
2332 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2334 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2336 for(Int_t it = 0; it < fTimeMax; it++){
2337 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2339 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2341 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2343 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2348 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2350 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2357 //____________Online trackling in AliTRDtrigger________________________________
2358 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2362 // Fill a simple average pulse height
2366 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2372 //____________Write_____________________________________________________
2373 //_____________________________________________________________________
2374 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2377 // Write infos to file
2381 if ( fDebugStreamer ) {
2382 delete fDebugStreamer;
2383 fDebugStreamer = 0x0;
2386 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2391 ,fNumberUsedPh[1]));
2393 TDirectory *backup = gDirectory;
2399 option = "recreate";
2401 TFile f(filename,option.Data());
2403 TStopwatch stopwatch;
2406 f.WriteTObject(fCalibraVector);
2411 f.WriteTObject(fCH2d);
2416 f.WriteTObject(fPH2d);
2421 f.WriteTObject(fPRF2d);
2424 if(fLinearFitterOn){
2425 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2426 f.WriteTObject(fLinearVdriftFit);
2431 if ( backup ) backup->cd();
2433 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2434 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2436 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2438 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2439 //___________________________________________probe the histos__________________________________________________
2440 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2443 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2444 // debug mode with 2 for TH2I and 3 for TProfile2D
2445 // It gives a pointer to a Double_t[7] with the info following...
2446 // [0] : number of calibration groups with entries
2447 // [1] : minimal number of entries found
2448 // [2] : calibration group number of the min
2449 // [3] : maximal number of entries found
2450 // [4] : calibration group number of the max
2451 // [5] : mean number of entries found
2452 // [6] : mean relative error
2455 Double_t *info = new Double_t[7];
2457 // Number of Xbins (detectors or groups of pads)
2458 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2459 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2462 Double_t nbwe = 0; //number of calibration groups with entries
2463 Double_t minentries = 0; //minimal number of entries found
2464 Double_t maxentries = 0; //maximal number of entries found
2465 Double_t placemin = 0; //calibration group number of the min
2466 Double_t placemax = -1; //calibration group number of the max
2467 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2468 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2470 Double_t counter = 0;
2473 TH1F *nbEntries = 0x0;//distribution of the number of entries
2474 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2475 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2477 // Beginning of the loop over the calibration groups
2478 for (Int_t idect = 0; idect < nbins; idect++) {
2480 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2481 projch->SetDirectory(0);
2483 // Number of entries for this calibration group
2484 Double_t nentries = 0.0;
2486 for (Int_t k = 0; k < nxbins; k++) {
2487 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2491 for (Int_t k = 0; k < nxbins; k++) {
2492 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2493 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2494 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2503 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2504 nbEntries->SetDirectory(0);
2505 nbEntries->Fill(nentries);
2506 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2507 nbEntriesPerGroup->SetDirectory(0);
2508 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2509 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));
2510 nbEntriesPerSp->SetDirectory(0);
2511 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2516 if(nentries > maxentries){
2517 maxentries = nentries;
2521 minentries = nentries;
2523 if(nentries < minentries){
2524 minentries = nentries;
2530 meanstats += nentries;
2532 }//calibration groups loop
2534 if(nbwe > 0) meanstats /= nbwe;
2535 if(counter > 0) meanrelativerror /= counter;
2537 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2538 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2539 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2540 AliInfo(Form("The mean number of entries is %f",meanstats));
2541 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2544 info[1] = minentries;
2546 info[3] = maxentries;
2548 info[5] = meanstats;
2549 info[6] = meanrelativerror;
2551 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2552 gStyle->SetPalette(1);
2553 gStyle->SetOptStat(1111);
2554 gStyle->SetPadBorderMode(0);
2555 gStyle->SetCanvasColor(10);
2556 gStyle->SetPadLeftMargin(0.13);
2557 gStyle->SetPadRightMargin(0.01);
2558 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2561 nbEntries->Draw("");
2563 nbEntriesPerSp->SetStats(0);
2564 nbEntriesPerSp->Draw("");
2565 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2567 nbEntriesPerGroup->SetStats(0);
2568 nbEntriesPerGroup->Draw("");
2574 //____________________________________________________________________________
2575 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2578 // Return a Int_t[4] with:
2579 // 0 Mean number of entries
2580 // 1 median of number of entries
2581 // 2 rms of number of entries
2582 // 3 number of group with entries
2585 Double_t *stat = new Double_t[4];
2588 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2590 Double_t *weight = new Double_t[nbofgroups];
2591 Double_t *nonul = new Double_t[nbofgroups];
2593 for(Int_t k = 0; k < nbofgroups; k++){
2594 if(fEntriesCH[k] > 0) {
2596 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2599 else weight[k] = 0.0;
2601 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2602 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2603 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2611 //____________________________________________________________________________
2612 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2615 // Return a Int_t[4] with:
2616 // 0 Mean number of entries
2617 // 1 median of number of entries
2618 // 2 rms of number of entries
2619 // 3 number of group with entries
2622 Double_t *stat = new Double_t[4];
2625 Int_t nbofgroups = 540;
2626 Double_t *weight = new Double_t[nbofgroups];
2627 Int_t *nonul = new Int_t[nbofgroups];
2629 for(Int_t k = 0; k < nbofgroups; k++){
2630 if(fEntriesLinearFitter[k] > 0) {
2632 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2635 else weight[k] = 0.0;
2637 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2638 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2639 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2647 //////////////////////////////////////////////////////////////////////////////////////
2649 //////////////////////////////////////////////////////////////////////////////////////
2650 //_____________________________________________________________________________
2651 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2654 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2655 // If fNgroupprf is zero then no binning in tnp
2659 name += fCalibraMode->GetNz(2);
2661 name += fCalibraMode->GetNrphi(2);
2665 if(fNgroupprf != 0){
2667 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2668 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2669 fPRF2d->SetYTitle("Det/pad groups");
2670 fPRF2d->SetXTitle("Position x/W [pad width units]");
2671 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2672 fPRF2d->SetStats(0);
2675 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2676 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2677 fPRF2d->SetYTitle("Det/pad groups");
2678 fPRF2d->SetXTitle("Position x/W [pad width units]");
2679 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2680 fPRF2d->SetStats(0);
2685 //_____________________________________________________________________________
2686 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2689 // Create the 2D histos
2692 TString name("Ver");
2693 name += fVersionVdriftUsed;
2695 name += fSubVersionVdriftUsed;
2697 name += fFirstRunVdrift;
2699 name += fCalibraMode->GetNz(1);
2701 name += fCalibraMode->GetNrphi(1);
2703 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2704 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2706 fPH2d->SetYTitle("Det/pad groups");
2707 fPH2d->SetXTitle("time [#mus]");
2708 fPH2d->SetZTitle("<PH> [a.u.]");
2712 //_____________________________________________________________________________
2713 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2716 // Create the 2D histos
2719 TString name("Ver");
2720 name += fVersionGainUsed;
2722 name += fSubVersionGainUsed;
2724 name += fFirstRunGain;
2726 name += fCalibraMode->GetNz(0);
2728 name += fCalibraMode->GetNrphi(0);
2730 fCH2d = new TH2I("CH2d",(const Char_t *) name
2731 ,fNumberBinCharge,0,300,nn,0,nn);
2732 fCH2d->SetYTitle("Det/pad groups");
2733 fCH2d->SetXTitle("charge deposit [a.u]");
2734 fCH2d->SetZTitle("counts");
2739 //////////////////////////////////////////////////////////////////////////////////
2740 // Set relative scale
2741 /////////////////////////////////////////////////////////////////////////////////
2742 //_____________________________________________________________________________
2743 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2746 // Set the factor that will divide the deposited charge
2747 // to fit in the histo range [0,300]
2750 if (RelativeScale > 0.0) {
2751 fRelativeScale = RelativeScale;
2754 AliInfo("RelativeScale must be strict positif!");
2758 //////////////////////////////////////////////////////////////////////////////////
2759 // Quick way to fill a histo
2760 //////////////////////////////////////////////////////////////////////////////////
2761 //_____________________________________________________________________
2762 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2765 // FillCH2d: Marian style
2768 //skip simply the value out of range
2769 if((y>=300.0) || (y<0.0)) return;
2771 //Calcul the y place
2772 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2773 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2776 fCH2d->GetArray()[place]++;
2780 //////////////////////////////////////////////////////////////////////////////////
2781 // Geometrical functions
2782 ///////////////////////////////////////////////////////////////////////////////////
2783 //_____________________________________________________________________________
2784 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2787 // Reconstruct the layer number from the detector number
2790 return ((Int_t) (d % 6));
2794 //_____________________________________________________________________________
2795 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2798 // Reconstruct the stack number from the detector number
2800 const Int_t kNlayer = 6;
2802 return ((Int_t) (d % 30) / kNlayer);
2806 //_____________________________________________________________________________
2807 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2810 // Reconstruct the sector number from the detector number
2814 return ((Int_t) (d / fg));
2817 ///////////////////////////////////////////////////////////////////////////////////
2818 // Getter functions for DAQ of the CH2d and the PH2d
2819 //////////////////////////////////////////////////////////////////////////////////
2820 //_____________________________________________________________________
2821 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2824 // return pointer to fPH2d TProfile2D
2825 // create a new TProfile2D if it doesn't exist allready
2831 fTimeMax = nbtimebin;
2832 fSf = samplefrequency;
2838 //_____________________________________________________________________
2839 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2842 // return pointer to fCH2d TH2I
2843 // create a new TH2I if it doesn't exist allready
2852 ////////////////////////////////////////////////////////////////////////////////////////////
2853 // Drift velocity calibration
2854 ///////////////////////////////////////////////////////////////////////////////////////////
2855 //_____________________________________________________________________
2856 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2859 // return pointer to TLinearFitter Calibration
2860 // if force is true create a new TLinearFitter if it doesn't exist allready
2863 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2864 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2867 // if we are forced and TLinearFitter doesn't yet exist create it
2869 // new TLinearFitter
2870 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2871 fLinearFitterArray.AddAt(linearfitter,detector);
2872 return linearfitter;
2875 //____________________________________________________________________________
2876 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2879 // Analyse array of linear fitter because can not be written
2880 // Store two arrays: one with the param the other one with the error param + number of entries
2883 for(Int_t k = 0; k < 540; k++){
2884 TLinearFitter *linearfitter = GetLinearFitter(k);
2885 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2886 TVectorD *par = new TVectorD(2);
2887 TVectorD pare = TVectorD(2);
2888 TVectorD *parE = new TVectorD(3);
2889 if((linearfitter->EvalRobust(0.8)==0)) {
2890 //linearfitter->Eval();
2891 linearfitter->GetParameters(*par);
2892 //linearfitter->GetErrors(pare);
2893 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2894 //(*parE)[0] = pare[0]*ppointError;
2895 //(*parE)[1] = pare[1]*ppointError;
2899 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2900 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2901 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);