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)
156 ,fTakeSnapshot(kTRUE)
159 ,fSubVersionGainUsed(0)
160 ,fFirstRunGainLocal(0)
161 ,fVersionGainLocalUsed(0)
162 ,fSubVersionGainLocalUsed(0)
164 ,fVersionVdriftUsed(0)
165 ,fSubVersionVdriftUsed(0)
168 ,fSubVersionExBUsed(0)
169 ,fCalibraMode(new AliTRDCalibraMode())
172 ,fDetectorPreviousTrack(-1)
176 ,fNumberClustersf(30)
177 ,fNumberClustersProcent(0.5)
178 ,fThresholdClustersDAQ(120.0)
186 ,fRangeHistoCharge(150)
187 ,fNumberBinCharge(50)
193 ,fGoodTracklet(kTRUE)
194 ,fLinearFitterTracklet(0x0)
196 ,fEntriesLinearFitter(0x0)
201 ,fLinearFitterArray(540)
202 ,fLinearVdriftFit(0x0)
208 // Default constructor
212 // Init some default values
215 fNumberUsedCh[0] = 0;
216 fNumberUsedCh[1] = 0;
217 fNumberUsedPh[0] = 0;
218 fNumberUsedPh[1] = 0;
220 fGeo = new AliTRDgeometry();
221 fCalibDB = AliTRDcalibDB::Instance();
224 //______________________________________________________________________________________
225 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
232 ,fPRF2dOn(c.fPRF2dOn)
233 ,fHisto2d(c.fHisto2d)
234 ,fVector2d(c.fVector2d)
235 ,fLinearFitterOn(c.fLinearFitterOn)
236 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
237 ,fExbAltFitOn(c.fExbAltFitOn)
238 ,fScaleWithTPCSignal(c.fScaleWithTPCSignal)
239 ,fRelativeScale(c.fRelativeScale)
240 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
241 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
242 ,fFillWithZero(c.fFillWithZero)
243 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
244 ,fMaxCluster(c.fMaxCluster)
245 ,fNbMaxCluster(c.fNbMaxCluster)
246 ,fCutWithVdriftCalib(c.fCutWithVdriftCalib)
247 ,fMinNbTRDtracklets(c.fMinNbTRDtracklets)
248 ,fMinTRDMomentum(c.fMinTRDMomentum)
249 ,fTakeSnapshot(c.fTakeSnapshot)
250 ,fFirstRunGain(c.fFirstRunGain)
251 ,fVersionGainUsed(c.fVersionGainUsed)
252 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
253 ,fFirstRunGainLocal(c.fFirstRunGainLocal)
254 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
255 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
256 ,fFirstRunVdrift(c.fFirstRunVdrift)
257 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
258 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
259 ,fFirstRunExB(c.fFirstRunExB)
260 ,fVersionExBUsed(c.fVersionExBUsed)
261 ,fSubVersionExBUsed(c.fSubVersionExBUsed)
264 ,fDebugLevel(c.fDebugLevel)
265 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
266 ,fMCMPrevious(c.fMCMPrevious)
267 ,fROBPrevious(c.fROBPrevious)
268 ,fNumberClusters(c.fNumberClusters)
269 ,fNumberClustersf(c.fNumberClustersf)
270 ,fNumberClustersProcent(c.fNumberClustersProcent)
271 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
272 ,fNumberRowDAQ(c.fNumberRowDAQ)
273 ,fNumberColDAQ(c.fNumberColDAQ)
274 ,fProcent(c.fProcent)
275 ,fDifference(c.fDifference)
276 ,fNumberTrack(c.fNumberTrack)
277 ,fTimeMax(c.fTimeMax)
279 ,fRangeHistoCharge(c.fRangeHistoCharge)
280 ,fNumberBinCharge(c.fNumberBinCharge)
281 ,fNumberBinPRF(c.fNumberBinPRF)
282 ,fNgroupprf(c.fNgroupprf)
286 ,fGoodTracklet(c.fGoodTracklet)
287 ,fLinearFitterTracklet(0x0)
289 ,fEntriesLinearFitter(0x0)
294 ,fLinearFitterArray(540)
295 ,fLinearVdriftFit(0x0)
303 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
304 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
306 fPH2d = new TProfile2D(*c.fPH2d);
307 fPH2d->SetDirectory(0);
310 fPRF2d = new TProfile2D(*c.fPRF2d);
311 fPRF2d->SetDirectory(0);
314 fCH2d = new TH2I(*c.fCH2d);
315 fCH2d->SetDirectory(0);
317 if(c.fLinearVdriftFit){
318 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
321 fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
324 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
325 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
330 fGeo = new AliTRDgeometry();
331 fCalibDB = AliTRDcalibDB::Instance();
333 fNumberUsedCh[0] = 0;
334 fNumberUsedCh[1] = 0;
335 fNumberUsedPh[0] = 0;
336 fNumberUsedPh[1] = 0;
340 //____________________________________________________________________________________
341 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
344 // AliTRDCalibraFillHisto destructor
348 if ( fDebugStreamer ) delete fDebugStreamer;
350 if ( fCalDetGain ) delete fCalDetGain;
351 if ( fCalROCGain ) delete fCalROCGain;
353 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
357 delete [] fEntriesCH;
358 delete [] fEntriesLinearFitter;
361 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
362 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
365 if(fLinearVdriftFit) delete fLinearVdriftFit;
366 if(fExbAltFit) delete fExbAltFit;
372 //_____________________________________________________________________________
373 void AliTRDCalibraFillHisto::Destroy()
384 //_____________________________________________________________________________
385 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
388 // Delete DebugStreamer
391 if ( fDebugStreamer ) delete fDebugStreamer;
392 fDebugStreamer = 0x0;
395 //_____________________________________________________________________________
396 void AliTRDCalibraFillHisto::ClearHistos()
416 //////////////////////////////////////////////////////////////////////////////////
417 // calibration with AliTRDtrackV1: Init, Update
418 //////////////////////////////////////////////////////////////////////////////////
419 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
420 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
423 // Init the histograms and stuff to be filled
428 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
430 AliInfo("Could not get calibDB");
433 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
435 AliInfo("Could not get CommonParam");
440 if(nboftimebin > 0) fTimeMax = nboftimebin;
441 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
442 if(fTimeMax <= 0) fTimeMax = 30;
443 printf("////////////////////////////////////////////\n");
444 printf("Number of time bins in calibration component %d\n",fTimeMax);
445 printf("////////////////////////////////////////////\n");
446 fSf = parCom->GetSamplingFrequency();
447 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
448 else fRelativeScale = 1.18;
449 fNumberClustersf = fTimeMax;
450 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
452 // Init linear fitter
453 if(!fLinearFitterTracklet) {
454 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
455 fLinearFitterTracklet->StoreData(kTRUE);
458 // Calcul Xbins Chambd0, Chamb2
459 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
460 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
461 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
463 // If vector method On initialised all the stuff
465 fCalibraVector = new AliTRDCalibraVector();
466 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
467 fCalibraVector->SetTimeMax(fTimeMax);
468 if(fNgroupprf != 0) {
469 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
470 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
473 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
474 fCalibraVector->SetPRFRange(1.5);
476 for(Int_t k = 0; k < 3; k++){
477 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
478 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
480 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
481 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
482 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
483 fCalibraVector->SetNbGroupPRF(fNgroupprf);
486 // Create the 2D histos corresponding to the pad groupCalibration mode
489 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
490 ,fCalibraMode->GetNz(0)
491 ,fCalibraMode->GetNrphi(0)));
493 // Create the 2D histo
498 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
499 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
503 fEntriesCH = new Int_t[ntotal0];
504 for(Int_t k = 0; k < ntotal0; k++){
511 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
512 ,fCalibraMode->GetNz(1)
513 ,fCalibraMode->GetNrphi(1)));
515 // Create the 2D histo
520 fPHPlace = new Short_t[fTimeMax];
521 for (Int_t k = 0; k < fTimeMax; k++) {
524 fPHValue = new Float_t[fTimeMax];
525 for (Int_t k = 0; k < fTimeMax; k++) {
529 if (fLinearFitterOn) {
530 if(fLinearFitterDebugOn) {
531 fLinearFitterArray.SetName("ArrayLinearFitters");
532 fEntriesLinearFitter = new Int_t[540];
533 for(Int_t k = 0; k < 540; k++){
534 fEntriesLinearFitter[k] = 0;
537 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
538 TString nameee("Ver");
539 nameee += fVersionExBUsed;
541 nameee += fSubVersionExBUsed;
542 nameee += "FirstRun";
543 nameee += fFirstRunExB;
545 fLinearVdriftFit->SetNameCalibUsed(nameee);
548 fExbAltFit = new AliTRDCalibraExbAltFit();
553 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
554 ,fCalibraMode->GetNz(2)
555 ,fCalibraMode->GetNrphi(2)));
556 // Create the 2D histo
558 CreatePRF2d(ntotal2);
565 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
566 Bool_t AliTRDCalibraFillHisto::InitCalDet()
569 // Init the Gain Cal Det
574 AliCDBEntry *entry = 0x0;
576 entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor");
579 entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
582 AliError("No gain det calibration entry found");
585 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
587 AliError("No calDet gain found");
593 fCalDetGain->~AliTRDCalDet();
594 new(fCalDetGain) AliTRDCalDet(*(calDet));
595 }else fCalDetGain = new AliTRDCalDet(*(calDet));
600 name += fVersionGainUsed;
602 name += fSubVersionGainUsed;
604 name += fFirstRunGain;
606 name += fCalibraMode->GetNz(0);
608 name += fCalibraMode->GetNrphi(0);
610 fCH2d->SetTitle(name);
613 TString namee("Ver");
614 namee += fVersionVdriftUsed;
616 namee += fSubVersionVdriftUsed;
618 namee += fFirstRunVdrift;
620 namee += fCalibraMode->GetNz(1);
622 namee += fCalibraMode->GetNrphi(1);
624 fPH2d->SetTitle(namee);
626 // title AliTRDCalibraVdriftLinearFit
627 TString nameee("Ver");
628 nameee += fVersionExBUsed;
630 nameee += fSubVersionExBUsed;
631 nameee += "FirstRun";
632 nameee += fFirstRunExB;
636 fLinearVdriftFit->SetNameCalibUsed(nameee);
643 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
644 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
647 // Init the Gain Cal Pad
652 AliCDBEntry *entry = 0x0;
654 entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor");
657 entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
660 AliError("No gain pad calibration entry found");
663 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
665 AliError("No calPad gain found");
668 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
670 AliError("No calRoc gain found");
675 fCalROCGain->~AliTRDCalROC();
676 new(fCalROCGain) AliTRDCalROC(*(calRoc));
677 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
686 //____________Offline tracking in the AliTRDtracker____________________________
687 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t,const AliESDtrack *esdtrack)
690 // Use AliTRDtrackV1 for the calibration
694 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
695 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
696 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
697 Bool_t newtr = kTRUE; // new track
701 // Cut on the number of TRD tracklets
703 Int_t numberoftrdtracklets = t->GetNumberOfTracklets();
704 if(numberoftrdtracklets < fMinNbTRDtracklets) return kFALSE;
706 Double_t tpcsignal = 1.0;
707 if(esdtrack) tpcsignal = esdtrack->GetTPCsignal()/50.0;
708 if(fScaleWithTPCSignal && tpcsignal <0.00001) return kFALSE;
712 AliInfo("Could not get calibDB");
717 ///////////////////////////
718 // loop over the tracklet
719 ///////////////////////////
720 for(Int_t itr = 0; itr < 6; itr++){
722 if(!(tracklet = t->GetTracklet(itr))) continue;
723 if(!tracklet->IsOK()) continue;
725 ResetfVariablestracklet();
726 Float_t momentum = t->GetMomentum(itr);
727 if(TMath::Abs(momentum) < fMinTRDMomentum) continue;
730 //////////////////////////////////////////
731 // localisation of the tracklet and dqdl
732 //////////////////////////////////////////
733 Int_t layer = tracklet->GetPlane();
735 while(!(cl = tracklet->GetClusters(ic++))) continue;
736 Int_t detector = cl->GetDetector();
737 if (detector != fDetectorPreviousTrack) {
738 // if not a new track
740 // don't use the rest of this track if in the same plane
741 if (layer == GetLayer(fDetectorPreviousTrack)) {
742 //printf("bad tracklet, same layer for detector %d\n",detector);
746 //Localise the detector bin
747 LocalisationDetectorXbins(detector);
749 if(!fIsHLT) InitCalPad(detector);
752 fDetectorPreviousTrack = detector;
756 ////////////////////////////
757 // loop over the clusters
758 ////////////////////////////
759 Double_t chargeQ = 0.0;
760 Int_t nbclusters = 0;
761 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
762 if(!(cl = tracklet->GetClusters(jc))) continue;
765 // Store the info bis of the tracklet
766 Int_t row = cl->GetPadRow();
767 Int_t col = cl->GetPadCol();
768 CheckGoodTrackletV1(cl);
769 Int_t group[2] = {0,0};
770 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
771 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
772 // Add the charge if shared cluster
773 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
775 //Scale with TPC signal or not
776 if(!fScaleWithTPCSignal) {
777 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
778 //printf("Do not scale now\n");
781 chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc)/tpcsignal,group,row,col,cls); //tracklet->GetdQdl(jc)
786 ////////////////////////////////////////
787 // Fill the stuffs if a good tracklet
788 ////////////////////////////////////////
791 // drift velocity unables to cut bad tracklets
792 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
794 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
798 if(fCutWithVdriftCalib) {
799 if(pass) FillTheInfoOfTheTrackCH(nbclusters);
801 FillTheInfoOfTheTrackCH(nbclusters);
807 if(fCutWithVdriftCalib) {
808 if(pass) FillTheInfoOfTheTrackPH();
811 FillTheInfoOfTheTrackPH();
815 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
818 /////////////////////////////////////////////////////////
820 ////////////////////////////////////////////////////////
823 if ( !fDebugStreamer ) {
825 TDirectory *backup = gDirectory;
826 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
827 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
830 Int_t stacke = AliTRDgeometry::GetStack(detector);
831 Int_t sme = AliTRDgeometry::GetSector(detector);
832 Int_t layere = AliTRDgeometry::GetLayer(detector);
834 Float_t b[2] = {0.0,0.0};
835 Float_t bCov[3] = {0.0,0.0,0.0};
836 if(esdtrack) esdtrack->GetImpactParameters(b,bCov);
837 if (bCov[0]<=0 || bCov[2]<=0) {
838 bCov[0]=0; bCov[2]=0;
840 Float_t dcaxy = b[0];
842 Int_t tpcnbclusters = 0;
843 if(esdtrack) tpcnbclusters = esdtrack->GetTPCclusters(0);
844 Double_t ttpcsignal = 0.0;
845 if(esdtrack) ttpcsignal = esdtrack->GetTPCsignal();
846 Int_t cutvdriftlinear = 0;
847 if(!pass) cutvdriftlinear = 1;
849 (* fDebugStreamer) << "FillCharge"<<
850 "detector="<<detector<<
856 "nbtpccls="<<tpcnbclusters<<
857 "tpcsignal="<<ttpcsignal<<
858 "cutvdriftlinear="<<cutvdriftlinear<<
860 "nbtrdtracklet="<<numberoftrdtracklets<<
866 } // if a good tracklet
872 ///////////////////////////////////////////////////////////////////////////////////
873 // Routine inside the update with AliTRDtrack
874 ///////////////////////////////////////////////////////////////////////////////////
875 //____________Offine tracking in the AliTRDtracker_____________________________
876 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
879 // Drift velocity calibration:
880 // Fit the clusters with a straight line
881 // From the slope find the drift velocity
884 ////////////////////////////////////////////////
885 //Number of points: if less than 3 return kFALSE
886 /////////////////////////////////////////////////
887 if(nbclusters <= 2) return kFALSE;
892 // results of the linear fit
893 Double_t dydt = 0.0; // dydt tracklet after straight line fit
894 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
895 Double_t pointError = 0.0; // error after straight line fit
896 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
897 Int_t crossrow = 0; // if it crosses a pad row
898 Int_t rowp = -1; // if it crosses a pad row
899 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
900 fLinearFitterTracklet->ClearPoints();
903 ///////////////////////////////////////////
904 // Take the parameters of the track
905 //////////////////////////////////////////
906 // take now the snp, tnp and tgl from the track
907 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
908 Double_t tnp = 0.0; // dy/dx at the end of the chamber
909 if( TMath::Abs(snp) < 1.){
910 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
912 Double_t tgl = tracklet->GetTgl(); // dz/dl
913 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
915 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
916 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
917 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
918 // at the end with correction due to linear fit
919 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
920 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
923 ////////////////////////////
924 // loop over the clusters
925 ////////////////////////////
927 AliTRDcluster *cl = 0x0;
928 //////////////////////////////
929 // Check no shared clusters
930 //////////////////////////////
931 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
932 cl = tracklet->GetClusters(icc);
935 //////////////////////////////////
937 //////////////////////////////////
939 Float_t sigArr[AliTRDfeeParam::GetNcol()];
940 memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
941 Int_t ncl=0, tbf=0, tbl=0;
943 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
944 if(!(cl = tracklet->GetClusters(ic))) continue;
949 Int_t col = cl->GetPadCol();
950 for(int ip=-1, jp=2; jp<5; ip++, jp++){
952 if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
953 sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
956 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
958 Double_t ycluster = cl->GetY();
959 Int_t time = cl->GetPadTime();
960 Double_t timeis = time/fSf;
961 //See if cross two pad rows
962 Int_t row = cl->GetPadRow();
963 if(rowp==-1) rowp = row;
964 if(row != rowp) crossrow = 1;
966 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
972 ////////////////////////////////////
973 // Do the straight line fit now
974 ///////////////////////////////////
976 fLinearFitterTracklet->ClearPoints();
980 fLinearFitterTracklet->Eval();
981 fLinearFitterTracklet->GetParameters(pars);
982 pointError = TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
983 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
985 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
986 fLinearFitterTracklet->ClearPoints();
988 ////////////////////////////////////
989 // Calc the projection of the clusters on the y direction
990 ///////////////////////////////////
992 Float_t signalSum(0.);
993 Float_t mean = 0.0, rms = 0.0;
994 Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
995 Float_t dz = dzdx*(tbl-tbf)/10;
997 for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
998 signalSum+=sigArr[ip];
1001 if(signalSum > 0.0) mean/=signalSum;
1003 for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
1004 rms+=sigArr[ip]*(ip-mean)*(ip-mean);
1006 if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
1008 rms -= TMath::Abs(dz*tilt);
1012 ////////////////////////////////
1014 ///////////////////////////////
1017 if(fDebugLevel > 1){
1018 if ( !fDebugStreamer ) {
1020 TDirectory *backup = gDirectory;
1021 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1022 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1025 float xcoord = tnp-dzdx*tnt;
1026 float pt = tracklet->GetPt();
1027 Int_t layer = GetLayer(fDetectorPreviousTrack);
1029 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1030 //"snpright="<<snpright<<
1032 "nbclusters="<<nbclusters<<
1033 "detector="<<fDetectorPreviousTrack<<
1041 "crossrow="<<crossrow<<
1042 "errorpar="<<errorpar<<
1043 "pointError="<<pointError<<
1055 /////////////////////////
1057 ////////////////////////
1059 if(nbclusters < fNumberClusters) return kFALSE;
1060 if(nbclusters > fNumberClustersf) return kFALSE;
1061 if(pointError >= 0.3) return kFALSE;
1062 if(crossrow == 1) return kTRUE;
1064 ///////////////////////
1066 //////////////////////
1068 if(fLinearFitterOn){
1069 //Add to the linear fitter of the detector
1070 if( TMath::Abs(snp) < 1.){
1071 Double_t x = tnp-dzdx*tnt;
1072 if(fLinearFitterDebugOn) {
1073 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1074 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1076 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1080 fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
1085 //____________Offine tracking in the AliTRDtracker_____________________________
1086 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1089 // PRF width calibration
1090 // Assume a Gaussian shape: determinate the position of the three pad clusters
1091 // Fit with a straight line
1092 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1093 // Fill the PRF as function of angle of the track
1097 //printf("begin\n");
1098 ///////////////////////////////////////////
1099 // Take the parameters of the track
1100 //////////////////////////////////////////
1101 // take now the snp, tnp and tgl from the track
1102 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1103 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1104 if( TMath::Abs(snp) < 1.){
1105 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1107 Double_t tgl = tracklet->GetTgl(); // dz/dl
1108 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1110 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1111 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1112 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1113 // at the end with correction due to linear fit
1114 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1115 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1117 ///////////////////////////////
1118 // Calculate tnp group shift
1119 ///////////////////////////////
1120 Bool_t echec = kFALSE;
1121 Double_t shift = 0.0;
1122 //Calculate the shift in x coresponding to this tnp
1123 if(fNgroupprf != 0.0){
1124 shift = -3.0*(fNgroupprf-1)-1.5;
1125 Double_t limithigh = -0.2*(fNgroupprf-1);
1126 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1128 while(tnp > limithigh){
1134 // do nothing if out of tnp range
1135 //printf("echec %d\n",(Int_t)echec);
1136 if(echec) return kFALSE;
1138 ///////////////////////
1140 //////////////////////
1142 Int_t nb3pc = 0; // number of three pads clusters used for fit
1143 // to see the difference between the fit and the 3 pad clusters position
1144 Double_t padPositions[AliTRDseedV1::kNtb];
1145 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1146 fLinearFitterTracklet->ClearPoints();
1148 //printf("loop clusters \n");
1149 ////////////////////////////
1150 // loop over the clusters
1151 ////////////////////////////
1152 AliTRDcluster *cl = 0x0;
1153 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1154 // reject shared clusters on pad row
1155 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1156 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1159 cl = tracklet->GetClusters(ic);
1161 Double_t time = cl->GetPadTime();
1162 if((time<=7) || (time>=21)) continue;
1163 Short_t *signals = cl->GetSignals();
1164 Float_t xcenter = 0.0;
1165 Bool_t echec1 = kTRUE;
1167 /////////////////////////////////////////////////////////////
1168 // Center 3 balanced: position with the center of the pad
1169 /////////////////////////////////////////////////////////////
1170 if ((((Float_t) signals[3]) > 0.0) &&
1171 (((Float_t) signals[2]) > 0.0) &&
1172 (((Float_t) signals[4]) > 0.0)) {
1174 // Security if the denomiateur is 0
1175 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1176 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1177 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1178 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1179 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1185 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1186 if(echec1) continue;
1188 ////////////////////////////////////////////////////////
1189 //if no echec1: calculate with the position of the pad
1190 // Position of the cluster
1191 // fill the linear fitter
1192 ///////////////////////////////////////////////////////
1193 Double_t padPosition = xcenter + cl->GetPadCol();
1194 padPositions[ic] = padPosition;
1196 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1201 //printf("Fin loop clusters \n");
1202 //////////////////////////////
1203 // fit with a straight line
1204 /////////////////////////////
1206 fLinearFitterTracklet->ClearPoints();
1209 fLinearFitterTracklet->Eval();
1211 fLinearFitterTracklet->GetParameters(line);
1212 Float_t pointError = -1.0;
1213 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1214 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1216 fLinearFitterTracklet->ClearPoints();
1218 //printf("PRF second loop \n");
1219 ////////////////////////////////////////////////
1220 // Fill the PRF: Second loop over clusters
1221 //////////////////////////////////////////////
1222 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1223 // reject shared clusters on pad row
1224 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1225 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1227 cl = tracklet->GetClusters(ic);
1230 Short_t *signals = cl->GetSignals(); // signal
1231 Double_t time = cl->GetPadTime(); // time bin
1232 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1233 Float_t padPos = cl->GetPadCol(); // middle pad
1234 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1235 Float_t ycenter = 0.0; // relative center charge
1236 Float_t ymin = 0.0; // relative left charge
1237 Float_t ymax = 0.0; // relative right charge
1239 ////////////////////////////////////////////////////////////////
1240 // Calculate ycenter, ymin and ymax even for two pad clusters
1241 ////////////////////////////////////////////////////////////////
1242 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1243 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1244 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1245 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1246 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1247 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1250 /////////////////////////
1251 // Calibration group
1252 ////////////////////////
1253 Int_t rowcl = cl->GetPadRow(); // row of cluster
1254 Int_t colcl = cl->GetPadCol(); // col of cluster
1255 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1256 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1257 Float_t xcl = cl->GetY(); // y cluster
1258 Float_t qcl = cl->GetQ(); // charge cluster
1259 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1260 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1261 Double_t xdiff = dpad; // reconstructed position constant
1262 Double_t x = dpad; // reconstructed position moved
1263 Float_t ep = pointError; // error of fit
1264 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1265 Float_t signal3 = (Float_t)signals[3]; // signal
1266 Float_t signal2 = (Float_t)signals[2]; // signal
1267 Float_t signal4 = (Float_t)signals[4]; // signal
1268 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1272 /////////////////////
1274 ////////////////////
1276 if(fDebugLevel > 1){
1277 if ( !fDebugStreamer ) {
1279 TDirectory *backup = gDirectory;
1280 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1281 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1286 Float_t y = ycenter;
1287 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1288 "caligroup="<<caligroup<<
1289 "detector="<<fDetectorPreviousTrack<<
1292 "npoints="<<nbclusters<<
1301 "padPosition="<<padPositions[ic]<<
1302 "padPosTracklet="<<padPosTracklet<<
1307 "signal1="<<signal1<<
1308 "signal2="<<signal2<<
1309 "signal3="<<signal3<<
1310 "signal4="<<signal4<<
1311 "signal5="<<signal5<<
1317 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1318 "caligroup="<<caligroup<<
1319 "detector="<<fDetectorPreviousTrack<<
1322 "npoints="<<nbclusters<<
1331 "padPosition="<<padPositions[ic]<<
1332 "padPosTracklet="<<padPosTracklet<<
1337 "signal1="<<signal1<<
1338 "signal2="<<signal2<<
1339 "signal3="<<signal3<<
1340 "signal4="<<signal4<<
1341 "signal5="<<signal5<<
1347 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1348 "caligroup="<<caligroup<<
1349 "detector="<<fDetectorPreviousTrack<<
1352 "npoints="<<nbclusters<<
1361 "padPosition="<<padPositions[ic]<<
1362 "padPosTracklet="<<padPosTracklet<<
1367 "signal1="<<signal1<<
1368 "signal2="<<signal2<<
1369 "signal3="<<signal3<<
1370 "signal4="<<signal4<<
1371 "signal5="<<signal5<<
1377 /////////////////////
1379 /////////////////////
1380 if(nbclusters < fNumberClusters) continue;
1381 if(nbclusters > fNumberClustersf) continue;
1382 if(nb3pc <= 5) continue;
1383 if((time >= 21) || (time < 7)) continue;
1384 if(TMath::Abs(qcl) < 80) continue;
1385 if( TMath::Abs(snp) > 1.) continue;
1388 ////////////////////////
1390 ///////////////////////
1392 if(TMath::Abs(dpad) < 1.5) {
1393 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1394 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1395 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1397 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1398 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1399 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1401 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1402 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1403 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1408 if(TMath::Abs(dpad) < 1.5) {
1409 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1410 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1412 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1413 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1414 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1416 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1417 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1418 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1421 } // second loop over clusters
1426 ///////////////////////////////////////////////////////////////////////////////////////
1427 // Pad row col stuff: see if masked or not
1428 ///////////////////////////////////////////////////////////////////////////////////////
1429 //_____________________________________________________________________________
1430 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1433 // See if we are not near a masked pad
1436 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1440 //_____________________________________________________________________________
1441 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1444 // See if we are not near a masked pad
1447 if (!IsPadOn(detector, col, row)) {
1448 fGoodTracklet = kFALSE;
1452 if (!IsPadOn(detector, col-1, row)) {
1453 fGoodTracklet = kFALSE;
1458 if (!IsPadOn(detector, col+1, row)) {
1459 fGoodTracklet = kFALSE;
1464 //_____________________________________________________________________________
1465 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1468 // Look in the choosen database if the pad is On.
1469 // If no the track will be "not good"
1472 // Get the parameter object
1473 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1475 AliInfo("Could not get calibDB");
1479 if (!cal->IsChamberGood(detector) ||
1480 cal->IsChamberNoData(detector) ||
1481 cal->IsPadMasked(detector,col,row)) {
1489 ///////////////////////////////////////////////////////////////////////////////////////
1490 // Calibration groups: calculate the number of groups, localise...
1491 ////////////////////////////////////////////////////////////////////////////////////////
1492 //_____________________________________________________________________________
1493 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1496 // Calculate the calibration group number for i
1499 // Row of the cluster and position in the pad groups
1501 if (fCalibraMode->GetNnZ(i) != 0) {
1502 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1506 // Col of the cluster and position in the pad groups
1508 if (fCalibraMode->GetNnRphi(i) != 0) {
1509 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1512 return posc*fCalibraMode->GetNfragZ(i)+posr;
1515 //____________________________________________________________________________________
1516 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1519 // Calculate the total number of calibration groups
1525 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1527 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1532 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1534 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1539 fCalibraMode->ModePadCalibration(2,i);
1540 fCalibraMode->ModePadFragmentation(0,2,0,i);
1541 fCalibraMode->SetDetChamb2(i);
1542 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1543 fCalibraMode->ModePadCalibration(0,i);
1544 fCalibraMode->ModePadFragmentation(0,0,0,i);
1545 fCalibraMode->SetDetChamb0(i);
1546 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1547 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1551 //_____________________________________________________________________________
1552 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1555 // Set the mode of calibration group in the z direction for the parameter i
1560 fCalibraMode->SetNz(i, Nz);
1563 AliInfo("You have to choose between 0 and 4");
1567 //_____________________________________________________________________________
1568 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1571 // Set the mode of calibration group in the rphi direction for the parameter i
1576 fCalibraMode->SetNrphi(i ,Nrphi);
1579 AliInfo("You have to choose between 0 and 6");
1584 //_____________________________________________________________________________
1585 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1588 // Set the mode of calibration group all together
1590 if(fVector2d == kTRUE) {
1591 AliInfo("Can not work with the vector method");
1594 fCalibraMode->SetAllTogether(i);
1598 //_____________________________________________________________________________
1599 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1602 // Set the mode of calibration group per supermodule
1604 if(fVector2d == kTRUE) {
1605 AliInfo("Can not work with the vector method");
1608 fCalibraMode->SetPerSuperModule(i);
1612 //____________Set the pad calibration variables for the detector_______________
1613 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1616 // For the detector calcul the first Xbins and set the number of row
1617 // and col pads per calibration groups, the number of calibration
1618 // groups in the detector.
1621 // first Xbins of the detector
1623 fCalibraMode->CalculXBins(detector,0);
1626 fCalibraMode->CalculXBins(detector,1);
1629 fCalibraMode->CalculXBins(detector,2);
1632 // fragmentation of idect
1633 for (Int_t i = 0; i < 3; i++) {
1634 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1635 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1636 , (Int_t) GetStack(detector)
1637 , (Int_t) GetSector(detector),i);
1643 //_____________________________________________________________________________
1644 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1647 // Should be between 0 and 6
1650 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1651 AliInfo("The number of groups must be between 0 and 6!");
1654 fNgroupprf = numberGroupsPRF;
1658 ///////////////////////////////////////////////////////////////////////////////////////////
1659 // Per tracklet: store or reset the info, fill the histos with the info
1660 //////////////////////////////////////////////////////////////////////////////////////////
1661 //_____________________________________________________________________________
1662 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)
1665 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1666 // Correct from the gain correction before
1667 // cls is shared cluster if any
1668 // Return the charge
1672 //printf("StoreInfoCHPHtrack\n");
1674 // time bin of the cluster not corrected
1675 Int_t time = cl->GetPadTime();
1676 Float_t charge = TMath::Abs(cl->GetQ());
1678 charge += TMath::Abs(cls->GetQ());
1679 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1682 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1684 //Correct for the gain coefficient used in the database for reconstruction
1685 Float_t correctthegain = 1.0;
1686 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1687 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1688 Float_t correction = 1.0;
1689 Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1690 // we divide with gain in AliTRDclusterizer::Transform...
1691 if( correctthegain > 0 ) normalisation /= correctthegain;
1694 // take dd/dl corrected from the angle of the track
1695 correction = dqdl / (normalisation);
1698 // Fill the fAmpTotal with the charge
1700 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1701 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1702 fAmpTotal[(Int_t) group[0]] += correction;
1706 // Fill the fPHPlace and value
1708 if (time>=0 && time<fTimeMax) {
1709 fPHPlace[time] = group[1];
1710 fPHValue[time] = charge;
1717 //____________Offine tracking in the AliTRDtracker_____________________________
1718 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1721 // Reset values per tracklet
1724 //Reset good tracklet
1725 fGoodTracklet = kTRUE;
1727 // Reset the fPHValue
1729 //Reset the fPHValue and fPHPlace
1730 for (Int_t k = 0; k < fTimeMax; k++) {
1736 // Reset the fAmpTotal where we put value
1738 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1743 //____________Offine tracking in the AliTRDtracker_____________________________
1744 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1747 // For the offline tracking or mcm tracklets
1748 // This function will be called in the functions UpdateHistogram...
1749 // to fill the info of a track for the relativ gain calibration
1752 Int_t nb = 0; // Nombre de zones traversees
1753 Int_t fd = -1; // Premiere zone non nulle
1754 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1756 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1758 if(nbclusters < fNumberClusters) return;
1759 if(nbclusters > fNumberClustersf) return;
1762 // Normalize with the number of clusters
1763 Double_t normalizeCst = fRelativeScale;
1764 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1766 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1768 // See if the track goes through different zones
1769 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1770 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1771 if (fAmpTotal[k] > 0.0) {
1772 totalcharge += fAmpTotal[k];
1780 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1786 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1788 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1789 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1792 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1796 if ((fAmpTotal[fd] > 0.0) &&
1797 (fAmpTotal[fd+1] > 0.0)) {
1798 // One of the two very big
1799 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1801 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1802 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1805 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1808 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1810 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1812 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1813 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1816 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1819 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1822 if (fCalibraMode->GetNfragZ(0) > 1) {
1823 if (fAmpTotal[fd] > 0.0) {
1824 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1825 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1826 // One of the two very big
1827 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1829 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1830 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1833 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1836 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1838 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1840 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1841 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1844 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1846 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1857 //____________Offine tracking in the AliTRDtracker_____________________________
1858 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1861 // For the offline tracking or mcm tracklets
1862 // This function will be called in the functions UpdateHistogram...
1863 // to fill the info of a track for the drift velocity calibration
1866 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1867 Int_t fd1 = -1; // Premiere zone non nulle
1868 Int_t fd2 = -1; // Deuxieme zone non nulle
1869 Int_t k1 = -1; // Debut de la premiere zone
1870 Int_t k2 = -1; // Debut de la seconde zone
1871 Int_t nbclusters = 0; // number of clusters
1875 // See if the track goes through different zones
1876 for (Int_t k = 0; k < fTimeMax; k++) {
1877 if (fPHValue[k] > 0.0) {
1883 if (fPHPlace[k] != fd1) {
1889 if (fPHPlace[k] != fd2) {
1896 // See if noise before and after
1897 if(fMaxCluster > 0) {
1898 if(fPHValue[0] > fMaxCluster) return;
1899 if(fTimeMax > fNbMaxCluster) {
1900 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1901 if(fPHValue[k] > fMaxCluster) return;
1906 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1908 if(nbclusters < fNumberClusters) return;
1909 if(nbclusters > fNumberClustersf) return;
1915 for (Int_t i = 0; i < fTimeMax; i++) {
1917 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1919 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1921 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1924 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1926 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1932 if ((fd1 == fd2+1) ||
1934 // One of the two fast all the think
1935 if (k2 > (k1+fDifference)) {
1936 //we choose to fill the fd1 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)+fd1)+0.5,(Float_t) fPHValue[i]);
1942 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1946 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1948 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1953 if ((k2+fDifference) < fTimeMax) {
1954 //we choose to fill the fd2 with all the values
1956 for (Int_t i = 0; i < fTimeMax; i++) {
1958 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1960 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1964 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1966 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1972 // Two zones voisines sinon rien!
1973 if (fCalibraMode->GetNfragZ(1) > 1) {
1975 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1976 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1977 // One of the two fast all the think
1978 if (k2 > (k1+fDifference)) {
1979 //we choose to fill the fd1 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)+fd1)+0.5,(Float_t) fPHValue[i]);
1985 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1989 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1991 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1996 if ((k2+fDifference) < fTimeMax) {
1997 //we choose to fill the fd2 with all the values
1999 for (Int_t i = 0; i < fTimeMax; i++) {
2001 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2003 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2007 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2009 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2016 // Two zones voisines sinon rien!
2018 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2019 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2020 // One of the two fast all the think
2021 if (k2 > (k1 + fDifference)) {
2022 //we choose to fill the fd1 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)+fd1)+0.5,(Float_t) fPHValue[i]);
2028 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2032 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2034 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2039 if ((k2+fDifference) < fTimeMax) {
2040 //we choose to fill the fd2 with all the values
2042 for (Int_t i = 0; i < fTimeMax; i++) {
2044 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2046 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2050 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2052 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2064 //////////////////////////////////////////////////////////////////////////////////////////
2065 // DAQ process functions
2066 /////////////////////////////////////////////////////////////////////////////////////////
2067 //_____________________________________________________________________
2068 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2071 // Event Processing loop - AliTRDrawStream
2073 // 0 timebin problem
2076 // Same algorithm as TestBeam but different reader
2079 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2081 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2082 digitsManager->CreateArrays();
2084 Int_t withInput = 1;
2086 Double_t phvalue[16][144][36];
2087 for(Int_t k = 0; k < 36; k++){
2088 for(Int_t j = 0; j < 16; j++){
2089 for(Int_t c = 0; c < 144; c++){
2090 phvalue[j][c][k] = 0.0;
2095 fDetectorPreviousTrack = -1;
2099 Int_t nbtimebin = 0;
2100 Int_t baseline = 10;
2106 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2108 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2109 // printf("there is ADC data on this chamber!\n");
2111 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2112 if (digits->HasData()) { //array
2114 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2115 if (indexes->IsAllocated() == kFALSE) {
2116 AliError("Indexes do not exist!");
2121 indexes->ResetCounters();
2123 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2124 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2125 //while (rawStream->Next()) {
2127 Int_t idetector = det; // current detector
2128 //Int_t imcm = rawStream->GetMCM(); // current MCM
2129 //Int_t irob = rawStream->GetROB(); // current ROB
2132 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2134 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2137 for(Int_t k = 0; k < 36; k++){
2138 for(Int_t j = 0; j < 16; j++){
2139 for(Int_t c = 0; c < 144; c++){
2140 phvalue[j][c][k] = 0.0;
2146 fDetectorPreviousTrack = idetector;
2147 //fMCMPrevious = imcm;
2148 //fROBPrevious = irob;
2150 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2151 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2152 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2153 baseline = digitParam->GetADCbaseline(det); // baseline
2155 if(nbtimebin == 0) return 0;
2156 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2157 fTimeMax = nbtimebin;
2159 fNumberClustersf = fTimeMax;
2160 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2163 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2164 // phvalue[row][col][itime] = signal[itime]-baseline;
2165 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2166 /*if(phvalue[iRow][iCol][itime] >= 20) {
2167 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2171 (Short_t)(digits->GetData(iRow,iCol,itime)),
2178 // fill the last one
2179 if(fDetectorPreviousTrack != -1){
2182 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2183 // printf("\n ---> withinput %d\n\n",withInput);
2185 for(Int_t k = 0; k < 36; k++){
2186 for(Int_t j = 0; j < 16; j++){
2187 for(Int_t c = 0; c < 144; c++){
2188 phvalue[j][c][k] = 0.0;
2196 digitsManager->ClearArrays(det);
2198 delete digitsManager;
2203 //_____________________________________________________________________
2204 //////////////////////////////////////////////////////////////////////////////
2205 // Routine inside the DAQ process
2206 /////////////////////////////////////////////////////////////////////////////
2207 //_______________________________________________________________________
2208 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2211 // Look for the maximum by collapsing over the time
2212 // Sum over four pad col and two pad row
2218 Int_t idect = fDetectorPreviousTrack;
2219 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2221 for(Int_t tb = 0; tb < 36; tb++){
2225 //fGoodTracklet = kTRUE;
2226 //fDetectorPreviousTrack = 0;
2229 ///////////////////////////
2231 /////////////////////////
2235 Double_t integralMax = -1;
2237 for (Int_t ir = 1; ir <= 15; ir++)
2239 for (Int_t ic = 2; ic <= 142; ic++)
2241 Double_t integral = 0;
2242 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2244 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2246 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2247 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2250 for(Int_t tb = 0; tb< fTimeMax; tb++){
2251 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2256 if (integralMax < integral)
2260 integralMax = integral;
2262 } // check max integral
2266 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2267 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2272 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2276 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2277 //if(!fGoodTracklet) used = 1;;
2279 // /////////////////////////////////////////////////////
2280 // sum ober 2 row and 4 pad cols for each time bins
2281 // ////////////////////////////////////////////////////
2285 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2287 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2289 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2290 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2292 for(Int_t it = 0; it < fTimeMax; it++){
2293 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2300 Double_t sumcharge = 0.0;
2301 for(Int_t it = 0; it < fTimeMax; it++){
2302 sumcharge += sum[it];
2303 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2307 /////////////////////////////////////////////////////////
2309 ////////////////////////////////////////////////////////
2310 if(fDebugLevel > 1){
2311 if ( !fDebugStreamer ) {
2313 TDirectory *backup = gDirectory;
2314 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2315 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2318 Double_t amph0 = sum[0];
2319 Double_t amphlast = sum[fTimeMax-1];
2320 Double_t rms = TMath::RMS(fTimeMax,sum);
2321 Int_t goodtracklet = (Int_t) fGoodTracklet;
2322 for(Int_t it = 0; it < fTimeMax; it++){
2323 Double_t clustera = sum[it];
2325 (* fDebugStreamer) << "FillDAQa"<<
2326 "ampTotal="<<sumcharge<<
2329 "detector="<<idect<<
2331 "amphlast="<<amphlast<<
2332 "goodtracklet="<<goodtracklet<<
2333 "clustera="<<clustera<<
2341 ////////////////////////////////////////////////////////
2343 ///////////////////////////////////////////////////////
2344 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2345 if(sum[0] > 100.0) used = 1;
2346 if(nbcl < fNumberClusters) used = 1;
2347 if(nbcl > fNumberClustersf) used = 1;
2349 //if(fDetectorPreviousTrack == 15){
2350 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2352 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2354 for(Int_t it = 0; it < fTimeMax; it++){
2355 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2357 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2359 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2361 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2366 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2368 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2375 //____________Online trackling in AliTRDtrigger________________________________
2376 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2380 // Fill a simple average pulse height
2384 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2390 //____________Write_____________________________________________________
2391 //_____________________________________________________________________
2392 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2395 // Write infos to file
2399 if ( fDebugStreamer ) {
2400 delete fDebugStreamer;
2401 fDebugStreamer = 0x0;
2404 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2409 ,fNumberUsedPh[1]));
2411 TDirectory *backup = gDirectory;
2417 option = "recreate";
2419 TFile f(filename,option.Data());
2421 TStopwatch stopwatch;
2424 f.WriteTObject(fCalibraVector);
2429 f.WriteTObject(fCH2d);
2434 f.WriteTObject(fPH2d);
2439 f.WriteTObject(fPRF2d);
2442 if(fLinearFitterOn){
2443 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2444 f.WriteTObject(fLinearVdriftFit);
2449 if ( backup ) backup->cd();
2451 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2452 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2454 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2456 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2457 //___________________________________________probe the histos__________________________________________________
2458 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2461 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2462 // debug mode with 2 for TH2I and 3 for TProfile2D
2463 // It gives a pointer to a Double_t[7] with the info following...
2464 // [0] : number of calibration groups with entries
2465 // [1] : minimal number of entries found
2466 // [2] : calibration group number of the min
2467 // [3] : maximal number of entries found
2468 // [4] : calibration group number of the max
2469 // [5] : mean number of entries found
2470 // [6] : mean relative error
2473 Double_t *info = new Double_t[7];
2475 // Number of Xbins (detectors or groups of pads)
2476 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2477 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2480 Double_t nbwe = 0; //number of calibration groups with entries
2481 Double_t minentries = 0; //minimal number of entries found
2482 Double_t maxentries = 0; //maximal number of entries found
2483 Double_t placemin = 0; //calibration group number of the min
2484 Double_t placemax = -1; //calibration group number of the max
2485 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2486 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2488 Double_t counter = 0;
2491 TH1F *nbEntries = 0x0;//distribution of the number of entries
2492 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2493 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2495 // Beginning of the loop over the calibration groups
2496 for (Int_t idect = 0; idect < nbins; idect++) {
2498 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2499 projch->SetDirectory(0);
2501 // Number of entries for this calibration group
2502 Double_t nentries = 0.0;
2504 for (Int_t k = 0; k < nxbins; k++) {
2505 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2509 for (Int_t k = 0; k < nxbins; k++) {
2510 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2511 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2512 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2521 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2522 nbEntries->SetDirectory(0);
2523 nbEntries->Fill(nentries);
2524 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2525 nbEntriesPerGroup->SetDirectory(0);
2526 nbEntriesPerGroup->Fill(idect+0.5,nentries);
2527 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));
2528 nbEntriesPerSp->SetDirectory(0);
2529 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2534 if(nentries > maxentries){
2535 maxentries = nentries;
2539 minentries = nentries;
2541 if(nentries < minentries){
2542 minentries = nentries;
2548 meanstats += nentries;
2550 }//calibration groups loop
2552 if(nbwe > 0) meanstats /= nbwe;
2553 if(counter > 0) meanrelativerror /= counter;
2555 AliInfo(Form("There are %f calibration groups with entries",nbwe));
2556 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2557 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2558 AliInfo(Form("The mean number of entries is %f",meanstats));
2559 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2562 info[1] = minentries;
2564 info[3] = maxentries;
2566 info[5] = meanstats;
2567 info[6] = meanrelativerror;
2569 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2570 gStyle->SetPalette(1);
2571 gStyle->SetOptStat(1111);
2572 gStyle->SetPadBorderMode(0);
2573 gStyle->SetCanvasColor(10);
2574 gStyle->SetPadLeftMargin(0.13);
2575 gStyle->SetPadRightMargin(0.01);
2576 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2579 nbEntries->Draw("");
2581 nbEntriesPerSp->SetStats(0);
2582 nbEntriesPerSp->Draw("");
2583 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2585 nbEntriesPerGroup->SetStats(0);
2586 nbEntriesPerGroup->Draw("");
2592 //____________________________________________________________________________
2593 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2596 // Return a Int_t[4] with:
2597 // 0 Mean number of entries
2598 // 1 median of number of entries
2599 // 2 rms of number of entries
2600 // 3 number of group with entries
2603 Double_t *stat = new Double_t[4];
2606 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
2608 Double_t *weight = new Double_t[nbofgroups];
2609 Double_t *nonul = new Double_t[nbofgroups];
2611 for(Int_t k = 0; k < nbofgroups; k++){
2612 if(fEntriesCH[k] > 0) {
2614 nonul[(Int_t)stat[3]] = fEntriesCH[k];
2617 else weight[k] = 0.0;
2619 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
2620 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
2621 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2629 //____________________________________________________________________________
2630 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2633 // Return a Int_t[4] with:
2634 // 0 Mean number of entries
2635 // 1 median of number of entries
2636 // 2 rms of number of entries
2637 // 3 number of group with entries
2640 Double_t *stat = new Double_t[4];
2643 Int_t nbofgroups = 540;
2644 Double_t *weight = new Double_t[nbofgroups];
2645 Int_t *nonul = new Int_t[nbofgroups];
2647 for(Int_t k = 0; k < nbofgroups; k++){
2648 if(fEntriesLinearFitter[k] > 0) {
2650 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2653 else weight[k] = 0.0;
2655 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2656 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2657 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
2665 //////////////////////////////////////////////////////////////////////////////////////
2667 //////////////////////////////////////////////////////////////////////////////////////
2668 //_____________________________________________________________________________
2669 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2672 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2673 // If fNgroupprf is zero then no binning in tnp
2677 name += fCalibraMode->GetNz(2);
2679 name += fCalibraMode->GetNrphi(2);
2683 if(fNgroupprf != 0){
2685 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2686 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2687 fPRF2d->SetYTitle("Det/pad groups");
2688 fPRF2d->SetXTitle("Position x/W [pad width units]");
2689 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2690 fPRF2d->SetStats(0);
2693 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2694 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2695 fPRF2d->SetYTitle("Det/pad groups");
2696 fPRF2d->SetXTitle("Position x/W [pad width units]");
2697 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2698 fPRF2d->SetStats(0);
2703 //_____________________________________________________________________________
2704 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2707 // Create the 2D histos
2710 TString name("Ver");
2711 name += fVersionVdriftUsed;
2713 name += fSubVersionVdriftUsed;
2715 name += fFirstRunVdrift;
2717 name += fCalibraMode->GetNz(1);
2719 name += fCalibraMode->GetNrphi(1);
2721 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2722 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2724 fPH2d->SetYTitle("Det/pad groups");
2725 fPH2d->SetXTitle("time [#mus]");
2726 fPH2d->SetZTitle("<PH> [a.u.]");
2730 //_____________________________________________________________________________
2731 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2734 // Create the 2D histos
2737 TString name("Ver");
2738 name += fVersionGainUsed;
2740 name += fSubVersionGainUsed;
2742 name += fFirstRunGain;
2744 name += fCalibraMode->GetNz(0);
2746 name += fCalibraMode->GetNrphi(0);
2748 fCH2d = new TH2I("CH2d",(const Char_t *) name
2749 ,(Int_t)fNumberBinCharge,0,fRangeHistoCharge,nn,0,nn);
2750 fCH2d->SetYTitle("Det/pad groups");
2751 fCH2d->SetXTitle("charge deposit [a.u]");
2752 fCH2d->SetZTitle("counts");
2757 //////////////////////////////////////////////////////////////////////////////////
2758 // Set relative scale
2759 /////////////////////////////////////////////////////////////////////////////////
2760 //_____________________________________________________________________________
2761 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2764 // Set the factor that will divide the deposited charge
2765 // to fit in the histo range [0,fRangeHistoCharge]
2768 if (RelativeScale > 0.0) {
2769 fRelativeScale = RelativeScale;
2772 AliInfo("RelativeScale must be strict positif!");
2776 //////////////////////////////////////////////////////////////////////////////////
2777 // Quick way to fill a histo
2778 //////////////////////////////////////////////////////////////////////////////////
2779 //_____________________________________________________________________
2780 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2783 // FillCH2d: Marian style
2786 //skip simply the value out of range
2787 if((y>=fRangeHistoCharge) || (y<0.0)) return;
2788 if(fRangeHistoCharge < 0.0) return;
2790 //Calcul the y place
2791 Int_t yplace = (Int_t) (fNumberBinCharge*y/fRangeHistoCharge)+1;
2792 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2795 fCH2d->GetArray()[place]++;
2799 //////////////////////////////////////////////////////////////////////////////////
2800 // Geometrical functions
2801 ///////////////////////////////////////////////////////////////////////////////////
2802 //_____________________________________________________________________________
2803 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2806 // Reconstruct the layer number from the detector number
2809 return ((Int_t) (d % 6));
2813 //_____________________________________________________________________________
2814 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2817 // Reconstruct the stack number from the detector number
2819 const Int_t kNlayer = 6;
2821 return ((Int_t) (d % 30) / kNlayer);
2825 //_____________________________________________________________________________
2826 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2829 // Reconstruct the sector number from the detector number
2833 return ((Int_t) (d / fg));
2836 ///////////////////////////////////////////////////////////////////////////////////
2837 // Getter functions for DAQ of the CH2d and the PH2d
2838 //////////////////////////////////////////////////////////////////////////////////
2839 //_____________________________________________________________________
2840 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2843 // return pointer to fPH2d TProfile2D
2844 // create a new TProfile2D if it doesn't exist allready
2850 fTimeMax = nbtimebin;
2851 fSf = samplefrequency;
2857 //_____________________________________________________________________
2858 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2861 // return pointer to fCH2d TH2I
2862 // create a new TH2I if it doesn't exist allready
2871 ////////////////////////////////////////////////////////////////////////////////////////////
2872 // Drift velocity calibration
2873 ///////////////////////////////////////////////////////////////////////////////////////////
2874 //_____________________________________________________________________
2875 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2878 // return pointer to TLinearFitter Calibration
2879 // if force is true create a new TLinearFitter if it doesn't exist allready
2882 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2883 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2886 // if we are forced and TLinearFitter doesn't yet exist create it
2888 // new TLinearFitter
2889 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2890 fLinearFitterArray.AddAt(linearfitter,detector);
2891 return linearfitter;
2894 //____________________________________________________________________________
2895 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2898 // Analyse array of linear fitter because can not be written
2899 // Store two arrays: one with the param the other one with the error param + number of entries
2902 for(Int_t k = 0; k < 540; k++){
2903 TLinearFitter *linearfitter = GetLinearFitter(k);
2904 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2905 TVectorD *par = new TVectorD(2);
2906 TVectorD pare = TVectorD(2);
2907 TVectorD *parE = new TVectorD(3);
2908 if((linearfitter->EvalRobust(0.8)==0)) {
2909 //linearfitter->Eval();
2910 linearfitter->GetParameters(*par);
2911 //linearfitter->GetErrors(pare);
2912 //Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2913 //(*parE)[0] = pare[0]*ppointError;
2914 //(*parE)[1] = pare[1]*ppointError;
2918 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2919 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2920 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);