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/AliTRDCalDet.h"
73 #include "AliTRDrawFastStream.h"
74 #include "AliTRDdigitsManager.h"
75 #include "AliTRDdigitsParam.h"
76 #include "AliTRDSignalIndex.h"
77 #include "AliTRDarrayADC.h"
84 ClassImp(AliTRDCalibraFillHisto)
86 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
87 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
89 //_____________singleton implementation_________________________________________________
90 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
93 // Singleton implementation
96 if (fgTerminated != kFALSE) {
100 if (fgInstance == 0) {
101 fgInstance = new AliTRDCalibraFillHisto();
108 //______________________________________________________________________________________
109 void AliTRDCalibraFillHisto::Terminate()
112 // Singleton implementation
113 // Deletes the instance of this class
116 fgTerminated = kTRUE;
118 if (fgInstance != 0) {
125 //______________________________________________________________________________________
126 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
135 ,fLinearFitterOn(kFALSE)
136 ,fLinearFitterDebugOn(kFALSE)
138 ,fThresholdClusterPRF2(15.0)
139 ,fLimitChargeIntegration(kFALSE)
140 ,fFillWithZero(kFALSE)
141 ,fNormalizeNbOfCluster(kFALSE)
144 ,fCalibraMode(new AliTRDCalibraMode())
147 ,fDetectorPreviousTrack(-1)
151 ,fNumberClustersf(30)
152 ,fNumberClustersProcent(0.5)
153 ,fThresholdClustersDAQ(120.0)
161 ,fNumberBinCharge(50)
167 ,fGoodTracklet(kTRUE)
168 ,fLinearFitterTracklet(0x0)
170 ,fEntriesLinearFitter(0x0)
175 ,fLinearFitterArray(540)
176 ,fLinearVdriftFit(0x0)
181 // Default constructor
185 // Init some default values
188 fNumberUsedCh[0] = 0;
189 fNumberUsedCh[1] = 0;
190 fNumberUsedPh[0] = 0;
191 fNumberUsedPh[1] = 0;
193 fGeo = new AliTRDgeometry();
197 //______________________________________________________________________________________
198 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
204 ,fPRF2dOn(c.fPRF2dOn)
205 ,fHisto2d(c.fHisto2d)
206 ,fVector2d(c.fVector2d)
207 ,fLinearFitterOn(c.fLinearFitterOn)
208 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
209 ,fRelativeScale(c.fRelativeScale)
210 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
211 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
212 ,fFillWithZero(c.fFillWithZero)
213 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
214 ,fMaxCluster(c.fMaxCluster)
215 ,fNbMaxCluster(c.fNbMaxCluster)
218 ,fDebugLevel(c.fDebugLevel)
219 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
220 ,fMCMPrevious(c.fMCMPrevious)
221 ,fROBPrevious(c.fROBPrevious)
222 ,fNumberClusters(c.fNumberClusters)
223 ,fNumberClustersf(c.fNumberClustersf)
224 ,fNumberClustersProcent(c.fNumberClustersProcent)
225 ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
226 ,fNumberRowDAQ(c.fNumberRowDAQ)
227 ,fNumberColDAQ(c.fNumberColDAQ)
228 ,fProcent(c.fProcent)
229 ,fDifference(c.fDifference)
230 ,fNumberTrack(c.fNumberTrack)
231 ,fTimeMax(c.fTimeMax)
233 ,fNumberBinCharge(c.fNumberBinCharge)
234 ,fNumberBinPRF(c.fNumberBinPRF)
235 ,fNgroupprf(c.fNgroupprf)
239 ,fGoodTracklet(c.fGoodTracklet)
240 ,fLinearFitterTracklet(0x0)
242 ,fEntriesLinearFitter(0x0)
247 ,fLinearFitterArray(540)
248 ,fLinearVdriftFit(0x0)
255 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
256 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
258 fPH2d = new TProfile2D(*c.fPH2d);
259 fPH2d->SetDirectory(0);
262 fPRF2d = new TProfile2D(*c.fPRF2d);
263 fPRF2d->SetDirectory(0);
266 fCH2d = new TH2I(*c.fCH2d);
267 fCH2d->SetDirectory(0);
269 if(c.fLinearVdriftFit){
270 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
273 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
274 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
279 fGeo = new AliTRDgeometry();
282 //____________________________________________________________________________________
283 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
286 // AliTRDCalibraFillHisto destructor
290 if ( fDebugStreamer ) delete fDebugStreamer;
292 if ( fCalDetGain ) delete fCalDetGain;
293 if ( fCalROCGain ) delete fCalROCGain;
295 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
299 delete [] fEntriesCH;
300 delete [] fEntriesLinearFitter;
303 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
304 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
307 if(fLinearVdriftFit) delete fLinearVdriftFit;
313 //_____________________________________________________________________________
314 void AliTRDCalibraFillHisto::Destroy()
325 //_____________________________________________________________________________
326 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
329 // Delete DebugStreamer
332 if ( fDebugStreamer ) delete fDebugStreamer;
333 fDebugStreamer = 0x0;
336 //_____________________________________________________________________________
337 void AliTRDCalibraFillHisto::ClearHistos()
357 //////////////////////////////////////////////////////////////////////////////////
358 // calibration with AliTRDtrackV1: Init, Update
359 //////////////////////////////////////////////////////////////////////////////////
360 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
361 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
364 // Init the histograms and stuff to be filled
369 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
371 AliInfo("Could not get calibDB");
374 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
376 AliInfo("Could not get CommonParam");
381 if(nboftimebin > 0) fTimeMax = nboftimebin;
382 else fTimeMax = cal->GetNumberOfTimeBinsDCS();
383 if(fTimeMax <= 0) fTimeMax = 30;
384 printf("////////////////////////////////////////////\n");
385 printf("Number of time bins in calibration component %d\n",fTimeMax);
386 printf("////////////////////////////////////////////\n");
387 fSf = parCom->GetSamplingFrequency();
388 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
389 else fRelativeScale = 1.18;
390 fNumberClustersf = fTimeMax;
391 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
393 // Init linear fitter
394 if(!fLinearFitterTracklet) {
395 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
396 fLinearFitterTracklet->StoreData(kTRUE);
399 //calib object from database used for reconstruction
401 fCalDetGain->~AliTRDCalDet();
402 new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
403 }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
405 // Calcul Xbins Chambd0, Chamb2
406 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
407 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
408 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
410 // If vector method On initialised all the stuff
412 fCalibraVector = new AliTRDCalibraVector();
413 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
414 fCalibraVector->SetTimeMax(fTimeMax);
415 if(fNgroupprf != 0) {
416 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
417 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
420 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
421 fCalibraVector->SetPRFRange(1.5);
423 for(Int_t k = 0; k < 3; k++){
424 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
425 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
427 fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
428 fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
429 fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
430 fCalibraVector->SetNbGroupPRF(fNgroupprf);
433 // Create the 2D histos corresponding to the pad groupCalibration mode
436 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
437 ,fCalibraMode->GetNz(0)
438 ,fCalibraMode->GetNrphi(0)));
440 // Create the 2D histo
445 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
446 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
450 fEntriesCH = new Int_t[ntotal0];
451 for(Int_t k = 0; k < ntotal0; k++){
458 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
459 ,fCalibraMode->GetNz(1)
460 ,fCalibraMode->GetNrphi(1)));
462 // Create the 2D histo
467 fPHPlace = new Short_t[fTimeMax];
468 for (Int_t k = 0; k < fTimeMax; k++) {
471 fPHValue = new Float_t[fTimeMax];
472 for (Int_t k = 0; k < fTimeMax; k++) {
476 if (fLinearFitterOn) {
477 if(fLinearFitterDebugOn) {
478 fLinearFitterArray.SetName("ArrayLinearFitters");
479 fEntriesLinearFitter = new Int_t[540];
480 for(Int_t k = 0; k < 540; k++){
481 fEntriesLinearFitter[k] = 0;
484 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
489 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
490 ,fCalibraMode->GetNz(2)
491 ,fCalibraMode->GetNrphi(2)));
492 // Create the 2D histo
494 CreatePRF2d(ntotal2);
501 //____________Offline tracking in the AliTRDtracker____________________________
502 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(const AliTRDtrack *t)
505 // Use AliTRDtrack for the calibration
509 AliTRDcluster *cl = 0x0; // pointeur to cluster
510 Int_t index0 = 0; // index of the first cluster in the current chamber
511 Int_t index1 = 0; // index of the current cluster in the current chamber
513 //////////////////////////////
514 // loop over the clusters
515 ///////////////////////////////
516 while((cl = t->GetCluster(index1))){
518 /////////////////////////////////////////////////////////////////////////
519 // Fill the infos for the previous clusters if not the same detector
520 ////////////////////////////////////////////////////////////////////////
521 Int_t detector = cl->GetDetector();
522 if ((detector != fDetectorPreviousTrack) &&
523 (index0 != index1)) {
527 //If the same track, then look if the previous detector is in
528 //the same plane, if yes: not a good track
529 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
533 // Fill only if the track doesn't touch a masked pad or doesn't
536 // drift velocity unables to cut bad tracklets
537 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
541 FillTheInfoOfTheTrackCH(index1-index0);
546 FillTheInfoOfTheTrackPH();
549 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
552 } // if a good tracklet
555 ResetfVariablestracklet();
558 } // Fill at the end the charge
561 /////////////////////////////////////////////////////////////
562 // Localise and take the calib gain object for the detector
563 ////////////////////////////////////////////////////////////
564 if (detector != fDetectorPreviousTrack) {
566 //Localise the detector
567 LocalisationDetectorXbins(detector);
570 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
572 AliInfo("Could not get calibDB");
579 fCalROCGain->~AliTRDCalROC();
580 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
582 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
586 // Reset the detectbjobsor
587 fDetectorPreviousTrack = detector;
589 ///////////////////////////////////////
590 // Store the info of the cluster
591 ///////////////////////////////////////
592 Int_t row = cl->GetPadRow();
593 Int_t col = cl->GetPadCol();
594 CheckGoodTrackletV1(cl);
595 Int_t group[2] = {0,0};
596 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
597 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
598 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
602 } // while on clusters
604 ///////////////////////////
605 // Fill the last plane
606 //////////////////////////
607 if( index0 != index1 ){
613 // drift velocity unables to cut bad tracklets
614 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
618 FillTheInfoOfTheTrackCH(index1-index0);
623 FillTheInfoOfTheTrackPH();
626 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
628 } // if a good tracklet
633 ResetfVariablestracklet();
638 //____________Offline tracking in the AliTRDtracker____________________________
639 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t)
642 // Use AliTRDtrackV1 for the calibration
646 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
647 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
648 AliTRDcluster *cls = 0x0; // shared cluster attached now to the tracklet
649 Bool_t newtr = kTRUE; // new track
652 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
654 AliInfo("Could not get calibDB");
658 ///////////////////////////
659 // loop over the tracklet
660 ///////////////////////////
661 for(Int_t itr = 0; itr < 6; itr++){
663 if(!(tracklet = t->GetTracklet(itr))) continue;
664 if(!tracklet->IsOK()) continue;
666 ResetfVariablestracklet();
668 //////////////////////////////////////////
669 // localisation of the tracklet and dqdl
670 //////////////////////////////////////////
671 Int_t layer = tracklet->GetPlane();
673 while(!(cl = tracklet->GetClusters(ic++))) continue;
674 Int_t detector = cl->GetDetector();
675 if (detector != fDetectorPreviousTrack) {
676 // if not a new track
678 // don't use the rest of this track if in the same plane
679 if (layer == GetLayer(fDetectorPreviousTrack)) {
680 //printf("bad tracklet, same layer for detector %d\n",detector);
684 //Localise the detector bin
685 LocalisationDetectorXbins(detector);
689 fCalROCGain->~AliTRDCalROC();
690 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
692 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
695 fDetectorPreviousTrack = detector;
699 ////////////////////////////
700 // loop over the clusters
701 ////////////////////////////
702 Int_t nbclusters = 0;
703 for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
704 if(!(cl = tracklet->GetClusters(jc))) continue;
707 // Store the info bis of the tracklet
708 Int_t row = cl->GetPadRow();
709 Int_t col = cl->GetPadCol();
710 CheckGoodTrackletV1(cl);
711 Int_t group[2] = {0,0};
712 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
713 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
714 // Add the charge if shared cluster
715 cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
717 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col,cls);
720 ////////////////////////////////////////
721 // Fill the stuffs if a good tracklet
722 ////////////////////////////////////////
725 // drift velocity unables to cut bad tracklets
726 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
728 //printf("pass %d and nbclusters %d\n",pass,nbclusters);
732 FillTheInfoOfTheTrackCH(nbclusters);
737 FillTheInfoOfTheTrackPH();
740 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
742 } // if a good tracklet
748 ///////////////////////////////////////////////////////////////////////////////////
749 // Routine inside the update with AliTRDtrack
750 ///////////////////////////////////////////////////////////////////////////////////
751 //____________Offine tracking in the AliTRDtracker_____________________________
752 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1)
755 // Drift velocity calibration:
756 // Fit the clusters with a straight line
757 // From the slope find the drift velocity
760 //Number of points: if less than 3 return kFALSE
761 Int_t npoints = index1-index0;
762 if(npoints <= 2) return kFALSE;
767 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
768 // parameters of the track
769 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
770 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
771 Double_t tnp = 0.0; // tan angle in the plan xy track
772 if( TMath::Abs(snp) < 1.){
773 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
775 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
776 // tilting pad and cross row
777 Int_t crossrow = 0; // if it crosses a pad row
778 Int_t rowp = -1; // if it crosses a pad row
779 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
780 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
781 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
783 fLinearFitterTracklet->ClearPoints();
784 Double_t dydt = 0.0; // dydt tracklet after straight line fit
785 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
786 Double_t pointError = 0.0; // error after straight line fit
787 Int_t nbli = 0; // number linear fitter points
789 //////////////////////////////
790 // loop over clusters
791 ////////////////////////////
792 for(Int_t k = 0; k < npoints; k++){
794 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
795 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
796 Double_t ycluster = cl->GetY();
797 Int_t time = cl->GetPadTime();
798 Double_t timeis = time/fSf;
799 //Double_t q = cl->GetQ();
800 //See if cross two pad rows
801 Int_t row = cl->GetPadRow();
803 if(row != rowp) crossrow = 1;
805 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
810 //////////////////////////////
812 //////////////////////////////
814 fLinearFitterTracklet->ClearPoints();
818 fLinearFitterTracklet->Eval();
819 fLinearFitterTracklet->GetParameters(pars);
820 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
821 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
823 fLinearFitterTracklet->ClearPoints();
825 /////////////////////////////
827 ////////////////////////////
829 if ( !fDebugStreamer ) {
831 TDirectory *backup = gDirectory;
832 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
833 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
837 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
838 //"snpright="<<snpright<<
839 "npoints="<<npoints<<
841 "detector="<<detector<<
848 "crossrow="<<crossrow<<
849 "errorpar="<<errorpar<<
850 "pointError="<<pointError<<
854 Int_t nbclusters = index1-index0;
855 Int_t layer = GetLayer(fDetectorPreviousTrack);
857 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
858 //"snpright="<<snpright<<
859 "nbclusters="<<nbclusters<<
860 "detector="<<fDetectorPreviousTrack<<
866 ///////////////////////////
868 ///////////////////////////
869 if(npoints < fNumberClusters) return kFALSE;
870 if(npoints > fNumberClustersf) return kFALSE;
871 if(pointError >= 0.1) return kFALSE;
872 if(crossrow == 1) return kFALSE;
874 ////////////////////////////
876 ////////////////////////////
878 //Add to the linear fitter of the detector
879 if( TMath::Abs(snp) < 1.){
880 Double_t x = tnp-dzdx*tnt;
881 if(fLinearFitterDebugOn) {
882 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
883 fEntriesLinearFitter[detector]++;
885 fLinearVdriftFit->Update(detector,x,pars[1]);
891 //____________Offine tracking in the AliTRDtracker_____________________________
892 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
895 // Drift velocity calibration:
896 // Fit the clusters with a straight line
897 // From the slope find the drift velocity
900 ////////////////////////////////////////////////
901 //Number of points: if less than 3 return kFALSE
902 /////////////////////////////////////////////////
903 if(nbclusters <= 2) return kFALSE;
908 // results of the linear fit
909 Double_t dydt = 0.0; // dydt tracklet after straight line fit
910 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
911 Double_t pointError = 0.0; // error after straight line fit
912 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
913 Int_t crossrow = 0; // if it crosses a pad row
914 Int_t rowp = -1; // if it crosses a pad row
915 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
916 fLinearFitterTracklet->ClearPoints();
919 ///////////////////////////////////////////
920 // Take the parameters of the track
921 //////////////////////////////////////////
922 // take now the snp, tnp and tgl from the track
923 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
924 Double_t tnp = 0.0; // dy/dx at the end of the chamber
925 if( TMath::Abs(snp) < 1.){
926 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
928 Double_t tgl = tracklet->GetTgl(); // dz/dl
929 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
931 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
932 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
933 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
934 // at the end with correction due to linear fit
935 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
936 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
939 ////////////////////////////
940 // loop over the clusters
941 ////////////////////////////
943 AliTRDcluster *cl = 0x0;
944 //////////////////////////////
945 // Check no shared clusters
946 //////////////////////////////
947 for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
948 if((cl = tracklet->GetClusters(icc))) crossrow = 1;
950 //////////////////////////////////
952 //////////////////////////////////
953 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
954 if(!(cl = tracklet->GetClusters(ic))) continue;
955 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
957 Double_t ycluster = cl->GetY();
958 Int_t time = cl->GetPadTime();
959 Double_t timeis = time/fSf;
960 //See if cross two pad rows
961 Int_t row = cl->GetPadRow();
962 if(rowp==-1) rowp = row;
963 if(row != rowp) crossrow = 1;
965 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
971 ////////////////////////////////////
972 // Do the straight line fit now
973 ///////////////////////////////////
975 fLinearFitterTracklet->ClearPoints();
979 fLinearFitterTracklet->Eval();
980 fLinearFitterTracklet->GetParameters(pars);
981 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
982 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
984 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
985 fLinearFitterTracklet->ClearPoints();
987 ////////////////////////////////
989 ///////////////////////////////
993 if ( !fDebugStreamer ) {
995 TDirectory *backup = gDirectory;
996 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
997 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1001 Int_t layer = GetLayer(fDetectorPreviousTrack);
1003 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1004 //"snpright="<<snpright<<
1006 "nbclusters="<<nbclusters<<
1007 "detector="<<fDetectorPreviousTrack<<
1015 "crossrow="<<crossrow<<
1016 "errorpar="<<errorpar<<
1017 "pointError="<<pointError<<
1022 /////////////////////////
1024 ////////////////////////
1026 if(nbclusters < fNumberClusters) return kFALSE;
1027 if(nbclusters > fNumberClustersf) return kFALSE;
1028 if(pointError >= 0.3) return kFALSE;
1029 if(crossrow == 1) return kFALSE;
1031 ///////////////////////
1033 //////////////////////
1035 if(fLinearFitterOn){
1036 //Add to the linear fitter of the detector
1037 if( TMath::Abs(snp) < 1.){
1038 Double_t x = tnp-dzdx*tnt;
1039 if(fLinearFitterDebugOn) {
1040 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1041 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1043 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1049 //____________Offine tracking in the AliTRDtracker_____________________________
1050 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(const AliTRDtrack *t, Int_t index0, Int_t index1)
1053 // PRF width calibration
1054 // Assume a Gaussian shape: determinate the position of the three pad clusters
1055 // Fit with a straight line
1056 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1057 // Fill the PRF as function of angle of the track
1062 //////////////////////////
1064 /////////////////////////
1065 Int_t npoints = index1-index0; // number of total points
1066 Int_t nb3pc = 0; // number of three pads clusters used for fit
1067 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1068 // To see the difference due to the fit
1069 Double_t *padPositions;
1070 padPositions = new Double_t[npoints];
1071 for(Int_t k = 0; k < npoints; k++){
1072 padPositions[k] = 0.0;
1074 // Take the tgl and snp with the track t now
1075 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1076 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1077 Float_t dzdx = 0.0; // dzdx
1079 if(TMath::Abs(snp) < 1.0){
1080 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1081 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1084 fLinearFitterTracklet->ClearPoints();
1086 ///////////////////////////
1087 // calcul the tnp group
1088 ///////////////////////////
1089 Bool_t echec = kFALSE;
1090 Double_t shift = 0.0;
1091 //Calculate the shift in x coresponding to this tnp
1092 if(fNgroupprf != 0.0){
1093 shift = -3.0*(fNgroupprf-1)-1.5;
1094 Double_t limithigh = -0.2*(fNgroupprf-1);
1095 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1097 while(tnp > limithigh){
1104 delete [] padPositions;
1108 //////////////////////
1110 /////////////////////
1111 for(Int_t k = 0; k < npoints; k++){
1113 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1115 Short_t *signals = cl->GetSignals();
1116 Double_t time = cl->GetPadTime();
1117 //Calculate x if possible
1118 Float_t xcenter = 0.0;
1119 Bool_t echec1 = kTRUE;
1120 if((time<=7) || (time>=21)) continue;
1121 // Center 3 balanced: position with the center of the pad
1122 if ((((Float_t) signals[3]) > 0.0) &&
1123 (((Float_t) signals[2]) > 0.0) &&
1124 (((Float_t) signals[4]) > 0.0)) {
1126 // Security if the denomiateur is 0
1127 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1128 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1129 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1130 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1131 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1137 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1139 //if no echec: calculate with the position of the pad
1140 // Position of the cluster
1141 Double_t padPosition = xcenter + cl->GetPadCol();
1142 padPositions[k] = padPosition;
1144 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1148 /////////////////////////////
1150 ////////////////////////////
1152 delete [] padPositions;
1153 fLinearFitterTracklet->ClearPoints();
1156 fLinearFitterTracklet->Eval();
1158 fLinearFitterTracklet->GetParameters(line);
1159 Float_t pointError = -1.0;
1160 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1161 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1163 fLinearFitterTracklet->ClearPoints();
1165 /////////////////////////////////////////////////////
1166 // Now fill the PRF: second loop over clusters
1167 /////////////////////////////////////////////////////
1168 for(Int_t k = 0; k < npoints; k++){
1170 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1171 Short_t *signals = cl->GetSignals(); // signal
1172 Double_t time = cl->GetPadTime(); // time bin
1173 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1174 Float_t padPos = cl->GetPadCol(); // middle pad
1175 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1176 Float_t ycenter = 0.0; // relative center charge
1177 Float_t ymin = 0.0; // relative left charge
1178 Float_t ymax = 0.0; // relative right charge
1180 //Requiere simply two pads clusters at least
1181 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1182 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1183 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1184 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1185 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1186 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1190 Int_t rowcl = cl->GetPadRow(); // row of cluster
1191 Int_t colcl = cl->GetPadCol(); // col of cluster
1192 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1193 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1194 Float_t xcl = cl->GetY(); // y cluster
1195 Float_t qcl = cl->GetQ(); // charge cluster
1196 Int_t layer = GetLayer(detector); // layer
1197 Int_t stack = GetStack(detector); // stack
1198 Double_t xdiff = dpad; // reconstructed position constant
1199 Double_t x = dpad; // reconstructed position moved
1200 Float_t ep = pointError; // error of fit
1201 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1202 Float_t signal3 = (Float_t)signals[3]; // signal
1203 Float_t signal2 = (Float_t)signals[2]; // signal
1204 Float_t signal4 = (Float_t)signals[4]; // signal
1205 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1207 //////////////////////////////
1209 /////////////////////////////
1210 if(fDebugLevel > 0){
1211 if ( !fDebugStreamer ) {
1213 TDirectory *backup = gDirectory;
1214 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1215 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1221 Float_t y = ycenter;
1222 (* fDebugStreamer) << "HandlePRFtracklet"<<
1223 "caligroup="<<caligroup<<
1224 "detector="<<detector<<
1227 "npoints="<<npoints<<
1236 "padPosition="<<padPositions[k]<<
1237 "padPosTracklet="<<padPosTracklet<<
1242 "signal1="<<signal1<<
1243 "signal2="<<signal2<<
1244 "signal3="<<signal3<<
1245 "signal4="<<signal4<<
1246 "signal5="<<signal5<<
1252 (* fDebugStreamer) << "HandlePRFtracklet"<<
1253 "caligroup="<<caligroup<<
1254 "detector="<<detector<<
1257 "npoints="<<npoints<<
1266 "padPosition="<<padPositions[k]<<
1267 "padPosTracklet="<<padPosTracklet<<
1272 "signal1="<<signal1<<
1273 "signal2="<<signal2<<
1274 "signal3="<<signal3<<
1275 "signal4="<<signal4<<
1276 "signal5="<<signal5<<
1282 (* fDebugStreamer) << "HandlePRFtracklet"<<
1283 "caligroup="<<caligroup<<
1284 "detector="<<detector<<
1287 "npoints="<<npoints<<
1296 "padPosition="<<padPositions[k]<<
1297 "padPosTracklet="<<padPosTracklet<<
1302 "signal1="<<signal1<<
1303 "signal2="<<signal2<<
1304 "signal3="<<signal3<<
1305 "signal4="<<signal4<<
1306 "signal5="<<signal5<<
1312 ////////////////////////////
1314 ///////////////////////////
1315 if(npoints < fNumberClusters) continue;
1316 if(npoints > fNumberClustersf) continue;
1317 if(nb3pc <= 5) continue;
1318 if((time >= 21) || (time < 7)) continue;
1319 if(TMath::Abs(snp) >= 1.0) continue;
1320 if(TMath::Abs(qcl) < 80) continue;
1322 ////////////////////////////
1324 ///////////////////////////
1326 if(TMath::Abs(dpad) < 1.5) {
1327 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1328 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1330 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1331 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1332 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1334 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1335 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1336 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1340 if(TMath::Abs(dpad) < 1.5) {
1341 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1342 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1344 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1345 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1346 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1348 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1349 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1350 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1354 delete [] padPositions;
1358 //____________Offine tracking in the AliTRDtracker_____________________________
1359 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1362 // PRF width calibration
1363 // Assume a Gaussian shape: determinate the position of the three pad clusters
1364 // Fit with a straight line
1365 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1366 // Fill the PRF as function of angle of the track
1370 //printf("begin\n");
1371 ///////////////////////////////////////////
1372 // Take the parameters of the track
1373 //////////////////////////////////////////
1374 // take now the snp, tnp and tgl from the track
1375 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1376 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1377 if( TMath::Abs(snp) < 1.){
1378 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1380 Double_t tgl = tracklet->GetTgl(); // dz/dl
1381 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1383 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1384 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1385 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1386 // at the end with correction due to linear fit
1387 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1388 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1390 ///////////////////////////////
1391 // Calculate tnp group shift
1392 ///////////////////////////////
1393 Bool_t echec = kFALSE;
1394 Double_t shift = 0.0;
1395 //Calculate the shift in x coresponding to this tnp
1396 if(fNgroupprf != 0.0){
1397 shift = -3.0*(fNgroupprf-1)-1.5;
1398 Double_t limithigh = -0.2*(fNgroupprf-1);
1399 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1401 while(tnp > limithigh){
1407 // do nothing if out of tnp range
1408 //printf("echec %d\n",(Int_t)echec);
1409 if(echec) return kFALSE;
1411 ///////////////////////
1413 //////////////////////
1415 Int_t nb3pc = 0; // number of three pads clusters used for fit
1416 // to see the difference between the fit and the 3 pad clusters position
1417 Double_t padPositions[AliTRDseedV1::kNtb];
1418 memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1419 fLinearFitterTracklet->ClearPoints();
1421 //printf("loop clusters \n");
1422 ////////////////////////////
1423 // loop over the clusters
1424 ////////////////////////////
1425 AliTRDcluster *cl = 0x0;
1426 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1427 // reject shared clusters on pad row
1428 if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1429 if((cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1431 if(!(cl = tracklet->GetClusters(ic))) continue;
1432 Double_t time = cl->GetPadTime();
1433 if((time<=7) || (time>=21)) continue;
1434 Short_t *signals = cl->GetSignals();
1435 Float_t xcenter = 0.0;
1436 Bool_t echec1 = kTRUE;
1438 /////////////////////////////////////////////////////////////
1439 // Center 3 balanced: position with the center of the pad
1440 /////////////////////////////////////////////////////////////
1441 if ((((Float_t) signals[3]) > 0.0) &&
1442 (((Float_t) signals[2]) > 0.0) &&
1443 (((Float_t) signals[4]) > 0.0)) {
1445 // Security if the denomiateur is 0
1446 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1447 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1448 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1449 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1450 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1456 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1457 if(echec1) continue;
1459 ////////////////////////////////////////////////////////
1460 //if no echec1: calculate with the position of the pad
1461 // Position of the cluster
1462 // fill the linear fitter
1463 ///////////////////////////////////////////////////////
1464 Double_t padPosition = xcenter + cl->GetPadCol();
1465 padPositions[ic] = padPosition;
1467 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1472 //printf("Fin loop clusters \n");
1473 //////////////////////////////
1474 // fit with a straight line
1475 /////////////////////////////
1477 fLinearFitterTracklet->ClearPoints();
1480 fLinearFitterTracklet->Eval();
1482 fLinearFitterTracklet->GetParameters(line);
1483 Float_t pointError = -1.0;
1484 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1485 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1487 fLinearFitterTracklet->ClearPoints();
1489 //printf("PRF second loop \n");
1490 ////////////////////////////////////////////////
1491 // Fill the PRF: Second loop over clusters
1492 //////////////////////////////////////////////
1493 for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1494 // reject shared clusters on pad row
1495 if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1497 if(!(cl = tracklet->GetClusters(ic))) continue;
1499 Short_t *signals = cl->GetSignals(); // signal
1500 Double_t time = cl->GetPadTime(); // time bin
1501 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1502 Float_t padPos = cl->GetPadCol(); // middle pad
1503 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1504 Float_t ycenter = 0.0; // relative center charge
1505 Float_t ymin = 0.0; // relative left charge
1506 Float_t ymax = 0.0; // relative right charge
1508 ////////////////////////////////////////////////////////////////
1509 // Calculate ycenter, ymin and ymax even for two pad clusters
1510 ////////////////////////////////////////////////////////////////
1511 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1512 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1513 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1514 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1515 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1516 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1519 /////////////////////////
1520 // Calibration group
1521 ////////////////////////
1522 Int_t rowcl = cl->GetPadRow(); // row of cluster
1523 Int_t colcl = cl->GetPadCol(); // col of cluster
1524 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1525 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1526 Float_t xcl = cl->GetY(); // y cluster
1527 Float_t qcl = cl->GetQ(); // charge cluster
1528 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1529 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1530 Double_t xdiff = dpad; // reconstructed position constant
1531 Double_t x = dpad; // reconstructed position moved
1532 Float_t ep = pointError; // error of fit
1533 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1534 Float_t signal3 = (Float_t)signals[3]; // signal
1535 Float_t signal2 = (Float_t)signals[2]; // signal
1536 Float_t signal4 = (Float_t)signals[4]; // signal
1537 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1541 /////////////////////
1543 ////////////////////
1545 if(fDebugLevel > 0){
1546 if ( !fDebugStreamer ) {
1548 TDirectory *backup = gDirectory;
1549 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1550 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1555 Float_t y = ycenter;
1556 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1557 "caligroup="<<caligroup<<
1558 "detector="<<fDetectorPreviousTrack<<
1561 "npoints="<<nbclusters<<
1570 "padPosition="<<padPositions[ic]<<
1571 "padPosTracklet="<<padPosTracklet<<
1576 "signal1="<<signal1<<
1577 "signal2="<<signal2<<
1578 "signal3="<<signal3<<
1579 "signal4="<<signal4<<
1580 "signal5="<<signal5<<
1586 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1587 "caligroup="<<caligroup<<
1588 "detector="<<fDetectorPreviousTrack<<
1591 "npoints="<<nbclusters<<
1600 "padPosition="<<padPositions[ic]<<
1601 "padPosTracklet="<<padPosTracklet<<
1606 "signal1="<<signal1<<
1607 "signal2="<<signal2<<
1608 "signal3="<<signal3<<
1609 "signal4="<<signal4<<
1610 "signal5="<<signal5<<
1616 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1617 "caligroup="<<caligroup<<
1618 "detector="<<fDetectorPreviousTrack<<
1621 "npoints="<<nbclusters<<
1630 "padPosition="<<padPositions[ic]<<
1631 "padPosTracklet="<<padPosTracklet<<
1636 "signal1="<<signal1<<
1637 "signal2="<<signal2<<
1638 "signal3="<<signal3<<
1639 "signal4="<<signal4<<
1640 "signal5="<<signal5<<
1646 /////////////////////
1648 /////////////////////
1649 if(nbclusters < fNumberClusters) continue;
1650 if(nbclusters > fNumberClustersf) continue;
1651 if(nb3pc <= 5) continue;
1652 if((time >= 21) || (time < 7)) continue;
1653 if(TMath::Abs(qcl) < 80) continue;
1654 if( TMath::Abs(snp) > 1.) continue;
1657 ////////////////////////
1659 ///////////////////////
1661 if(TMath::Abs(dpad) < 1.5) {
1662 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1663 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1664 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1666 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1667 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1668 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1670 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1671 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1672 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1677 if(TMath::Abs(dpad) < 1.5) {
1678 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1679 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1681 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1682 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1683 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1685 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1686 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1687 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1690 } // second loop over clusters
1695 ///////////////////////////////////////////////////////////////////////////////////////
1696 // Pad row col stuff: see if masked or not
1697 ///////////////////////////////////////////////////////////////////////////////////////
1698 //_____________________________________________________________________________
1699 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1702 // See if we are not near a masked pad
1705 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1709 //_____________________________________________________________________________
1710 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1713 // See if we are not near a masked pad
1716 if (!IsPadOn(detector, col, row)) {
1717 fGoodTracklet = kFALSE;
1721 if (!IsPadOn(detector, col-1, row)) {
1722 fGoodTracklet = kFALSE;
1727 if (!IsPadOn(detector, col+1, row)) {
1728 fGoodTracklet = kFALSE;
1733 //_____________________________________________________________________________
1734 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1737 // Look in the choosen database if the pad is On.
1738 // If no the track will be "not good"
1741 // Get the parameter object
1742 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1744 AliInfo("Could not get calibDB");
1748 if (!cal->IsChamberInstalled(detector) ||
1749 cal->IsChamberMasked(detector) ||
1750 cal->IsPadMasked(detector,col,row)) {
1758 ///////////////////////////////////////////////////////////////////////////////////////
1759 // Calibration groups: calculate the number of groups, localise...
1760 ////////////////////////////////////////////////////////////////////////////////////////
1761 //_____________________________________________________________________________
1762 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1765 // Calculate the calibration group number for i
1768 // Row of the cluster and position in the pad groups
1770 if (fCalibraMode->GetNnZ(i) != 0) {
1771 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1775 // Col of the cluster and position in the pad groups
1777 if (fCalibraMode->GetNnRphi(i) != 0) {
1778 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1781 return posc*fCalibraMode->GetNfragZ(i)+posr;
1784 //____________________________________________________________________________________
1785 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1788 // Calculate the total number of calibration groups
1794 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1796 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1801 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1803 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1808 fCalibraMode->ModePadCalibration(2,i);
1809 fCalibraMode->ModePadFragmentation(0,2,0,i);
1810 fCalibraMode->SetDetChamb2(i);
1811 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1812 fCalibraMode->ModePadCalibration(0,i);
1813 fCalibraMode->ModePadFragmentation(0,0,0,i);
1814 fCalibraMode->SetDetChamb0(i);
1815 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1816 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1820 //_____________________________________________________________________________
1821 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1824 // Set the mode of calibration group in the z direction for the parameter i
1829 fCalibraMode->SetNz(i, Nz);
1832 AliInfo("You have to choose between 0 and 4");
1836 //_____________________________________________________________________________
1837 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1840 // Set the mode of calibration group in the rphi direction for the parameter i
1845 fCalibraMode->SetNrphi(i ,Nrphi);
1848 AliInfo("You have to choose between 0 and 6");
1853 //_____________________________________________________________________________
1854 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1857 // Set the mode of calibration group all together
1859 if(fVector2d == kTRUE) {
1860 AliInfo("Can not work with the vector method");
1863 fCalibraMode->SetAllTogether(i);
1867 //_____________________________________________________________________________
1868 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1871 // Set the mode of calibration group per supermodule
1873 if(fVector2d == kTRUE) {
1874 AliInfo("Can not work with the vector method");
1877 fCalibraMode->SetPerSuperModule(i);
1881 //____________Set the pad calibration variables for the detector_______________
1882 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1885 // For the detector calcul the first Xbins and set the number of row
1886 // and col pads per calibration groups, the number of calibration
1887 // groups in the detector.
1890 // first Xbins of the detector
1892 fCalibraMode->CalculXBins(detector,0);
1895 fCalibraMode->CalculXBins(detector,1);
1898 fCalibraMode->CalculXBins(detector,2);
1901 // fragmentation of idect
1902 for (Int_t i = 0; i < 3; i++) {
1903 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1904 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1905 , (Int_t) GetStack(detector)
1906 , (Int_t) GetSector(detector),i);
1912 //_____________________________________________________________________________
1913 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1916 // Should be between 0 and 6
1919 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1920 AliInfo("The number of groups must be between 0 and 6!");
1923 fNgroupprf = numberGroupsPRF;
1927 ///////////////////////////////////////////////////////////////////////////////////////////
1928 // Per tracklet: store or reset the info, fill the histos with the info
1929 //////////////////////////////////////////////////////////////////////////////////////////
1930 //_____________________________________________________________________________
1931 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)
1934 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1935 // Correct from the gain correction before
1936 // cls is shared cluster if any
1939 //printf("StoreInfoCHPHtrack\n");
1941 // time bin of the cluster not corrected
1942 Int_t time = cl->GetPadTime();
1943 Float_t charge = TMath::Abs(cl->GetQ());
1945 charge += TMath::Abs(cls->GetQ());
1946 //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1949 //printf("Store::time %d, amplitude %f\n",time,dqdl);
1951 //Correct for the gain coefficient used in the database for reconstruction
1952 Float_t correctthegain = 1.0;
1953 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1954 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1955 Float_t correction = 1.0;
1956 Float_t normalisation = 6.67;
1957 // we divide with gain in AliTRDclusterizer::Transform...
1958 if( correctthegain > 0 ) normalisation /= correctthegain;
1961 // take dd/dl corrected from the angle of the track
1962 correction = dqdl / (normalisation);
1965 // Fill the fAmpTotal with the charge
1967 if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1968 //printf("Store::group %d, amplitude %f\n",group[0],correction);
1969 fAmpTotal[(Int_t) group[0]] += correction;
1973 // Fill the fPHPlace and value
1975 fPHPlace[time] = group[1];
1976 fPHValue[time] = charge;
1980 //____________Offine tracking in the AliTRDtracker_____________________________
1981 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1984 // Reset values per tracklet
1987 //Reset good tracklet
1988 fGoodTracklet = kTRUE;
1990 // Reset the fPHValue
1992 //Reset the fPHValue and fPHPlace
1993 for (Int_t k = 0; k < fTimeMax; k++) {
1999 // Reset the fAmpTotal where we put value
2001 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2006 //____________Offine tracking in the AliTRDtracker_____________________________
2007 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
2010 // For the offline tracking or mcm tracklets
2011 // This function will be called in the functions UpdateHistogram...
2012 // to fill the info of a track for the relativ gain calibration
2015 Int_t nb = 0; // Nombre de zones traversees
2016 Int_t fd = -1; // Premiere zone non nulle
2017 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
2019 //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2021 if(nbclusters < fNumberClusters) return;
2022 if(nbclusters > fNumberClustersf) return;
2025 // Normalize with the number of clusters
2026 Double_t normalizeCst = fRelativeScale;
2027 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
2029 //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
2031 // See if the track goes through different zones
2032 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2033 //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
2034 if (fAmpTotal[k] > 0.0) {
2035 totalcharge += fAmpTotal[k];
2043 //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
2049 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2051 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2052 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2055 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2059 if ((fAmpTotal[fd] > 0.0) &&
2060 (fAmpTotal[fd+1] > 0.0)) {
2061 // One of the two very big
2062 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2064 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2065 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2068 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2071 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2073 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2075 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2076 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2079 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2082 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2085 if (fCalibraMode->GetNfragZ(0) > 1) {
2086 if (fAmpTotal[fd] > 0.0) {
2087 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2088 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2089 // One of the two very big
2090 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2092 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2093 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2096 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2099 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2101 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2103 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2104 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2107 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2109 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2120 //____________Offine tracking in the AliTRDtracker_____________________________
2121 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2124 // For the offline tracking or mcm tracklets
2125 // This function will be called in the functions UpdateHistogram...
2126 // to fill the info of a track for the drift velocity calibration
2129 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2130 Int_t fd1 = -1; // Premiere zone non nulle
2131 Int_t fd2 = -1; // Deuxieme zone non nulle
2132 Int_t k1 = -1; // Debut de la premiere zone
2133 Int_t k2 = -1; // Debut de la seconde zone
2134 Int_t nbclusters = 0; // number of clusters
2138 // See if the track goes through different zones
2139 for (Int_t k = 0; k < fTimeMax; k++) {
2140 if (fPHValue[k] > 0.0) {
2146 if (fPHPlace[k] != fd1) {
2152 if (fPHPlace[k] != fd2) {
2159 // See if noise before and after
2160 if(fMaxCluster > 0) {
2161 if(fPHValue[0] > fMaxCluster) return;
2162 if(fTimeMax > fNbMaxCluster) {
2163 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2164 if(fPHValue[k] > fMaxCluster) return;
2169 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2171 if(nbclusters < fNumberClusters) return;
2172 if(nbclusters > fNumberClustersf) return;
2178 for (Int_t i = 0; i < fTimeMax; i++) {
2180 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2182 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2184 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2187 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2189 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2195 if ((fd1 == fd2+1) ||
2197 // One of the two fast all the think
2198 if (k2 > (k1+fDifference)) {
2199 //we choose to fill the fd1 with all the values
2201 for (Int_t i = 0; i < fTimeMax; i++) {
2203 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2205 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2209 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2211 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2216 if ((k2+fDifference) < fTimeMax) {
2217 //we choose to fill the fd2 with all the values
2219 for (Int_t i = 0; i < fTimeMax; i++) {
2221 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2223 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2227 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2229 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2235 // Two zones voisines sinon rien!
2236 if (fCalibraMode->GetNfragZ(1) > 1) {
2238 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2239 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2240 // One of the two fast all the think
2241 if (k2 > (k1+fDifference)) {
2242 //we choose to fill the fd1 with all the values
2244 for (Int_t i = 0; i < fTimeMax; i++) {
2246 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2248 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2252 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2254 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2259 if ((k2+fDifference) < fTimeMax) {
2260 //we choose to fill the fd2 with all the values
2262 for (Int_t i = 0; i < fTimeMax; i++) {
2264 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2266 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2270 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2272 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2279 // Two zones voisines sinon rien!
2281 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2282 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2283 // One of the two fast all the think
2284 if (k2 > (k1 + fDifference)) {
2285 //we choose to fill the fd1 with all the values
2287 for (Int_t i = 0; i < fTimeMax; i++) {
2289 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2291 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2295 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2297 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2302 if ((k2+fDifference) < fTimeMax) {
2303 //we choose to fill the fd2 with all the values
2305 for (Int_t i = 0; i < fTimeMax; i++) {
2307 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2309 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2313 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2315 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2327 //////////////////////////////////////////////////////////////////////////////////////////
2328 // DAQ process functions
2329 /////////////////////////////////////////////////////////////////////////////////////////
2330 //_____________________________________________________________________
2331 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2334 // Event Processing loop - AliTRDrawStreamBase
2335 // TestBeam 2007 version
2336 // 0 timebin problem
2339 // Same algorithm as TestBeam but different reader
2342 rawStream->SetSharedPadReadout(kFALSE);
2344 Int_t withInput = 1;
2346 Double_t phvalue[16][144][36];
2347 for(Int_t k = 0; k < 36; k++){
2348 for(Int_t j = 0; j < 16; j++){
2349 for(Int_t c = 0; c < 144; c++){
2350 phvalue[j][c][k] = 0.0;
2355 fDetectorPreviousTrack = -1;
2358 Int_t nbtimebin = 0;
2359 Int_t baseline = 10;
2360 //printf("------------Detector\n");
2366 while (rawStream->Next()) {
2368 Int_t idetector = rawStream->GetDet(); // current detector
2369 Int_t imcm = rawStream->GetMCM(); // current MCM
2370 Int_t irob = rawStream->GetROB(); // current ROB
2372 //printf("Detector %d\n",idetector);
2374 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2377 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2381 for(Int_t k = 0; k < 36; k++){
2382 for(Int_t j = 0; j < 16; j++){
2383 for(Int_t c = 0; c < 144; c++){
2384 phvalue[j][c][k] = 0.0;
2390 fDetectorPreviousTrack = idetector;
2391 fMCMPrevious = imcm;
2392 fROBPrevious = irob;
2394 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2395 if(nbtimebin == 0) return 0;
2396 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2397 fTimeMax = nbtimebin;
2399 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2400 fNumberClustersf = fTimeMax;
2401 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2404 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2405 Int_t col = rawStream->GetCol();
2406 Int_t row = rawStream->GetRow();
2409 // printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2412 for(Int_t itime = 0; itime < nbtimebin; itime++){
2413 phvalue[row][col][itime] = signal[itime]-baseline;
2417 // fill the last one
2418 if(fDetectorPreviousTrack != -1){
2421 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2424 for(Int_t k = 0; k < 36; k++){
2425 for(Int_t j = 0; j < 16; j++){
2426 for(Int_t c = 0; c < 144; c++){
2427 phvalue[j][c][k] = 0.0;
2436 while (rawStream->Next()) { //iddetecte
2438 Int_t idetector = rawStream->GetDet(); // current detector
2439 Int_t imcm = rawStream->GetMCM(); // current MCM
2440 Int_t irob = rawStream->GetROB(); // current ROB
2442 //printf("Detector %d\n",idetector);
2444 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2447 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2450 for(Int_t k = 0; k < 36; k++){
2451 for(Int_t j = 0; j < 16; j++){
2452 for(Int_t c = 0; c < 144; c++){
2453 phvalue[j][c][k] = 0.0;
2459 fDetectorPreviousTrack = idetector;
2460 fMCMPrevious = imcm;
2461 fROBPrevious = irob;
2463 //baseline = rawStream->GetCommonAdditive(); // common baseline
2465 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2466 fNumberClustersf = fTimeMax;
2467 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2468 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2469 Int_t col = rawStream->GetCol();
2470 Int_t row = rawStream->GetRow();
2473 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2475 for(Int_t itime = 0; itime < fTimeMax; itime++){
2476 phvalue[row][col][itime] = signal[itime]-baseline;
2477 /*if(phvalue[row][col][itime] >= 20) {
2478 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2488 // fill the last one
2489 if(fDetectorPreviousTrack != -1){
2492 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2495 for(Int_t k = 0; k < 36; k++){
2496 for(Int_t j = 0; j < 16; j++){
2497 for(Int_t c = 0; c < 144; c++){
2498 phvalue[j][c][k] = 0.0;
2508 //_____________________________________________________________________
2509 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2512 // Event processing loop - AliRawReader
2513 // Testbeam 2007 version
2516 AliTRDrawStreamBase rawStream(rawReader);
2518 rawReader->Select("TRD");
2520 return ProcessEventDAQ(&rawStream, nocheck);
2523 //_________________________________________________________________________
2524 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2526 const eventHeaderStruct *event,
2529 const eventHeaderStruct* /*event*/,
2536 // process date event
2537 // Testbeam 2007 version
2540 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2541 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2545 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2550 //_____________________________________________________________________
2551 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ2(AliRawReader *rawReader)
2554 // Event Processing loop - AliTRDrawFastStream
2556 // 0 timebin problem
2559 // Same algorithm as TestBeam but different reader
2562 // AliTRDrawFastStream *rawStream = AliTRDrawFastStream::GetRawStream(rawReader);
2563 AliTRDrawFastStream *rawStream = new AliTRDrawFastStream(rawReader);
2564 rawStream->SetNoErrorWarning();
2565 rawStream->SetSharedPadReadout(kFALSE);
2567 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2568 digitsManager->CreateArrays();
2570 Int_t withInput = 1;
2572 Double_t phvalue[16][144][36];
2573 for(Int_t k = 0; k < 36; k++){
2574 for(Int_t j = 0; j < 16; j++){
2575 for(Int_t c = 0; c < 144; c++){
2576 phvalue[j][c][k] = 0.0;
2581 fDetectorPreviousTrack = -1;
2585 Int_t nbtimebin = 0;
2586 Int_t baseline = 10;
2592 while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2594 if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2595 // printf("there is ADC data on this chamber!\n");
2597 AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2598 if (digits->HasData()) { //array
2600 AliTRDSignalIndex *indexes = digitsManager->GetIndexes(det);
2601 if (indexes->IsAllocated() == kFALSE) {
2602 AliError("Indexes do not exist!");
2607 indexes->ResetCounters();
2609 while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2610 //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2611 //while (rawStream->Next()) {
2613 Int_t idetector = det; // current detector
2614 //Int_t imcm = rawStream->GetMCM(); // current MCM
2615 //Int_t irob = rawStream->GetROB(); // current ROB
2618 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2620 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2623 for(Int_t k = 0; k < 36; k++){
2624 for(Int_t j = 0; j < 16; j++){
2625 for(Int_t c = 0; c < 144; c++){
2626 phvalue[j][c][k] = 0.0;
2632 fDetectorPreviousTrack = idetector;
2633 //fMCMPrevious = imcm;
2634 //fROBPrevious = irob;
2636 // nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2637 AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2638 nbtimebin = digitParam->GetNTimeBins(); // number of time bins read from data
2639 baseline = digitParam->GetADCbaseline(); // baseline
2641 if(nbtimebin == 0) return 0;
2642 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2643 fTimeMax = nbtimebin;
2645 fNumberClustersf = fTimeMax;
2646 fNumberClusters = (Int_t)(fNumberClustersProcent*fTimeMax);
2649 for(Int_t itime = 0; itime < nbtimebin; itime++) {
2650 // phvalue[row][col][itime] = signal[itime]-baseline;
2651 phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2652 /*if(phvalue[iRow][iCol][itime] >= 20) {
2653 printf("----------> phvalue[%d][%d][%d] %d baseline %d \n",
2657 (Short_t)(digits->GetData(iRow,iCol,itime)),
2664 // fill the last one
2665 if(fDetectorPreviousTrack != -1){
2668 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2669 // printf("\n ---> withinput %d\n\n",withInput);
2671 for(Int_t k = 0; k < 36; k++){
2672 for(Int_t j = 0; j < 16; j++){
2673 for(Int_t c = 0; c < 144; c++){
2674 phvalue[j][c][k] = 0.0;
2682 digitsManager->ClearArrays(det);
2684 delete digitsManager;
2689 //_____________________________________________________________________
2690 //////////////////////////////////////////////////////////////////////////////
2691 // Routine inside the DAQ process
2692 /////////////////////////////////////////////////////////////////////////////
2693 //_______________________________________________________________________
2694 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2697 // Look for the maximum by collapsing over the time
2698 // Sum over four pad col and two pad row
2704 Int_t idect = fDetectorPreviousTrack;
2705 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2707 for(Int_t tb = 0; tb < 36; tb++){
2711 //fGoodTracklet = kTRUE;
2712 //fDetectorPreviousTrack = 0;
2715 ///////////////////////////
2717 /////////////////////////
2721 Double_t integralMax = -1;
2723 for (Int_t ir = 1; ir <= 15; ir++)
2725 for (Int_t ic = 2; ic <= 142; ic++)
2727 Double_t integral = 0;
2728 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2730 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2732 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2733 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2736 for(Int_t tb = 0; tb< fTimeMax; tb++){
2737 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2742 if (integralMax < integral)
2746 integralMax = integral;
2748 } // check max integral
2752 // printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2753 //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2758 if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2762 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2763 //if(!fGoodTracklet) used = 1;;
2765 // /////////////////////////////////////////////////////
2766 // sum ober 2 row and 4 pad cols for each time bins
2767 // ////////////////////////////////////////////////////
2771 for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2773 for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2775 if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2776 imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2778 for(Int_t it = 0; it < fTimeMax; it++){
2779 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2786 Double_t sumcharge = 0.0;
2787 for(Int_t it = 0; it < fTimeMax; it++){
2788 sumcharge += sum[it];
2789 if(sum[it] > fThresholdClustersDAQ) nbcl++;
2793 /////////////////////////////////////////////////////////
2795 ////////////////////////////////////////////////////////
2796 if(fDebugLevel > 0){
2797 if ( !fDebugStreamer ) {
2799 TDirectory *backup = gDirectory;
2800 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2801 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2804 Double_t amph0 = sum[0];
2805 Double_t amphlast = sum[fTimeMax-1];
2806 Double_t rms = TMath::RMS(fTimeMax,sum);
2807 Int_t goodtracklet = (Int_t) fGoodTracklet;
2808 for(Int_t it = 0; it < fTimeMax; it++){
2809 Double_t clustera = sum[it];
2811 (* fDebugStreamer) << "FillDAQa"<<
2812 "ampTotal="<<sumcharge<<
2815 "detector="<<idect<<
2817 "amphlast="<<amphlast<<
2818 "goodtracklet="<<goodtracklet<<
2819 "clustera="<<clustera<<
2827 ////////////////////////////////////////////////////////
2829 ///////////////////////////////////////////////////////
2830 //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2831 if(sum[0] > 100.0) used = 1;
2832 if(nbcl < fNumberClusters) used = 1;
2833 if(nbcl > fNumberClustersf) used = 1;
2835 //if(fDetectorPreviousTrack == 15){
2836 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2838 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2840 for(Int_t it = 0; it < fTimeMax; it++){
2841 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2843 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2845 //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2847 // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2852 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2854 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2861 //____________Online trackling in AliTRDtrigger________________________________
2862 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2866 // Fill a simple average pulse height
2870 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2876 //____________Write_____________________________________________________
2877 //_____________________________________________________________________
2878 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2881 // Write infos to file
2885 if ( fDebugStreamer ) {
2886 delete fDebugStreamer;
2887 fDebugStreamer = 0x0;
2890 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2895 ,fNumberUsedPh[1]));
2897 TDirectory *backup = gDirectory;
2903 option = "recreate";
2905 TFile f(filename,option.Data());
2907 TStopwatch stopwatch;
2910 f.WriteTObject(fCalibraVector);
2915 f.WriteTObject(fCH2d);
2920 f.WriteTObject(fPH2d);
2925 f.WriteTObject(fPRF2d);
2928 if(fLinearFitterOn){
2929 if(fLinearFitterDebugOn) AnalyseLinearFitter();
2930 f.WriteTObject(fLinearVdriftFit);
2935 if ( backup ) backup->cd();
2937 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2938 ,stopwatch.RealTime(),stopwatch.CpuTime()));
2940 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2942 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2943 //___________________________________________probe the histos__________________________________________________
2944 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2947 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2948 // debug mode with 2 for TH2I and 3 for TProfile2D
2949 // It gives a pointer to a Double_t[7] with the info following...
2950 // [0] : number of calibration groups with entries
2951 // [1] : minimal number of entries found
2952 // [2] : calibration group number of the min
2953 // [3] : maximal number of entries found
2954 // [4] : calibration group number of the max
2955 // [5] : mean number of entries found
2956 // [6] : mean relative error
2959 Double_t *info = new Double_t[7];
2961 // Number of Xbins (detectors or groups of pads)
2962 Int_t nbins = h->GetNbinsY(); //number of calibration groups
2963 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
2966 Double_t nbwe = 0; //number of calibration groups with entries
2967 Double_t minentries = 0; //minimal number of entries found
2968 Double_t maxentries = 0; //maximal number of entries found
2969 Double_t placemin = 0; //calibration group number of the min
2970 Double_t placemax = -1; //calibration group number of the max
2971 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2972 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2974 Double_t counter = 0;
2977 TH1F *nbEntries = 0x0;//distribution of the number of entries
2978 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2979 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2981 // Beginning of the loop over the calibration groups
2982 for (Int_t idect = 0; idect < nbins; idect++) {
2984 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2985 projch->SetDirectory(0);
2987 // Number of entries for this calibration group
2988 Double_t nentries = 0.0;
2990 for (Int_t k = 0; k < nxbins; k++) {
2991 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2995 for (Int_t k = 0; k < nxbins; k++) {
2996 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2997 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2998 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3006 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3007 nbEntries = new TH1F("Number of entries","Number of entries"
3008 ,100,(Int_t)nentries/2,nentries*2);
3009 nbEntries->SetDirectory(0);
3010 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3012 nbEntriesPerGroup->SetDirectory(0);
3013 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3014 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3015 nbEntriesPerSp->SetDirectory(0);
3018 if(nentries > 0) nbEntries->Fill(nentries);
3019 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3020 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3025 if(nentries > maxentries){
3026 maxentries = nentries;
3030 minentries = nentries;
3032 if(nentries < minentries){
3033 minentries = nentries;
3039 meanstats += nentries;
3041 }//calibration groups loop
3043 if(nbwe > 0) meanstats /= nbwe;
3044 if(counter > 0) meanrelativerror /= counter;
3046 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3047 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3048 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3049 AliInfo(Form("The mean number of entries is %f",meanstats));
3050 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3053 info[1] = minentries;
3055 info[3] = maxentries;
3057 info[5] = meanstats;
3058 info[6] = meanrelativerror;
3061 gStyle->SetPalette(1);
3062 gStyle->SetOptStat(1111);
3063 gStyle->SetPadBorderMode(0);
3064 gStyle->SetCanvasColor(10);
3065 gStyle->SetPadLeftMargin(0.13);
3066 gStyle->SetPadRightMargin(0.01);
3067 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3070 nbEntries->Draw("");
3072 nbEntriesPerSp->SetStats(0);
3073 nbEntriesPerSp->Draw("");
3074 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3076 nbEntriesPerGroup->SetStats(0);
3077 nbEntriesPerGroup->Draw("");
3083 //____________________________________________________________________________
3084 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3087 // Return a Int_t[4] with:
3088 // 0 Mean number of entries
3089 // 1 median of number of entries
3090 // 2 rms of number of entries
3091 // 3 number of group with entries
3094 Double_t *stat = new Double_t[4];
3097 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3098 Double_t *weight = new Double_t[nbofgroups];
3099 Int_t *nonul = new Int_t[nbofgroups];
3101 for(Int_t k = 0; k < nbofgroups; k++){
3102 if(fEntriesCH[k] > 0) {
3104 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3107 else weight[k] = 0.0;
3109 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3110 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3111 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3116 //____________________________________________________________________________
3117 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3120 // Return a Int_t[4] with:
3121 // 0 Mean number of entries
3122 // 1 median of number of entries
3123 // 2 rms of number of entries
3124 // 3 number of group with entries
3127 Double_t *stat = new Double_t[4];
3130 Int_t nbofgroups = 540;
3131 Double_t *weight = new Double_t[nbofgroups];
3132 Int_t *nonul = new Int_t[nbofgroups];
3134 for(Int_t k = 0; k < nbofgroups; k++){
3135 if(fEntriesLinearFitter[k] > 0) {
3137 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3140 else weight[k] = 0.0;
3142 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3143 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3144 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3149 //////////////////////////////////////////////////////////////////////////////////////
3151 //////////////////////////////////////////////////////////////////////////////////////
3152 //_____________________________________________________________________________
3153 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3156 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3157 // If fNgroupprf is zero then no binning in tnp
3161 name += fCalibraMode->GetNz(2);
3163 name += fCalibraMode->GetNrphi(2);
3167 if(fNgroupprf != 0){
3169 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3170 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3171 fPRF2d->SetYTitle("Det/pad groups");
3172 fPRF2d->SetXTitle("Position x/W [pad width units]");
3173 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3174 fPRF2d->SetStats(0);
3177 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3178 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3179 fPRF2d->SetYTitle("Det/pad groups");
3180 fPRF2d->SetXTitle("Position x/W [pad width units]");
3181 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3182 fPRF2d->SetStats(0);
3187 //_____________________________________________________________________________
3188 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3191 // Create the 2D histos
3195 name += fCalibraMode->GetNz(1);
3197 name += fCalibraMode->GetNrphi(1);
3199 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3200 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3202 fPH2d->SetYTitle("Det/pad groups");
3203 fPH2d->SetXTitle("time [#mus]");
3204 fPH2d->SetZTitle("<PH> [a.u.]");
3208 //_____________________________________________________________________________
3209 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3212 // Create the 2D histos
3216 name += fCalibraMode->GetNz(0);
3218 name += fCalibraMode->GetNrphi(0);
3220 fCH2d = new TH2I("CH2d",(const Char_t *) name
3221 ,fNumberBinCharge,0,300,nn,0,nn);
3222 fCH2d->SetYTitle("Det/pad groups");
3223 fCH2d->SetXTitle("charge deposit [a.u]");
3224 fCH2d->SetZTitle("counts");
3229 //////////////////////////////////////////////////////////////////////////////////
3230 // Set relative scale
3231 /////////////////////////////////////////////////////////////////////////////////
3232 //_____________________________________________________________________________
3233 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3236 // Set the factor that will divide the deposited charge
3237 // to fit in the histo range [0,300]
3240 if (RelativeScale > 0.0) {
3241 fRelativeScale = RelativeScale;
3244 AliInfo("RelativeScale must be strict positif!");
3248 //////////////////////////////////////////////////////////////////////////////////
3249 // Quick way to fill a histo
3250 //////////////////////////////////////////////////////////////////////////////////
3251 //_____________________________________________________________________
3252 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3255 // FillCH2d: Marian style
3258 //skip simply the value out of range
3259 if((y>=300.0) || (y<0.0)) return;
3261 //Calcul the y place
3262 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3263 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3266 fCH2d->GetArray()[place]++;
3270 //////////////////////////////////////////////////////////////////////////////////
3271 // Geometrical functions
3272 ///////////////////////////////////////////////////////////////////////////////////
3273 //_____________________________________________________________________________
3274 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3277 // Reconstruct the layer number from the detector number
3280 return ((Int_t) (d % 6));
3284 //_____________________________________________________________________________
3285 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3288 // Reconstruct the stack number from the detector number
3290 const Int_t kNlayer = 6;
3292 return ((Int_t) (d % 30) / kNlayer);
3296 //_____________________________________________________________________________
3297 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3300 // Reconstruct the sector number from the detector number
3304 return ((Int_t) (d / fg));
3307 ///////////////////////////////////////////////////////////////////////////////////
3308 // Getter functions for DAQ of the CH2d and the PH2d
3309 //////////////////////////////////////////////////////////////////////////////////
3310 //_____________________________________________________________________
3311 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3314 // return pointer to fPH2d TProfile2D
3315 // create a new TProfile2D if it doesn't exist allready
3321 fTimeMax = nbtimebin;
3322 fSf = samplefrequency;
3328 //_____________________________________________________________________
3329 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3332 // return pointer to fCH2d TH2I
3333 // create a new TH2I if it doesn't exist allready
3342 ////////////////////////////////////////////////////////////////////////////////////////////
3343 // Drift velocity calibration
3344 ///////////////////////////////////////////////////////////////////////////////////////////
3345 //_____________________________________________________________________
3346 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3349 // return pointer to TLinearFitter Calibration
3350 // if force is true create a new TLinearFitter if it doesn't exist allready
3353 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3354 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3357 // if we are forced and TLinearFitter doesn't yet exist create it
3359 // new TLinearFitter
3360 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3361 fLinearFitterArray.AddAt(linearfitter,detector);
3362 return linearfitter;
3365 //____________________________________________________________________________
3366 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3369 // Analyse array of linear fitter because can not be written
3370 // Store two arrays: one with the param the other one with the error param + number of entries
3373 for(Int_t k = 0; k < 540; k++){
3374 TLinearFitter *linearfitter = GetLinearFitter(k);
3375 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3376 TVectorD *par = new TVectorD(2);
3377 TVectorD pare = TVectorD(2);
3378 TVectorD *parE = new TVectorD(3);
3379 linearfitter->Eval();
3380 linearfitter->GetParameters(*par);
3381 linearfitter->GetErrors(pare);
3382 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3383 (*parE)[0] = pare[0]*ppointError;
3384 (*parE)[1] = pare[1]*ppointError;
3385 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3386 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3387 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);