1 //**************************************************************************
2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 // * Author: The ALICE Off-line Project. *
5 // * Contributors are mentioned in the code where appropriate. *
7 // * Permission to use, copy, modify and distribute this software and its *
8 // * documentation strictly for non-commercial purposes is hereby granted *
9 // * without fee, provided that the above copyright notice appears in all *
10 // * copies and that both the copyright notice and this permission notice *
11 //* appear in the supporting documentation. The authors make no claims *
12 //* about the suitability of this software for any purpose. It is *
13 //* provided "as is" without express or implied warranty. *
14 //**************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////
20 // AliTRDCalibraFillHisto
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fills histos or vectors.
24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
25 // The user has to choose with the functions SetNz and SetNrphi the precision of the calibration (see AliTRDCalibraMode).
26 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
27 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
28 // in the function "FollowBackProlongation" (AliTRDtracker)
29 // Per default the functions to fill are off.
32 // R. Bailhache (R.Bailhache@gsi.de, rbailhache@ikf.uni-frankfurt.de)
33 // J. Book (jbook@ikf.uni-frankfurt.de)
35 //////////////////////////////////////////////////////////////////////////////////////
37 #include <TProfile2D.h>
42 #include <TObjArray.h>
47 #include <TStopwatch.h>
49 #include <TDirectory.h>
50 #include <TTreeStream.h>
52 #include <TLinearFitter.h>
56 #include "AliTRDCalibraFillHisto.h"
57 #include "AliTRDCalibraMode.h"
58 #include "AliTRDCalibraVector.h"
59 #include "AliTRDCalibraVdriftLinearFit.h"
60 #include "AliTRDcalibDB.h"
61 #include "AliTRDCommonParam.h"
62 #include "AliTRDpadPlane.h"
63 #include "AliTRDcluster.h"
64 #include "AliTRDtrack.h"
65 #include "AliTRDtrackV1.h"
66 #include "AliTRDrawStreamBase.h"
67 #include "AliRawReader.h"
68 #include "AliRawReaderDate.h"
69 #include "AliTRDgeometry.h"
70 #include "./Cal/AliTRDCalROC.h"
71 #include "./Cal/AliTRDCalPad.h"
72 #include "./Cal/AliTRDCalDet.h"
74 #include "AliTRDrawFastStream.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)
145 ,fThresholdClusterPRF2(15.0)
146 ,fLimitChargeIntegration(kFALSE)
147 ,fFillWithZero(kFALSE)
148 ,fNormalizeNbOfCluster(kFALSE)
152 ,fSubVersionGainUsed(0)
153 ,fVersionGainLocalUsed(0)
154 ,fSubVersionGainLocalUsed(0)
155 ,fVersionVdriftUsed(0)
156 ,fSubVersionVdriftUsed(0)
157 ,fCalibraMode(new AliTRDCalibraMode())
160 ,fDetectorPreviousTrack(-1)
164 ,fNumberClustersf(30)
165 ,fNumberClustersProcent(0.5)
166 ,fThresholdClustersDAQ(120.0)
174 ,fNumberBinCharge(50)
180 ,fGoodTracklet(kTRUE)
181 ,fLinearFitterTracklet(0x0)
183 ,fEntriesLinearFitter(0x0)
188 ,fLinearFitterArray(540)
189 ,fLinearVdriftFit(0x0)
194 // Default constructor
198 // Init some default values
201 fNumberUsedCh[0] = 0;
202 fNumberUsedCh[1] = 0;
203 fNumberUsedPh[0] = 0;
204 fNumberUsedPh[1] = 0;
206 fGeo = new AliTRDgeometry();
207 fCalibDB = AliTRDcalibDB::Instance();
210 //______________________________________________________________________________________
211 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
218 ,fPRF2dOn(c.fPRF2dOn)
219 ,fHisto2d(c.fHisto2d)
220 ,fVector2d(c.fVector2d)
221 ,fLinearFitterOn(c.fLinearFitterOn)
222 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
223 ,fRelativeScale(c.fRelativeScale)
224 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
225 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
226 ,fFillWithZero(c.fFillWithZero)
227 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
228 ,fMaxCluster(c.fMaxCluster)
229 ,fNbMaxCluster(c.fNbMaxCluster)
230 ,fVersionGainUsed(c.fVersionGainUsed)
231 ,fSubVersionGainUsed(c.fSubVersionGainUsed)
232 ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
233 ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
234 ,fVersionVdriftUsed(c.fVersionVdriftUsed)
235 ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
238 ,fDebugLevel(c.fDebugLevel)
239 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
240 ,fMCMPrevious(c.fMCMPrevious)
241 ,fROBPrevious(c.fROBPrevious)
242 ,fNumberClusters(c.fNumberClusters)
243 ,fNumberClustersf(c.fNumberClustersf)
244 ,fNumberClustersProcent(c.fNumberClustersProcent)
245 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
246 ,fNumberRowDAQ(c.fNumberRowDAQ)
247 ,fNumberColDAQ(c.fNumberColDAQ)
248 ,fProcent(c.fProcent)
249 ,fDifference(c.fDifference)
250 ,fNumberTrack(c.fNumberTrack)
251 ,fTimeMax(c.fTimeMax)
253 ,fNumberBinCharge(c.fNumberBinCharge)
254 ,fNumberBinPRF(c.fNumberBinPRF)
255 ,fNgroupprf(c.fNgroupprf)
259 ,fGoodTracklet(c.fGoodTracklet)
260 ,fLinearFitterTracklet(0x0)
262 ,fEntriesLinearFitter(0x0)
267 ,fLinearFitterArray(540)
268 ,fLinearVdriftFit(0x0)
275 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
276 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
278 fPH2d = new TProfile2D(*c.fPH2d);
279 fPH2d->SetDirectory(0);
282 fPRF2d = new TProfile2D(*c.fPRF2d);
283 fPRF2d->SetDirectory(0);
286 fCH2d = new TH2I(*c.fCH2d);
287 fCH2d->SetDirectory(0);
289 if(c.fLinearVdriftFit){
290 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
293 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
294 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
299 fGeo = new AliTRDgeometry();
300 fCalibDB = AliTRDcalibDB::Instance();
302 fNumberUsedCh[0] = 0;
303 fNumberUsedCh[1] = 0;
304 fNumberUsedPh[0] = 0;
305 fNumberUsedPh[1] = 0;
309 //____________________________________________________________________________________
310 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
313 // AliTRDCalibraFillHisto destructor
317 if ( fDebugStreamer ) delete fDebugStreamer;
319 if ( fCalDetGain ) delete fCalDetGain;
320 if ( fCalROCGain ) delete fCalROCGain;
322 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
326 delete [] fEntriesCH;
327 delete [] fEntriesLinearFitter;
330 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
331 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
334 if(fLinearVdriftFit) delete fLinearVdriftFit;
340 //_____________________________________________________________________________
341 void AliTRDCalibraFillHisto::Destroy()
352 //_____________________________________________________________________________
353 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
356 // Delete DebugStreamer
359 if ( fDebugStreamer ) delete fDebugStreamer;
360 fDebugStreamer = 0x0;
363 //_____________________________________________________________________________
364 void AliTRDCalibraFillHisto::ClearHistos()
384 //////////////////////////////////////////////////////////////////////////////////
385 // calibration with AliTRDtrackV1: Init, Update
386 //////////////////////////////////////////////////////////////////////////////////
387 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
388 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
391 // Init the histograms and stuff to be filled
396 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
398 AliInfo("Could not get calibDB");
401 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
403 AliInfo("Could not get CommonParam");
408 if(nboftimebin > 0) fTimeMax = nboftimebin;
409 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
410 if(fTimeMax <= 0) fTimeMax = 30;
411 printf("////////////////////////////////////////////\n");
412 printf("Number of time bins in calibration component %d\n",fTimeMax);
413 printf("////////////////////////////////////////////\n");
414 fSf = parCom->GetSamplingFrequency();
415 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
416 else fRelativeScale = 1.18;
417 fNumberClustersf = fTimeMax;
418 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
420 // Init linear fitter
421 if(!fLinearFitterTracklet) {
422 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
423 fLinearFitterTracklet->StoreData(kTRUE);
426 // Calcul Xbins Chambd0, Chamb2
427 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
428 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
429 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
431 // If vector method On initialised all the stuff
433 fCalibraVector = new AliTRDCalibraVector();
434 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
435 fCalibraVector->SetTimeMax(fTimeMax);
436 if(fNgroupprf != 0) {
437 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
438 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
441 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
442 fCalibraVector->SetPRFRange(1.5);
444 for(Int_t k = 0; k < 3; k++){
445 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
446 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
448 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
449 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
450 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
451 fCalibraVector->SetNbGroupPRF(fNgroupprf);
454 // Create the 2D histos corresponding to the pad groupCalibration mode
457 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
458 ,fCalibraMode->GetNz(0)
459 ,fCalibraMode->GetNrphi(0)));
461 // Create the 2D histo
466 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
467 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
471 fEntriesCH = new Int_t[ntotal0];
472 for(Int_t k = 0; k < ntotal0; k++){
479 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
480 ,fCalibraMode->GetNz(1)
481 ,fCalibraMode->GetNrphi(1)));
483 // Create the 2D histo
488 fPHPlace = new Short_t[fTimeMax];
489 for (Int_t k = 0; k < fTimeMax; k++) {
492 fPHValue = new Float_t[fTimeMax];
493 for (Int_t k = 0; k < fTimeMax; k++) {
497 if (fLinearFitterOn) {
498 if(fLinearFitterDebugOn) {
499 fLinearFitterArray.SetName("ArrayLinearFitters");
500 fEntriesLinearFitter = new Int_t[540];
501 for(Int_t k = 0; k < 540; k++){
502 fEntriesLinearFitter[k] = 0;
505 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
510 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
511 ,fCalibraMode->GetNz(2)
512 ,fCalibraMode->GetNrphi(2)));
513 // Create the 2D histo
515 CreatePRF2d(ntotal2);
522 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
523 Bool_t AliTRDCalibraFillHisto::InitCalDet()
526 // Init the Gain Cal Det
531 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainUsed,fSubVersionGainUsed);
533 AliError("No gain det calibration entry found");
536 AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
538 AliError("No calDet gain found");
544 fCalDetGain->~AliTRDCalDet();
545 new(fCalDetGain) AliTRDCalDet(*(calDet));
546 }else fCalDetGain = new AliTRDCalDet(*(calDet));
551 name += fVersionGainUsed;
553 name += fSubVersionGainUsed;
555 name += fCalibraMode->GetNz(0);
557 name += fCalibraMode->GetNrphi(0);
559 fCH2d->SetTitle(name);
562 TString namee("Ver");
563 namee += fVersionVdriftUsed;
565 namee += fSubVersionVdriftUsed;
567 namee += fCalibraMode->GetNz(1);
569 namee += fCalibraMode->GetNrphi(1);
571 fPH2d->SetTitle(namee);
577 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
578 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
581 // Init the Gain Cal Pad
586 AliCDBEntry *entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",AliCDBManager::Instance()->GetRun(),fVersionGainLocalUsed,fSubVersionGainLocalUsed);
588 AliError("No gain pad calibration entry found");
591 AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
593 AliError("No calPad gain found");
596 AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
598 AliError("No calRoc gain found");
603 fCalROCGain->~AliTRDCalROC();
604 new(fCalROCGain) AliTRDCalROC(*(calRoc));
605 }else fCalROCGain = new AliTRDCalROC(*(calRoc));
614 //____________Offline tracking in the AliTRDtracker____________________________
615 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(const AliTRDtrack *t)
618 // Use AliTRDtrack for the calibration
622 AliTRDcluster *cl = 0x0; // pointeur to cluster
623 Int_t index0 = 0; // index of the first cluster in the current chamber
624 Int_t index1 = 0; // index of the current cluster in the current chamber
626 //////////////////////////////
627 // loop over the clusters
628 ///////////////////////////////
629 while((cl = t->GetCluster(index1))){
631 /////////////////////////////////////////////////////////////////////////
632 // Fill the infos for the previous clusters if not the same detector
633 ////////////////////////////////////////////////////////////////////////
634 Int_t detector = cl->GetDetector();
635 if ((detector != fDetectorPreviousTrack) &&
636 (index0 != index1)) {
640 //If the same track, then look if the previous detector is in
641 //the same plane, if yes: not a good track
642 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
646 // Fill only if the track doesn't touch a masked pad or doesn't
649 // drift velocity unables to cut bad tracklets
650 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
654 FillTheInfoOfTheTrackCH(index1-index0);
659 FillTheInfoOfTheTrackPH();
662 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
665 } // if a good tracklet
668 ResetfVariablestracklet();
671 } // Fill at the end the charge
674 /////////////////////////////////////////////////////////////
675 // Localise and take the calib gain object for the detector
676 ////////////////////////////////////////////////////////////
677 if (detector != fDetectorPreviousTrack) {
679 //Localise the detector
680 LocalisationDetectorXbins(detector);
683 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
685 AliInfo("Could not get calibDB");
690 if(!fIsHLT) InitCalPad(detector);
694 // Reset the detectbjobsor
695 fDetectorPreviousTrack = detector;
697 ///////////////////////////////////////
698 // Store the info of the cluster
699 ///////////////////////////////////////
700 Int_t row = cl->GetPadRow();
701 Int_t col = cl->GetPadCol();
702 CheckGoodTrackletV1(cl);
703 Int_t group[2] = {0,0};
704 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
705 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
706 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
710 } // while on clusters
712 ///////////////////////////
713 // Fill the last plane
714 //////////////////////////
715 if( index0 != index1 ){
721 // drift velocity unables to cut bad tracklets
722 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
726 FillTheInfoOfTheTrackCH(index1-index0);
731 FillTheInfoOfTheTrackPH();
734 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
736 } // if a good tracklet
741 ResetfVariablestracklet();
746 //____________Offline tracking in the AliTRDtracker____________________________
747 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t)
750 // Use AliTRDtrackV1 for the calibration
754 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
755 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
756 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
757 Bool_t newtr = kTRUE; // new track
760 // AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
763 AliInfo("Could not get calibDB");
768 AliInfo("Could not get calibDB");
773 ///////////////////////////
774 // loop over the tracklet
775 ///////////////////////////
776 for(Int_t itr = 0; itr < 6; itr++){
778 if(!(tracklet = t->GetTracklet(itr))) continue;
779 if(!tracklet->IsOK()) continue;
781 ResetfVariablestracklet();
783 //////////////////////////////////////////
784 // localisation of the tracklet and dqdl
785 //////////////////////////////////////////
786 Int_t layer = tracklet->GetPlane();
788 while(!(cl = tracklet->GetClusters(ic++))) continue;
789 Int_t detector = cl->GetDetector();
790 if (detector != fDetectorPreviousTrack) {
791 // if not a new track
793 // don't use the rest of this track if in the same plane
794 if (layer == GetLayer(fDetectorPreviousTrack)) {
795 //printf("bad tracklet, same layer for detector %d\n",detector);
799 //Localise the detector bin
800 LocalisationDetectorXbins(detector);
802 if(!fIsHLT) InitCalPad(detector);
805 fDetectorPreviousTrack = detector;
809 ////////////////////////////
810 // loop over the clusters
811 ////////////////////////////
812 Int_t nbclusters = 0;
813 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
814 if(!(cl = tracklet->GetClusters(jc))) continue;
817 // Store the info bis of the tracklet
818 Int_t row = cl->GetPadRow();
819 Int_t col = cl->GetPadCol();
820 CheckGoodTrackletV1(cl);
821 Int_t group[2] = {0,0};
822 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
823 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
824 // Add the charge if shared cluster
825 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
827 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col,cls);
830 ////////////////////////////////////////
831 // Fill the stuffs if a good tracklet
832 ////////////////////////////////////////
835 // drift velocity unables to cut bad tracklets
836 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
838 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
842 FillTheInfoOfTheTrackCH(nbclusters);
847 FillTheInfoOfTheTrackPH();
850 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
852 } // if a good tracklet
858 ///////////////////////////////////////////////////////////////////////////////////
859 // Routine inside the update with AliTRDtrack
860 ///////////////////////////////////////////////////////////////////////////////////
861 //____________Offine tracking in the AliTRDtracker_____________________________
862 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1)
865 // Drift velocity calibration:
866 // Fit the clusters with a straight line
867 // From the slope find the drift velocity
870 //Number of points: if less than 3 return kFALSE
871 Int_t npoints = index1-index0;
872 if(npoints <= 2) return kFALSE;
877 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
878 // parameters of the track
879 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
880 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
881 Double_t tnp = 0.0; // tan angle in the plan xy track
882 if( TMath::Abs(snp) < 1.){
883 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
885 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
886 // tilting pad and cross row
887 Int_t crossrow = 0; // if it crosses a pad row
888 Int_t rowp = -1; // if it crosses a pad row
889 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
890 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
891 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
893 fLinearFitterTracklet->ClearPoints();
894 Double_t dydt = 0.0; // dydt tracklet after straight line fit
895 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
896 Double_t pointError = 0.0; // error after straight line fit
897 Int_t nbli = 0; // number linear fitter points
899 //////////////////////////////
900 // loop over clusters
901 ////////////////////////////
902 for(Int_t k = 0; k < npoints; k++){
904 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
905 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
906 Double_t ycluster = cl->GetY();
907 Int_t time = cl->GetPadTime();
908 Double_t timeis = time/fSf;
909 //Double_t q = cl->GetQ();
910 //See if cross two pad rows
911 Int_t row = cl->GetPadRow();
913 if(row != rowp) crossrow = 1;
915 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
920 //////////////////////////////
922 //////////////////////////////
924 fLinearFitterTracklet->ClearPoints();
928 fLinearFitterTracklet->Eval();
929 fLinearFitterTracklet->GetParameters(pars);
930 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
931 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
933 fLinearFitterTracklet->ClearPoints();
935 /////////////////////////////
937 ////////////////////////////
939 if ( !fDebugStreamer ) {
941 TDirectory *backup = gDirectory;
942 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
943 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
947 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
948 //"snpright="<<snpright<<
949 "npoints="<<npoints<<
951 "detector="<<detector<<
958 "crossrow="<<crossrow<<
959 "errorpar="<<errorpar<<
960 "pointError="<<pointError<<
964 Int_t nbclusters = index1-index0;
965 Int_t layer = GetLayer(fDetectorPreviousTrack);
967 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
968 //"snpright="<<snpright<<
969 "nbclusters="<<nbclusters<<
970 "detector="<<fDetectorPreviousTrack<<
976 ///////////////////////////
978 ///////////////////////////
979 if(npoints < fNumberClusters) return kFALSE;
980 if(npoints > fNumberClustersf) return kFALSE;
981 if(pointError >= 0.1) return kFALSE;
982 if(crossrow == 1) return kFALSE;
984 ////////////////////////////
986 ////////////////////////////
988 //Add to the linear fitter of the detector
989 if( TMath::Abs(snp) < 1.){
990 Double_t x = tnp-dzdx*tnt;
991 if(fLinearFitterDebugOn) {
992 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
993 fEntriesLinearFitter[detector]++;
995 fLinearVdriftFit->Update(detector,x,pars[1]);
1001 //____________Offine tracking in the AliTRDtracker_____________________________
1002 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1005 // Drift velocity calibration:
1006 // Fit the clusters with a straight line
1007 // From the slope find the drift velocity
1010 ////////////////////////////////////////////////
1011 //Number of points: if less than 3 return kFALSE
1012 /////////////////////////////////////////////////
1013 if(nbclusters <= 2) return kFALSE;
1018 // results of the linear fit
1019 Double_t dydt = 0.0; // dydt tracklet after straight line fit
1020 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
1021 Double_t pointError = 0.0; // error after straight line fit
1022 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
1023 Int_t crossrow = 0; // if it crosses a pad row
1024 Int_t rowp = -1; // if it crosses a pad row
1025 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
1026 fLinearFitterTracklet->ClearPoints();
1029 ///////////////////////////////////////////
1030 // Take the parameters of the track
1031 //////////////////////////////////////////
1032 // take now the snp, tnp and tgl from the track
1033 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1034 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1035 if( TMath::Abs(snp) < 1.){
1036 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1038 Double_t tgl = tracklet->GetTgl(); // dz/dl
1039 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1041 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1042 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1043 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1044 // at the end with correction due to linear fit
1045 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1046 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1049 ////////////////////////////
1050 // loop over the clusters
1051 ////////////////////////////
1053 AliTRDcluster *cl = 0x0;
1054 //////////////////////////////
1055 // Check no shared clusters
1056 //////////////////////////////
1057 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
1058 cl = tracklet->GetClusters(icc);
1059 if(cl) crossrow = 1;
1061 //////////////////////////////////
1063 //////////////////////////////////
1064 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1065 if(!(cl = tracklet->GetClusters(ic))) continue;
1066 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
1068 Double_t ycluster = cl->GetY();
1069 Int_t time = cl->GetPadTime();
1070 Double_t timeis = time/fSf;
1071 //See if cross two pad rows
1072 Int_t row = cl->GetPadRow();
1073 if(rowp==-1) rowp = row;
1074 if(row != rowp) crossrow = 1;
1076 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
1082 ////////////////////////////////////
1083 // Do the straight line fit now
1084 ///////////////////////////////////
1086 fLinearFitterTracklet->ClearPoints();
1090 fLinearFitterTracklet->Eval();
1091 fLinearFitterTracklet->GetParameters(pars);
1092 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
1093 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
1095 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
1096 fLinearFitterTracklet->ClearPoints();
1098 ////////////////////////////////
1100 ///////////////////////////////
1103 if(fDebugLevel > 0){
1104 if ( !fDebugStreamer ) {
1106 TDirectory *backup = gDirectory;
1107 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1108 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1112 Int_t layer = GetLayer(fDetectorPreviousTrack);
1114 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1115 //"snpright="<<snpright<<
1117 "nbclusters="<<nbclusters<<
1118 "detector="<<fDetectorPreviousTrack<<
1126 "crossrow="<<crossrow<<
1127 "errorpar="<<errorpar<<
1128 "pointError="<<pointError<<
1133 /////////////////////////
1135 ////////////////////////
1137 if(nbclusters < fNumberClusters) return kFALSE;
1138 if(nbclusters > fNumberClustersf) return kFALSE;
1139 if(pointError >= 0.3) return kFALSE;
1140 if(crossrow == 1) return kFALSE;
1142 ///////////////////////
1144 //////////////////////
1146 if(fLinearFitterOn){
1147 //Add to the linear fitter of the detector
1148 if( TMath::Abs(snp) < 1.){
1149 Double_t x = tnp-dzdx*tnt;
1150 if(fLinearFitterDebugOn) {
1151 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1152 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1154 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1160 //____________Offine tracking in the AliTRDtracker_____________________________
1161 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1)
1164 // PRF width calibration
1165 // Assume a Gaussian shape: determinate the position of the three pad clusters
1166 // Fit with a straight line
1167 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1168 // Fill the PRF as function of angle of the track
1173 //////////////////////////
1175 /////////////////////////
1176 Int_t npoints = index1-index0; // number of total points
1177 Int_t nb3pc = 0; // number of three pads clusters used for fit
1178 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1179 // To see the difference due to the fit
1180 Double_t *padPositions;
1181 padPositions = new Double_t[npoints];
1182 for(Int_t k = 0; k < npoints; k++){
1183 padPositions[k] = 0.0;
1185 // Take the tgl and snp with the track t now
1186 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1187 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1188 Float_t dzdx = 0.0; // dzdx
1190 if(TMath::Abs(snp) < 1.0){
1191 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1192 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1195 fLinearFitterTracklet->ClearPoints();
1197 ///////////////////////////
1198 // calcul the tnp group
1199 ///////////////////////////
1200 Bool_t echec = kFALSE;
1201 Double_t shift = 0.0;
1202 //Calculate the shift in x coresponding to this tnp
1203 if(fNgroupprf != 0.0){
1204 shift = -3.0*(fNgroupprf-1)-1.5;
1205 Double_t limithigh = -0.2*(fNgroupprf-1);
1206 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1208 while(tnp > limithigh){
1215 delete [] padPositions;
1219 //////////////////////
1221 /////////////////////
1222 for(Int_t k = 0; k < npoints; k++){
1224 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1226 Short_t *signals = cl->GetSignals();
1227 Double_t time = cl->GetPadTime();
1228 //Calculate x if possible
1229 Float_t xcenter = 0.0;
1230 Bool_t echec1 = kTRUE;
1231 if((time<=7) || (time>=21)) continue;
1232 // Center 3 balanced: position with the center of the pad
1233 if ((((Float_t) signals[3]) > 0.0) &&
1234 (((Float_t) signals[2]) > 0.0) &&
1235 (((Float_t) signals[4]) > 0.0)) {
1237 // Security if the denomiateur is 0
1238 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1239 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1240 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1241 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1242 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1248 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1250 //if no echec: calculate with the position of the pad
1251 // Position of the cluster
1252 Double_t padPosition = xcenter + cl->GetPadCol();
1253 padPositions[k] = padPosition;
1255 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1259 /////////////////////////////
1261 ////////////////////////////
1263 delete [] padPositions;
1264 fLinearFitterTracklet->ClearPoints();
1267 fLinearFitterTracklet->Eval();
1269 fLinearFitterTracklet->GetParameters(line);
1270 Float_t pointError = -1.0;
1271 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1272 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1274 fLinearFitterTracklet->ClearPoints();
1276 /////////////////////////////////////////////////////
1277 // Now fill the PRF: second loop over clusters
1278 /////////////////////////////////////////////////////
1279 for(Int_t k = 0; k < npoints; k++){
1281 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1282 Short_t *signals = cl->GetSignals(); // signal
1283 Double_t time = cl->GetPadTime(); // time bin
1284 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1285 Float_t padPos = cl->GetPadCol(); // middle pad
1286 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1287 Float_t ycenter = 0.0; // relative center charge
1288 Float_t ymin = 0.0; // relative left charge
1289 Float_t ymax = 0.0; // relative right charge
1291 //Requiere simply two pads clusters at least
1292 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1293 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1294 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1295 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1296 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1297 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1301 Int_t rowcl = cl->GetPadRow(); // row of cluster
1302 Int_t colcl = cl->GetPadCol(); // col of cluster
1303 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1304 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1305 Float_t xcl = cl->GetY(); // y cluster
1306 Float_t qcl = cl->GetQ(); // charge cluster
1307 Int_t layer = GetLayer(detector); // layer
1308 Int_t stack = GetStack(detector); // stack
1309 Double_t xdiff = dpad; // reconstructed position constant
1310 Double_t x = dpad; // reconstructed position moved
1311 Float_t ep = pointError; // error of fit
1312 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1313 Float_t signal3 = (Float_t)signals[3]; // signal
1314 Float_t signal2 = (Float_t)signals[2]; // signal
1315 Float_t signal4 = (Float_t)signals[4]; // signal
1316 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1318 //////////////////////////////
1320 /////////////////////////////
1321 if(fDebugLevel > 0){
1322 if ( !fDebugStreamer ) {
1324 TDirectory *backup = gDirectory;
1325 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1326 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1332 Float_t y = ycenter;
1333 (* fDebugStreamer) << "HandlePRFtracklet"<<
1334 "caligroup="<<caligroup<<
1335 "detector="<<detector<<
1338 "npoints="<<npoints<<
1347 "padPosition="<<padPositions[k]<<
1348 "padPosTracklet="<<padPosTracklet<<
1353 "signal1="<<signal1<<
1354 "signal2="<<signal2<<
1355 "signal3="<<signal3<<
1356 "signal4="<<signal4<<
1357 "signal5="<<signal5<<
1363 (* fDebugStreamer) << "HandlePRFtracklet"<<
1364 "caligroup="<<caligroup<<
1365 "detector="<<detector<<
1368 "npoints="<<npoints<<
1377 "padPosition="<<padPositions[k]<<
1378 "padPosTracklet="<<padPosTracklet<<
1383 "signal1="<<signal1<<
1384 "signal2="<<signal2<<
1385 "signal3="<<signal3<<
1386 "signal4="<<signal4<<
1387 "signal5="<<signal5<<
1393 (* fDebugStreamer) << "HandlePRFtracklet"<<
1394 "caligroup="<<caligroup<<
1395 "detector="<<detector<<
1398 "npoints="<<npoints<<
1407 "padPosition="<<padPositions[k]<<
1408 "padPosTracklet="<<padPosTracklet<<
1413 "signal1="<<signal1<<
1414 "signal2="<<signal2<<
1415 "signal3="<<signal3<<
1416 "signal4="<<signal4<<
1417 "signal5="<<signal5<<
1423 ////////////////////////////
1425 ///////////////////////////
1426 if(npoints < fNumberClusters) continue;
1427 if(npoints > fNumberClustersf) continue;
1428 if(nb3pc <= 5) continue;
1429 if((time >= 21) || (time < 7)) continue;
1430 if(TMath::Abs(snp) >= 1.0) continue;
1431 if(TMath::Abs(qcl) < 80) continue;
1433 ////////////////////////////
1435 ///////////////////////////
1437 if(TMath::Abs(dpad) < 1.5) {
1438 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1439 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1441 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1442 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1443 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1445 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1446 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1447 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1451 if(TMath::Abs(dpad) < 1.5) {
1452 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1453 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1455 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1456 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1457 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1459 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1460 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1461 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1465 delete [] padPositions;
1469 //____________Offine tracking in the AliTRDtracker_____________________________
1470 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1473 // PRF width calibration
1474 // Assume a Gaussian shape: determinate the position of the three pad clusters
1475 // Fit with a straight line
1476 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1477 // Fill the PRF as function of angle of the track
1481 //printf("begin\n");
1482 ///////////////////////////////////////////
1483 // Take the parameters of the track
1484 //////////////////////////////////////////
1485 // take now the snp, tnp and tgl from the track
1486 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1487 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1488 if( TMath::Abs(snp) < 1.){
1489 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1491 Double_t tgl = tracklet->GetTgl(); // dz/dl
1492 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1494 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1495 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1496 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1497 // at the end with correction due to linear fit
1498 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1499 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1501 ///////////////////////////////
1502 // Calculate tnp group shift
1503 ///////////////////////////////
1504 Bool_t echec = kFALSE;
1505 Double_t shift = 0.0;
1506 //Calculate the shift in x coresponding to this tnp
1507 if(fNgroupprf != 0.0){
1508 shift = -3.0*(fNgroupprf-1)-1.5;
1509 Double_t limithigh = -0.2*(fNgroupprf-1);
1510 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1512 while(tnp > limithigh){
1518 // do nothing if out of tnp range
1519 //printf("echec %d\n",(Int_t)echec);
1520 if(echec) return kFALSE;
1522 ///////////////////////
1524 //////////////////////
1526 Int_t nb3pc = 0; // number of three pads clusters used for fit
1527 // to see the difference between the fit and the 3 pad clusters position
1528 Double_t padPositions[AliTRDseedV1::kNtb];
1529 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1530 fLinearFitterTracklet->ClearPoints();
1532 //printf("loop clusters \n");
1533 ////////////////////////////
1534 // loop over the clusters
1535 ////////////////////////////
1536 AliTRDcluster *cl = 0x0;
1537 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1538 // reject shared clusters on pad row
1539 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1540 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1543 cl = tracklet->GetClusters(ic);
1545 Double_t time = cl->GetPadTime();
1546 if((time<=7) || (time>=21)) continue;
1547 Short_t *signals = cl->GetSignals();
1548 Float_t xcenter = 0.0;
1549 Bool_t echec1 = kTRUE;
1551 /////////////////////////////////////////////////////////////
1552 // Center 3 balanced: position with the center of the pad
1553 /////////////////////////////////////////////////////////////
1554 if ((((Float_t) signals[3]) > 0.0) &&
1555 (((Float_t) signals[2]) > 0.0) &&
1556 (((Float_t) signals[4]) > 0.0)) {
1558 // Security if the denomiateur is 0
1559 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1560 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1561 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1562 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1563 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1569 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1570 if(echec1) continue;
1572 ////////////////////////////////////////////////////////
1573 //if no echec1: calculate with the position of the pad
1574 // Position of the cluster
1575 // fill the linear fitter
1576 ///////////////////////////////////////////////////////
1577 Double_t padPosition = xcenter + cl->GetPadCol();
1578 padPositions[ic] = padPosition;
1580 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1585 //printf("Fin loop clusters \n");
1586 //////////////////////////////
1587 // fit with a straight line
1588 /////////////////////////////
1590 fLinearFitterTracklet->ClearPoints();
1593 fLinearFitterTracklet->Eval();
1595 fLinearFitterTracklet->GetParameters(line);
1596 Float_t pointError = -1.0;
1597 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1598 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1600 fLinearFitterTracklet->ClearPoints();
1602 //printf("PRF second loop \n");
1603 ////////////////////////////////////////////////
1604 // Fill the PRF: Second loop over clusters
1605 //////////////////////////////////////////////
1606 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1607 // reject shared clusters on pad row
1608 cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1609 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1611 cl = tracklet->GetClusters(ic);
1614 Short_t *signals = cl->GetSignals(); // signal
1615 Double_t time = cl->GetPadTime(); // time bin
1616 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1617 Float_t padPos = cl->GetPadCol(); // middle pad
1618 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1619 Float_t ycenter = 0.0; // relative center charge
1620 Float_t ymin = 0.0; // relative left charge
1621 Float_t ymax = 0.0; // relative right charge
1623 ////////////////////////////////////////////////////////////////
1624 // Calculate ycenter, ymin and ymax even for two pad clusters
1625 ////////////////////////////////////////////////////////////////
1626 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1627 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1628 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1629 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1630 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1631 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1634 /////////////////////////
1635 // Calibration group
1636 ////////////////////////
1637 Int_t rowcl = cl->GetPadRow(); // row of cluster
1638 Int_t colcl = cl->GetPadCol(); // col of cluster
1639 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1640 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1641 Float_t xcl = cl->GetY(); // y cluster
1642 Float_t qcl = cl->GetQ(); // charge cluster
1643 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1644 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1645 Double_t xdiff = dpad; // reconstructed position constant
1646 Double_t x = dpad; // reconstructed position moved
1647 Float_t ep = pointError; // error of fit
1648 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1649 Float_t signal3 = (Float_t)signals[3]; // signal
1650 Float_t signal2 = (Float_t)signals[2]; // signal
1651 Float_t signal4 = (Float_t)signals[4]; // signal
1652 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1656 /////////////////////
1658 ////////////////////
1660 if(fDebugLevel > 0){
1661 if ( !fDebugStreamer ) {
1663 TDirectory *backup = gDirectory;
1664 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1665 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1670 Float_t y = ycenter;
1671 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1672 "caligroup="<<caligroup<<
1673 "detector="<<fDetectorPreviousTrack<<
1676 "npoints="<<nbclusters<<
1685 "padPosition="<<padPositions[ic]<<
1686 "padPosTracklet="<<padPosTracklet<<
1691 "signal1="<<signal1<<
1692 "signal2="<<signal2<<
1693 "signal3="<<signal3<<
1694 "signal4="<<signal4<<
1695 "signal5="<<signal5<<
1701 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1702 "caligroup="<<caligroup<<
1703 "detector="<<fDetectorPreviousTrack<<
1706 "npoints="<<nbclusters<<
1715 "padPosition="<<padPositions[ic]<<
1716 "padPosTracklet="<<padPosTracklet<<
1721 "signal1="<<signal1<<
1722 "signal2="<<signal2<<
1723 "signal3="<<signal3<<
1724 "signal4="<<signal4<<
1725 "signal5="<<signal5<<
1731 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1732 "caligroup="<<caligroup<<
1733 "detector="<<fDetectorPreviousTrack<<
1736 "npoints="<<nbclusters<<
1745 "padPosition="<<padPositions[ic]<<
1746 "padPosTracklet="<<padPosTracklet<<
1751 "signal1="<<signal1<<
1752 "signal2="<<signal2<<
1753 "signal3="<<signal3<<
1754 "signal4="<<signal4<<
1755 "signal5="<<signal5<<
1761 /////////////////////
1763 /////////////////////
1764 if(nbclusters < fNumberClusters) continue;
1765 if(nbclusters > fNumberClustersf) continue;
1766 if(nb3pc <= 5) continue;
1767 if((time >= 21) || (time < 7)) continue;
1768 if(TMath::Abs(qcl) < 80) continue;
1769 if( TMath::Abs(snp) > 1.) continue;
1772 ////////////////////////
1774 ///////////////////////
1776 if(TMath::Abs(dpad) < 1.5) {
1777 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1778 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1779 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1781 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1782 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1783 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1785 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1786 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1787 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1792 if(TMath::Abs(dpad) < 1.5) {
1793 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1794 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1796 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1797 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1798 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1800 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1801 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1802 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1805 } // second loop over clusters
1810 ///////////////////////////////////////////////////////////////////////////////////////
1811 // Pad row col stuff: see if masked or not
1812 ///////////////////////////////////////////////////////////////////////////////////////
1813 //_____________________________________________________________________________
1814 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1817 // See if we are not near a masked pad
1820 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1824 //_____________________________________________________________________________
1825 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1828 // See if we are not near a masked pad
1831 if (!IsPadOn(detector, col, row)) {
1832 fGoodTracklet = kFALSE;
1836 if (!IsPadOn(detector, col-1, row)) {
1837 fGoodTracklet = kFALSE;
1842 if (!IsPadOn(detector, col+1, row)) {
1843 fGoodTracklet = kFALSE;
1848 //_____________________________________________________________________________
1849 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1852 // Look in the choosen database if the pad is On.
1853 // If no the track will be "not good"
1856 // Get the parameter object
1857 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1859 AliInfo("Could not get calibDB");
1863 if (!cal->IsChamberInstalled(detector) ||
1864 cal->IsChamberMasked(detector) ||
1865 cal->IsPadMasked(detector,col,row)) {
1873 ///////////////////////////////////////////////////////////////////////////////////////
1874 // Calibration groups: calculate the number of groups, localise...
1875 ////////////////////////////////////////////////////////////////////////////////////////
1876 //_____________________________________________________________________________
1877 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1880 // Calculate the calibration group number for i
1883 // Row of the cluster and position in the pad groups
1885 if (fCalibraMode->GetNnZ(i) != 0) {
1886 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1890 // Col of the cluster and position in the pad groups
1892 if (fCalibraMode->GetNnRphi(i) != 0) {
1893 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1896 return posc*fCalibraMode->GetNfragZ(i)+posr;
1899 //____________________________________________________________________________________
1900 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1903 // Calculate the total number of calibration groups
1909 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1911 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1916 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1918 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1923 fCalibraMode->ModePadCalibration(2,i);
1924 fCalibraMode->ModePadFragmentation(0,2,0,i);
1925 fCalibraMode->SetDetChamb2(i);
1926 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1927 fCalibraMode->ModePadCalibration(0,i);
1928 fCalibraMode->ModePadFragmentation(0,0,0,i);
1929 fCalibraMode->SetDetChamb0(i);
1930 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1931 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1935 //_____________________________________________________________________________
1936 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1939 // Set the mode of calibration group in the z direction for the parameter i
1944 fCalibraMode->SetNz(i, Nz);
1947 AliInfo("You have to choose between 0 and 4");
1951 //_____________________________________________________________________________
1952 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1955 // Set the mode of calibration group in the rphi direction for the parameter i
1960 fCalibraMode->SetNrphi(i ,Nrphi);
1963 AliInfo("You have to choose between 0 and 6");
1968 //_____________________________________________________________________________
1969 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1972 // Set the mode of calibration group all together
1974 if(fVector2d == kTRUE) {
1975 AliInfo("Can not work with the vector method");
1978 fCalibraMode->SetAllTogether(i);
1982 //_____________________________________________________________________________
1983 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1986 // Set the mode of calibration group per supermodule
1988 if(fVector2d == kTRUE) {
1989 AliInfo("Can not work with the vector method");
1992 fCalibraMode->SetPerSuperModule(i);
1996 //____________Set the pad calibration variables for the detector_______________
1997 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
2000 // For the detector calcul the first Xbins and set the number of row
2001 // and col pads per calibration groups, the number of calibration
2002 // groups in the detector.
2005 // first Xbins of the detector
2007 fCalibraMode->CalculXBins(detector,0);
2010 fCalibraMode->CalculXBins(detector,1);
2013 fCalibraMode->CalculXBins(detector,2);
2016 // fragmentation of idect
2017 for (Int_t i = 0; i < 3; i++) {
2018 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
2019 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
2020 , (Int_t) GetStack(detector)
2021 , (Int_t) GetSector(detector),i);
2027 //_____________________________________________________________________________
2028 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
2031 // Should be between 0 and 6
2034 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
2035 AliInfo("The number of groups must be between 0 and 6!");
2038 fNgroupprf = numberGroupsPRF;
2042 ///////////////////////////////////////////////////////////////////////////////////////////
2043 // Per tracklet: store or reset the info, fill the histos with the info
2044 //////////////////////////////////////////////////////////////////////////////////////////
2045 //_____________________________________________________________________________
2046 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(const AliTRDcluster *cl,const Double_t dqdl,const Int_t *group,const Int_t row,const Int_t col,const AliTRDcluster *cls)
2049 // Store the infos in fAmpTotal, fPHPlace and fPHValue
2050 // Correct from the gain correction before
2051 // cls is shared cluster if any
2054 //printf("StoreInfoCHPHtrack\n");
2056 // time bin of the cluster not corrected
2057 Int_t time = cl->GetPadTime();
2058 Float_t charge = TMath::Abs(cl->GetQ());
2060 charge += TMath::Abs(cls->GetQ());
2061 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
2064 //printf("Store::time %d, amplitude %f\n",time,dqdl);
2066 //Correct for the gain coefficient used in the database for reconstruction
2067 Float_t correctthegain = 1.0;
2068 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
2069 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
2070 Float_t correction = 1.0;
2071 Float_t normalisation = 6.67;
2072 // we divide with gain in AliTRDclusterizer::Transform...
2073 if( correctthegain > 0 ) normalisation /= correctthegain;
2076 // take dd/dl corrected from the angle of the track
2077 correction = dqdl / (normalisation);
2080 // Fill the fAmpTotal with the charge
2082 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
2083 //printf("Store::group %d, amplitude %f\n",group[0],correction);
2084 fAmpTotal[(Int_t) group[0]] += correction;
2088 // Fill the fPHPlace and value
2090 fPHPlace[time] = group[1];
2091 fPHValue[time] = charge;
2095 //____________Offine tracking in the AliTRDtracker_____________________________
2096 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
2099 // Reset values per tracklet
2102 //Reset good tracklet
2103 fGoodTracklet = kTRUE;
2105 // Reset the fPHValue
2107 //Reset the fPHValue and fPHPlace
2108 for (Int_t k = 0; k < fTimeMax; k++) {
2114 // Reset the fAmpTotal where we put value
2116 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2121 //____________Offine tracking in the AliTRDtracker_____________________________
2122 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
2125 // For the offline tracking or mcm tracklets
2126 // This function will be called in the functions UpdateHistogram...
2127 // to fill the info of a track for the relativ gain calibration
2130 Int_t nb = 0; // Nombre de zones traversees
2131 Int_t fd = -1; // Premiere zone non nulle
2132 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
2134 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2136 if(nbclusters < fNumberClusters) return;
2137 if(nbclusters > fNumberClustersf) return;
2140 // Normalize with the number of clusters
2141 Double_t normalizeCst = fRelativeScale;
2142 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
2144 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
2146 // See if the track goes through different zones
2147 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2148 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
2149 if (fAmpTotal[k] > 0.0) {
2150 totalcharge += fAmpTotal[k];
2158 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
2164 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2166 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2167 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2170 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2174 if ((fAmpTotal[fd] > 0.0) &&
2175 (fAmpTotal[fd+1] > 0.0)) {
2176 // One of the two very big
2177 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2179 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2180 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2183 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2186 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2188 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2190 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2191 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2194 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2197 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2200 if (fCalibraMode->GetNfragZ(0) > 1) {
2201 if (fAmpTotal[fd] > 0.0) {
2202 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2203 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2204 // One of the two very big
2205 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2207 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2208 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2211 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2214 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2216 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2218 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2219 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2222 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2224 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2235 //____________Offine tracking in the AliTRDtracker_____________________________
2236 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2239 // For the offline tracking or mcm tracklets
2240 // This function will be called in the functions UpdateHistogram...
2241 // to fill the info of a track for the drift velocity calibration
2244 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2245 Int_t fd1 = -1; // Premiere zone non nulle
2246 Int_t fd2 = -1; // Deuxieme zone non nulle
2247 Int_t k1 = -1; // Debut de la premiere zone
2248 Int_t k2 = -1; // Debut de la seconde zone
2249 Int_t nbclusters = 0; // number of clusters
2253 // See if the track goes through different zones
2254 for (Int_t k = 0; k < fTimeMax; k++) {
2255 if (fPHValue[k] > 0.0) {
2261 if (fPHPlace[k] != fd1) {
2267 if (fPHPlace[k] != fd2) {
2274 // See if noise before and after
2275 if(fMaxCluster > 0) {
2276 if(fPHValue[0] > fMaxCluster) return;
2277 if(fTimeMax > fNbMaxCluster) {
2278 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2279 if(fPHValue[k] > fMaxCluster) return;
2284 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2286 if(nbclusters < fNumberClusters) return;
2287 if(nbclusters > fNumberClustersf) return;
2293 for (Int_t i = 0; i < fTimeMax; i++) {
2295 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2297 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2299 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2302 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2304 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2310 if ((fd1 == fd2+1) ||
2312 // One of the two fast all the think
2313 if (k2 > (k1+fDifference)) {
2314 //we choose to fill the fd1 with all the values
2316 for (Int_t i = 0; i < fTimeMax; i++) {
2318 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2320 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2324 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2326 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2331 if ((k2+fDifference) < fTimeMax) {
2332 //we choose to fill the fd2 with all the values
2334 for (Int_t i = 0; i < fTimeMax; i++) {
2336 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2338 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2342 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2344 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2350 // Two zones voisines sinon rien!
2351 if (fCalibraMode->GetNfragZ(1) > 1) {
2353 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2354 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2355 // One of the two fast all the think
2356 if (k2 > (k1+fDifference)) {
2357 //we choose to fill the fd1 with all the values
2359 for (Int_t i = 0; i < fTimeMax; i++) {
2361 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2363 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2367 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2369 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2374 if ((k2+fDifference) < fTimeMax) {
2375 //we choose to fill the fd2 with all the values
2377 for (Int_t i = 0; i < fTimeMax; i++) {
2379 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2381 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2385 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2387 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2394 // Two zones voisines sinon rien!
2396 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2397 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2398 // One of the two fast all the think
2399 if (k2 > (k1 + fDifference)) {
2400 //we choose to fill the fd1 with all the values
2402 for (Int_t i = 0; i < fTimeMax; i++) {
2404 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2406 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2410 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2412 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2417 if ((k2+fDifference) < fTimeMax) {
2418 //we choose to fill the fd2 with all the values
2420 for (Int_t i = 0; i < fTimeMax; i++) {
2422 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2424 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2428 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2430 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2442 //////////////////////////////////////////////////////////////////////////////////////////
2443 // DAQ process functions
2444 /////////////////////////////////////////////////////////////////////////////////////////
2445 //_____________________________________________________________________
2446 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2449 // Event Processing loop - AliTRDrawStreamBase
2450 // TestBeam 2007 version
2451 // 0 timebin problem
2454 // Same algorithm as TestBeam but different reader
2457 rawStream->SetSharedPadReadout(kFALSE);
2459 Int_t withInput = 1;
2461 Double_t phvalue[16][144][36];
2462 for(Int_t k = 0; k < 36; k++){
2463 for(Int_t j = 0; j < 16; j++){
2464 for(Int_t c = 0; c < 144; c++){
2465 phvalue[j][c][k] = 0.0;
2470 fDetectorPreviousTrack = -1;
2473 Int_t nbtimebin = 0;
2474 Int_t baseline = 10;
2475 //printf("------------Detector\n");
2481 while (rawStream->Next()) {
2483 Int_t idetector = rawStream->GetDet(); // current detector
2484 Int_t imcm = rawStream->GetMCM(); // current MCM
2485 Int_t irob = rawStream->GetROB(); // current ROB
2487 //printf("Detector %d\n",idetector);
2489 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2492 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2496 for(Int_t k = 0; k < 36; k++){
2497 for(Int_t j = 0; j < 16; j++){
2498 for(Int_t c = 0; c < 144; c++){
2499 phvalue[j][c][k] = 0.0;
2505 fDetectorPreviousTrack = idetector;
2506 fMCMPrevious = imcm;
2507 fROBPrevious = irob;
2509 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2510 if(nbtimebin == 0) return 0;
2511 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2512 fTimeMax = nbtimebin;
2514 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2515 fNumberClustersf = fTimeMax;
2516 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2519 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2520 Int_t col = rawStream->GetCol();
2521 Int_t row = rawStream->GetRow();
2524 // printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2527 for(Int_t itime = 0; itime < nbtimebin; itime++){
2528 phvalue[row][col][itime] = signal[itime]-baseline;
2532 // fill the last one
2533 if(fDetectorPreviousTrack != -1){
2536 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2539 for(Int_t k = 0; k < 36; k++){
2540 for(Int_t j = 0; j < 16; j++){
2541 for(Int_t c = 0; c < 144; c++){
2542 phvalue[j][c][k] = 0.0;
2551 while (rawStream->Next()) { //iddetecte
2553 Int_t idetector = rawStream->GetDet(); // current detector
2554 Int_t imcm = rawStream->GetMCM(); // current MCM
2555 Int_t irob = rawStream->GetROB(); // current ROB
2557 //printf("Detector %d\n",idetector);
2559 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2562 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2565 for(Int_t k = 0; k < 36; k++){
2566 for(Int_t j = 0; j < 16; j++){
2567 for(Int_t c = 0; c < 144; c++){
2568 phvalue[j][c][k] = 0.0;
2574 fDetectorPreviousTrack = idetector;
2575 fMCMPrevious = imcm;
2576 fROBPrevious = irob;
2578 //baseline = rawStream->GetCommonAdditive(); // common baseline
2580 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2581 fNumberClustersf = fTimeMax;
2582 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2583 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2584 Int_t col = rawStream->GetCol();
2585 Int_t row = rawStream->GetRow();
2588 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2590 for(Int_t itime = 0; itime < fTimeMax; itime++){
2591 phvalue[row][col][itime] = signal[itime]-baseline;
2592 /*if(phvalue[row][col][itime] >= 20) {
2593 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2603 // fill the last one
2604 if(fDetectorPreviousTrack != -1){
2607 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2610 for(Int_t k = 0; k < 36; k++){
2611 for(Int_t j = 0; j < 16; j++){
2612 for(Int_t c = 0; c < 144; c++){
2613 phvalue[j][c][k] = 0.0;
2623 //_____________________________________________________________________
2624 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2627 // Event processing loop - AliRawReader
2628 // Testbeam 2007 version
2631 AliTRDrawStreamBase rawStream(rawReader);
2633 rawReader->Select("TRD");
2635 return ProcessEventDAQ(&rawStream, nocheck);
2638 //_________________________________________________________________________
2639 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2641 const eventHeaderStruct *event,
2644 const eventHeaderStruct* /*event*/,
2651 // process date event
2652 // Testbeam 2007 version
2655 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2656 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2660 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2665 //_____________________________________________________________________
2666 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ2(AliRawReader *rawReader)
2669 // Event Processing loop - AliTRDrawFastStream
2671 // 0 timebin problem
2674 // Same algorithm as TestBeam but different reader
2677 // AliTRDrawFastStream *rawStream = AliTRDrawFastStream::GetRawStream(rawReader);
2678 AliTRDrawFastStream *rawStream = new AliTRDrawFastStream(rawReader);
2679 rawStream->SetNoErrorWarning();
2680 rawStream->SetSharedPadReadout(kFALSE);
2682 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2683 digitsManager->CreateArrays();
2685 Int_t withInput = 1;
2687 Double_t phvalue[16][144][36];
2688 for(Int_t k = 0; k < 36; k++){
2689 for(Int_t j = 0; j < 16; j++){
2690 for(Int_t c = 0; c < 144; c++){
2691 phvalue[j][c][k] = 0.0;
2696 fDetectorPreviousTrack = -1;
2700 Int_t nbtimebin = 0;
2701 Int_t baseline = 10;
2707 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2709 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2710 // printf("there is ADC data on this chamber!\n");
2712 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2713 if (digits->HasData()) { //array
2715 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2716 if (indexes->IsAllocated() == kFALSE) {
2717 AliError("Indexes do not exist!");
2722 indexes->ResetCounters();
2724 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2725 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2726 //while (rawStream->Next()) {
2728 Int_t idetector = det; // current detector
2729 //Int_t imcm = rawStream->GetMCM(); // current MCM
2730 //Int_t irob = rawStream->GetROB(); // current ROB
2733 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2735 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2738 for(Int_t k = 0; k < 36; k++){
2739 for(Int_t j = 0; j < 16; j++){
2740 for(Int_t c = 0; c < 144; c++){
2741 phvalue[j][c][k] = 0.0;
2747 fDetectorPreviousTrack = idetector;
2748 //fMCMPrevious = imcm;
2749 //fROBPrevious = irob;
2751 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2752 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2753 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2754 baseline = digitParam->GetADCbaseline(det); // baseline
2756 if(nbtimebin == 0) return 0;
2757 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2758 fTimeMax = nbtimebin;
2760 fNumberClustersf = fTimeMax;
2761 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2764 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2765 // phvalue[row][col][itime] = signal[itime]-baseline;
2766 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2767 /*if(phvalue[iRow][iCol][itime] >= 20) {
2768 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2772 (Short_t)(digits->GetData(iRow,iCol,itime)),
2779 // fill the last one
2780 if(fDetectorPreviousTrack != -1){
2783 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2784 // printf("\n ---> withinput %d\n\n",withInput);
2786 for(Int_t k = 0; k < 36; k++){
2787 for(Int_t j = 0; j < 16; j++){
2788 for(Int_t c = 0; c < 144; c++){
2789 phvalue[j][c][k] = 0.0;
2797 digitsManager->ClearArrays(det);
2799 delete digitsManager;
2804 //_____________________________________________________________________
2805 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ3(AliRawReader *rawReader)
2808 // Event Processing loop - AliTRDrawStream
2810 // 0 timebin problem
2813 // Same algorithm as TestBeam but different reader
2816 // AliTRDrawFastStream *rawStream = AliTRDrawFastStream::GetRawStream(rawReader);
2817 AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2818 rawStream->SetNoErrorWarning();
2819 rawStream->SetSharedPadReadout(kFALSE);
2821 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2822 digitsManager->CreateArrays();
2824 Int_t withInput = 1;
2826 Double_t phvalue[16][144][36];
2827 for(Int_t k = 0; k < 36; k++){
2828 for(Int_t j = 0; j < 16; j++){
2829 for(Int_t c = 0; c < 144; c++){
2830 phvalue[j][c][k] = 0.0;
2835 fDetectorPreviousTrack = -1;
2839 Int_t nbtimebin = 0;
2840 Int_t baseline = 10;
2846 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2848 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2849 // printf("there is ADC data on this chamber!\n");
2851 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2852 if (digits->HasData()) { //array
2854 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2855 if (indexes->IsAllocated() == kFALSE) {
2856 AliError("Indexes do not exist!");
2861 indexes->ResetCounters();
2863 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2864 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2865 //while (rawStream->Next()) {
2867 Int_t idetector = det; // current detector
2868 //Int_t imcm = rawStream->GetMCM(); // current MCM
2869 //Int_t irob = rawStream->GetROB(); // current ROB
2872 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2874 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2877 for(Int_t k = 0; k < 36; k++){
2878 for(Int_t j = 0; j < 16; j++){
2879 for(Int_t c = 0; c < 144; c++){
2880 phvalue[j][c][k] = 0.0;
2886 fDetectorPreviousTrack = idetector;
2887 //fMCMPrevious = imcm;
2888 //fROBPrevious = irob;
2890 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2891 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2892 nbtimebin = digitParam->GetNTimeBins(det); // number of time bins read from data
2893 baseline = digitParam->GetADCbaseline(det); // baseline
2895 if(nbtimebin == 0) return 0;
2896 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2897 fTimeMax = nbtimebin;
2899 fNumberClustersf = fTimeMax;
2900 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2903 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2904 // phvalue[row][col][itime] = signal[itime]-baseline;
2905 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2906 /*if(phvalue[iRow][iCol][itime] >= 20) {
2907 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2911 (Short_t)(digits->GetData(iRow,iCol,itime)),
2918 // fill the last one
2919 if(fDetectorPreviousTrack != -1){
2922 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2923 // printf("\n ---> withinput %d\n\n",withInput);
2925 for(Int_t k = 0; k < 36; k++){
2926 for(Int_t j = 0; j < 16; j++){
2927 for(Int_t c = 0; c < 144; c++){
2928 phvalue[j][c][k] = 0.0;
2936 digitsManager->ClearArrays(det);
2938 delete digitsManager;
2943 //_____________________________________________________________________
2944 //////////////////////////////////////////////////////////////////////////////
2945 // Routine inside the DAQ process
2946 /////////////////////////////////////////////////////////////////////////////
2947 //_______________________________________________________________________
2948 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2951 // Look for the maximum by collapsing over the time
2952 // Sum over four pad col and two pad row
2958 Int_t idect = fDetectorPreviousTrack;
2959 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2961 for(Int_t tb = 0; tb < 36; tb++){
2965 //fGoodTracklet = kTRUE;
2966 //fDetectorPreviousTrack = 0;
2969 ///////////////////////////
2971 /////////////////////////
2975 Double_t integralMax = -1;
2977 for (Int_t ir = 1; ir <= 15; ir++)
2979 for (Int_t ic = 2; ic <= 142; ic++)
2981 Double_t integral = 0;
2982 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2984 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2986 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2987 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2990 for(Int_t tb = 0; tb< fTimeMax; tb++){
2991 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2996 if (integralMax < integral)
3000 integralMax = integral;
3002 } // check max integral
3006 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
3007 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
3012 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
3016 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
3017 //if(!fGoodTracklet) used = 1;;
3019 // /////////////////////////////////////////////////////
3020 // sum ober 2 row and 4 pad cols for each time bins
3021 // ////////////////////////////////////////////////////
3025 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
3027 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
3029 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
3030 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
3032 for(Int_t it = 0; it < fTimeMax; it++){
3033 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
3040 Double_t sumcharge = 0.0;
3041 for(Int_t it = 0; it < fTimeMax; it++){
3042 sumcharge += sum[it];
3043 if(sum[it] > fThresholdClustersDAQ) nbcl++;
3047 /////////////////////////////////////////////////////////
3049 ////////////////////////////////////////////////////////
3050 if(fDebugLevel > 0){
3051 if ( !fDebugStreamer ) {
3053 TDirectory *backup = gDirectory;
3054 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
3055 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3058 Double_t amph0 = sum[0];
3059 Double_t amphlast = sum[fTimeMax-1];
3060 Double_t rms = TMath::RMS(fTimeMax,sum);
3061 Int_t goodtracklet = (Int_t) fGoodTracklet;
3062 for(Int_t it = 0; it < fTimeMax; it++){
3063 Double_t clustera = sum[it];
3065 (* fDebugStreamer) << "FillDAQa"<<
3066 "ampTotal="<<sumcharge<<
3069 "detector="<<idect<<
3071 "amphlast="<<amphlast<<
3072 "goodtracklet="<<goodtracklet<<
3073 "clustera="<<clustera<<
3081 ////////////////////////////////////////////////////////
3083 ///////////////////////////////////////////////////////
3084 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
3085 if(sum[0] > 100.0) used = 1;
3086 if(nbcl < fNumberClusters) used = 1;
3087 if(nbcl > fNumberClustersf) used = 1;
3089 //if(fDetectorPreviousTrack == 15){
3090 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
3092 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
3094 for(Int_t it = 0; it < fTimeMax; it++){
3095 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
3097 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
3099 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
3101 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
3106 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
3108 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
3115 //____________Online trackling in AliTRDtrigger________________________________
3116 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
3120 // Fill a simple average pulse height
3124 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
3130 //____________Write_____________________________________________________
3131 //_____________________________________________________________________
3132 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
3135 // Write infos to file
3139 if ( fDebugStreamer ) {
3140 delete fDebugStreamer;
3141 fDebugStreamer = 0x0;
3144 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
3149 ,fNumberUsedPh[1]));
3151 TDirectory *backup = gDirectory;
3157 option = "recreate";
3159 TFile f(filename,option.Data());
3161 TStopwatch stopwatch;
3164 f.WriteTObject(fCalibraVector);
3169 f.WriteTObject(fCH2d);
3174 f.WriteTObject(fPH2d);
3179 f.WriteTObject(fPRF2d);
3182 if(fLinearFitterOn){
3183 if(fLinearFitterDebugOn) AnalyseLinearFitter();
3184 f.WriteTObject(fLinearVdriftFit);
3189 if ( backup ) backup->cd();
3191 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
3192 ,stopwatch.RealTime(),stopwatch.CpuTime()));
3194 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3196 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3197 //___________________________________________probe the histos__________________________________________________
3198 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
3201 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
3202 // debug mode with 2 for TH2I and 3 for TProfile2D
3203 // It gives a pointer to a Double_t[7] with the info following...
3204 // [0] : number of calibration groups with entries
3205 // [1] : minimal number of entries found
3206 // [2] : calibration group number of the min
3207 // [3] : maximal number of entries found
3208 // [4] : calibration group number of the max
3209 // [5] : mean number of entries found
3210 // [6] : mean relative error
3213 Double_t *info = new Double_t[7];
3215 // Number of Xbins (detectors or groups of pads)
3216 Int_t nbins = h->GetNbinsY(); //number of calibration groups
3217 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
3220 Double_t nbwe = 0; //number of calibration groups with entries
3221 Double_t minentries = 0; //minimal number of entries found
3222 Double_t maxentries = 0; //maximal number of entries found
3223 Double_t placemin = 0; //calibration group number of the min
3224 Double_t placemax = -1; //calibration group number of the max
3225 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
3226 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
3228 Double_t counter = 0;
3231 TH1F *nbEntries = 0x0;//distribution of the number of entries
3232 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
3233 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
3235 // Beginning of the loop over the calibration groups
3236 for (Int_t idect = 0; idect < nbins; idect++) {
3238 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
3239 projch->SetDirectory(0);
3241 // Number of entries for this calibration group
3242 Double_t nentries = 0.0;
3244 for (Int_t k = 0; k < nxbins; k++) {
3245 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
3249 for (Int_t k = 0; k < nxbins; k++) {
3250 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3251 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3252 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3261 if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
3262 nbEntries->SetDirectory(0);
3263 nbEntries->Fill(nentries);
3264 if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
3265 nbEntriesPerGroup->SetDirectory(0);
3266 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3267 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));
3268 nbEntriesPerSp->SetDirectory(0);
3269 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3274 if(nentries > maxentries){
3275 maxentries = nentries;
3279 minentries = nentries;
3281 if(nentries < minentries){
3282 minentries = nentries;
3288 meanstats += nentries;
3290 }//calibration groups loop
3292 if(nbwe > 0) meanstats /= nbwe;
3293 if(counter > 0) meanrelativerror /= counter;
3295 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3296 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3297 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3298 AliInfo(Form("The mean number of entries is %f",meanstats));
3299 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3302 info[1] = minentries;
3304 info[3] = maxentries;
3306 info[5] = meanstats;
3307 info[6] = meanrelativerror;
3309 if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
3310 gStyle->SetPalette(1);
3311 gStyle->SetOptStat(1111);
3312 gStyle->SetPadBorderMode(0);
3313 gStyle->SetCanvasColor(10);
3314 gStyle->SetPadLeftMargin(0.13);
3315 gStyle->SetPadRightMargin(0.01);
3316 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3319 nbEntries->Draw("");
3321 nbEntriesPerSp->SetStats(0);
3322 nbEntriesPerSp->Draw("");
3323 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3325 nbEntriesPerGroup->SetStats(0);
3326 nbEntriesPerGroup->Draw("");
3332 //____________________________________________________________________________
3333 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3336 // Return a Int_t[4] with:
3337 // 0 Mean number of entries
3338 // 1 median of number of entries
3339 // 2 rms of number of entries
3340 // 3 number of group with entries
3343 Double_t *stat = new Double_t[4];
3346 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3348 Double_t *weight = new Double_t[nbofgroups];
3349 Double_t *nonul = new Double_t[nbofgroups];
3351 for(Int_t k = 0; k < nbofgroups; k++){
3352 if(fEntriesCH[k] > 0) {
3354 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3357 else weight[k] = 0.0;
3359 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3360 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3361 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3369 //____________________________________________________________________________
3370 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3373 // Return a Int_t[4] with:
3374 // 0 Mean number of entries
3375 // 1 median of number of entries
3376 // 2 rms of number of entries
3377 // 3 number of group with entries
3380 Double_t *stat = new Double_t[4];
3383 Int_t nbofgroups = 540;
3384 Double_t *weight = new Double_t[nbofgroups];
3385 Int_t *nonul = new Int_t[nbofgroups];
3387 for(Int_t k = 0; k < nbofgroups; k++){
3388 if(fEntriesLinearFitter[k] > 0) {
3390 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3393 else weight[k] = 0.0;
3395 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3396 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3397 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3405 //////////////////////////////////////////////////////////////////////////////////////
3407 //////////////////////////////////////////////////////////////////////////////////////
3408 //_____________________________________________________________________________
3409 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3412 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3413 // If fNgroupprf is zero then no binning in tnp
3417 name += fCalibraMode->GetNz(2);
3419 name += fCalibraMode->GetNrphi(2);
3423 if(fNgroupprf != 0){
3425 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3426 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3427 fPRF2d->SetYTitle("Det/pad groups");
3428 fPRF2d->SetXTitle("Position x/W [pad width units]");
3429 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3430 fPRF2d->SetStats(0);
3433 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3434 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3435 fPRF2d->SetYTitle("Det/pad groups");
3436 fPRF2d->SetXTitle("Position x/W [pad width units]");
3437 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3438 fPRF2d->SetStats(0);
3443 //_____________________________________________________________________________
3444 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3447 // Create the 2D histos
3450 TString name("Ver");
3451 name += fVersionVdriftUsed;
3453 name += fSubVersionVdriftUsed;
3455 name += fCalibraMode->GetNz(1);
3457 name += fCalibraMode->GetNrphi(1);
3459 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3460 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3462 fPH2d->SetYTitle("Det/pad groups");
3463 fPH2d->SetXTitle("time [#mus]");
3464 fPH2d->SetZTitle("<PH> [a.u.]");
3468 //_____________________________________________________________________________
3469 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3472 // Create the 2D histos
3475 TString name("Ver");
3476 name += fVersionGainUsed;
3478 name += fSubVersionGainUsed;
3480 name += fCalibraMode->GetNz(0);
3482 name += fCalibraMode->GetNrphi(0);
3484 fCH2d = new TH2I("CH2d",(const Char_t *) name
3485 ,fNumberBinCharge,0,300,nn,0,nn);
3486 fCH2d->SetYTitle("Det/pad groups");
3487 fCH2d->SetXTitle("charge deposit [a.u]");
3488 fCH2d->SetZTitle("counts");
3493 //////////////////////////////////////////////////////////////////////////////////
3494 // Set relative scale
3495 /////////////////////////////////////////////////////////////////////////////////
3496 //_____________________________________________________________________________
3497 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3500 // Set the factor that will divide the deposited charge
3501 // to fit in the histo range [0,300]
3504 if (RelativeScale > 0.0) {
3505 fRelativeScale = RelativeScale;
3508 AliInfo("RelativeScale must be strict positif!");
3512 //////////////////////////////////////////////////////////////////////////////////
3513 // Quick way to fill a histo
3514 //////////////////////////////////////////////////////////////////////////////////
3515 //_____________________________________________________________________
3516 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3519 // FillCH2d: Marian style
3522 //skip simply the value out of range
3523 if((y>=300.0) || (y<0.0)) return;
3525 //Calcul the y place
3526 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3527 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3530 fCH2d->GetArray()[place]++;
3534 //////////////////////////////////////////////////////////////////////////////////
3535 // Geometrical functions
3536 ///////////////////////////////////////////////////////////////////////////////////
3537 //_____________________________________________________________________________
3538 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3541 // Reconstruct the layer number from the detector number
3544 return ((Int_t) (d % 6));
3548 //_____________________________________________________________________________
3549 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3552 // Reconstruct the stack number from the detector number
3554 const Int_t kNlayer = 6;
3556 return ((Int_t) (d % 30) / kNlayer);
3560 //_____________________________________________________________________________
3561 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3564 // Reconstruct the sector number from the detector number
3568 return ((Int_t) (d / fg));
3571 ///////////////////////////////////////////////////////////////////////////////////
3572 // Getter functions for DAQ of the CH2d and the PH2d
3573 //////////////////////////////////////////////////////////////////////////////////
3574 //_____________________________________________________________________
3575 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3578 // return pointer to fPH2d TProfile2D
3579 // create a new TProfile2D if it doesn't exist allready
3585 fTimeMax = nbtimebin;
3586 fSf = samplefrequency;
3592 //_____________________________________________________________________
3593 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3596 // return pointer to fCH2d TH2I
3597 // create a new TH2I if it doesn't exist allready
3606 ////////////////////////////////////////////////////////////////////////////////////////////
3607 // Drift velocity calibration
3608 ///////////////////////////////////////////////////////////////////////////////////////////
3609 //_____________________________________________________________________
3610 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3613 // return pointer to TLinearFitter Calibration
3614 // if force is true create a new TLinearFitter if it doesn't exist allready
3617 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3618 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3621 // if we are forced and TLinearFitter doesn't yet exist create it
3623 // new TLinearFitter
3624 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3625 fLinearFitterArray.AddAt(linearfitter,detector);
3626 return linearfitter;
3629 //____________________________________________________________________________
3630 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3633 // Analyse array of linear fitter because can not be written
3634 // Store two arrays: one with the param the other one with the error param + number of entries
3637 for(Int_t k = 0; k < 540; k++){
3638 TLinearFitter *linearfitter = GetLinearFitter(k);
3639 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3640 TVectorD *par = new TVectorD(2);
3641 TVectorD pare = TVectorD(2);
3642 TVectorD *parE = new TVectorD(3);
3643 linearfitter->Eval();
3644 linearfitter->GetParameters(*par);
3645 linearfitter->GetErrors(pare);
3646 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3647 (*parE)[0] = pare[0]*ppointError;
3648 (*parE)[1] = pare[1]*ppointError;
3649 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3650 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3651 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);