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)
34 //////////////////////////////////////////////////////////////////////////////////////
36 #include <TProfile2D.h>
41 #include <TObjArray.h>
46 #include <TStopwatch.h>
48 #include <TDirectory.h>
49 #include <TTreeStream.h>
54 #include "AliTRDCalibraFillHisto.h"
55 #include "AliTRDCalibraMode.h"
56 #include "AliTRDCalibraVector.h"
57 #include "AliTRDCalibraVdriftLinearFit.h"
58 #include "AliTRDcalibDB.h"
59 #include "AliTRDCommonParam.h"
60 #include "AliTRDmcmTracklet.h"
61 #include "AliTRDmcm.h"
62 #include "AliTRDtrigParam.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcluster.h"
65 #include "AliTRDtrack.h"
66 #include "AliTRDtrackV1.h"
67 #include "AliTRDrawStreamBase.h"
68 #include "AliRawReader.h"
69 #include "AliRawReaderDate.h"
70 #include "AliTRDgeometry.h"
71 #include "./Cal/AliTRDCalROC.h"
72 #include "./Cal/AliTRDCalDet.h"
79 ClassImp(AliTRDCalibraFillHisto)
81 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
82 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
84 //_____________singleton implementation_________________________________________________
85 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
88 // Singleton implementation
91 if (fgTerminated != kFALSE) {
95 if (fgInstance == 0) {
96 fgInstance = new AliTRDCalibraFillHisto();
103 //______________________________________________________________________________________
104 void AliTRDCalibraFillHisto::Terminate()
107 // Singleton implementation
108 // Deletes the instance of this class
111 fgTerminated = kTRUE;
113 if (fgInstance != 0) {
120 //______________________________________________________________________________________
121 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
125 ,fMcmCorrectAngle(kFALSE)
131 ,fLinearFitterOn(kFALSE)
132 ,fLinearFitterDebugOn(kFALSE)
134 ,fThresholdClusterPRF2(15.0)
135 ,fLimitChargeIntegration(kFALSE)
136 ,fFillWithZero(kFALSE)
137 ,fNormalizeNbOfCluster(kFALSE)
140 ,fCalibraMode(new AliTRDCalibraMode())
143 ,fDetectorPreviousTrack(-1)
147 ,fNumberClustersf(30)
153 ,fNumberBinCharge(100)
159 ,fGoodTracklet(kTRUE)
160 ,fLinearFitterTracklet(0x0)
162 ,fEntriesLinearFitter(0x0)
167 ,fLinearFitterArray(540)
168 ,fLinearVdriftFit(0x0)
173 // Default constructor
177 // Init some default values
180 fNumberUsedCh[0] = 0;
181 fNumberUsedCh[1] = 0;
182 fNumberUsedPh[0] = 0;
183 fNumberUsedPh[1] = 0;
185 fGeo = new AliTRDgeometry();
189 //______________________________________________________________________________________
190 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
194 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
197 ,fPRF2dOn(c.fPRF2dOn)
198 ,fHisto2d(c.fHisto2d)
199 ,fVector2d(c.fVector2d)
200 ,fLinearFitterOn(c.fLinearFitterOn)
201 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
202 ,fRelativeScale(c.fRelativeScale)
203 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
204 ,fLimitChargeIntegration(c.fLimitChargeIntegration)
205 ,fFillWithZero(c.fFillWithZero)
206 ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
207 ,fMaxCluster(c.fMaxCluster)
208 ,fNbMaxCluster(c.fNbMaxCluster)
211 ,fDebugLevel(c.fDebugLevel)
212 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
213 ,fMCMPrevious(c.fMCMPrevious)
214 ,fROBPrevious(c.fROBPrevious)
215 ,fNumberClusters(c.fNumberClusters)
216 ,fNumberClustersf(c.fNumberClustersf)
217 ,fProcent(c.fProcent)
218 ,fDifference(c.fDifference)
219 ,fNumberTrack(c.fNumberTrack)
220 ,fTimeMax(c.fTimeMax)
222 ,fNumberBinCharge(c.fNumberBinCharge)
223 ,fNumberBinPRF(c.fNumberBinPRF)
224 ,fNgroupprf(c.fNgroupprf)
228 ,fGoodTracklet(c.fGoodTracklet)
229 ,fLinearFitterTracklet(0x0)
231 ,fEntriesLinearFitter(0x0)
236 ,fLinearFitterArray(540)
237 ,fLinearVdriftFit(0x0)
244 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
245 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
247 fPH2d = new TProfile2D(*c.fPH2d);
248 fPH2d->SetDirectory(0);
251 fPRF2d = new TProfile2D(*c.fPRF2d);
252 fPRF2d->SetDirectory(0);
255 fCH2d = new TH2I(*c.fCH2d);
256 fCH2d->SetDirectory(0);
258 if(c.fLinearVdriftFit){
259 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
262 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
263 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
268 fGeo = new AliTRDgeometry();
271 //____________________________________________________________________________________
272 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
275 // AliTRDCalibraFillHisto destructor
279 if ( fDebugStreamer ) delete fDebugStreamer;
281 if ( fCalDetGain ) delete fCalDetGain;
282 if ( fCalROCGain ) delete fCalROCGain;
284 if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
288 delete [] fEntriesCH;
289 delete [] fEntriesLinearFitter;
292 for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
293 TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
301 //_____________________________________________________________________________
302 void AliTRDCalibraFillHisto::Destroy()
313 //_____________________________________________________________________________
314 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
317 // Delete DebugStreamer
320 if ( fDebugStreamer ) delete fDebugStreamer;
321 fDebugStreamer = 0x0;
324 //_____________________________________________________________________________
325 void AliTRDCalibraFillHisto::ClearHistos()
345 //////////////////////////////////////////////////////////////////////////////////
346 // calibration with AliTRDtrackV1: Init, Update
347 //////////////////////////////////////////////////////////////////////////////////
348 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
349 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
352 // Init the histograms and stuff to be filled
357 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
359 AliInfo("Could not get calibDB");
362 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
364 AliInfo("Could not get CommonParam");
369 fTimeMax = cal->GetNumberOfTimeBins();
370 fSf = parCom->GetSamplingFrequency();
371 if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
372 else fRelativeScale = 1.18;
373 fNumberClustersf = fTimeMax;
374 fNumberClusters = (Int_t)(0.5*fTimeMax);
376 // Init linear fitter
377 if(!fLinearFitterTracklet) {
378 fLinearFitterTracklet = new TLinearFitter(2,"pol1");
379 fLinearFitterTracklet->StoreData(kTRUE);
382 //calib object from database used for reconstruction
384 fCalDetGain->~AliTRDCalDet();
385 new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
386 }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
388 // Calcul Xbins Chambd0, Chamb2
389 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
390 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
391 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
393 // If vector method On initialised all the stuff
395 fCalibraVector = new AliTRDCalibraVector();
396 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
397 fCalibraVector->SetTimeMax(fTimeMax);
398 if(fNgroupprf != 0) {
399 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
400 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
403 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
404 fCalibraVector->SetPRFRange(1.5);
406 for(Int_t k = 0; k < 3; k++){
407 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
408 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
410 TString namech("Nz");
411 namech += fCalibraMode->GetNz(0);
413 namech += fCalibraMode->GetNrphi(0);
414 fCalibraVector->SetNameCH((const char* ) namech);
415 TString nameph("Nz");
416 nameph += fCalibraMode->GetNz(1);
418 nameph += fCalibraMode->GetNrphi(1);
419 fCalibraVector->SetNamePH((const char* ) nameph);
420 TString nameprf("Nz");
421 nameprf += fCalibraMode->GetNz(2);
423 nameprf += fCalibraMode->GetNrphi(2);
425 nameprf += fNgroupprf;
426 fCalibraVector->SetNamePRF((const char* ) nameprf);
429 // Create the 2D histos corresponding to the pad groupCalibration mode
432 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
433 ,fCalibraMode->GetNz(0)
434 ,fCalibraMode->GetNrphi(0)));
436 // Create the 2D histo
441 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
442 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
446 fEntriesCH = new Int_t[ntotal0];
447 for(Int_t k = 0; k < ntotal0; k++){
454 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
455 ,fCalibraMode->GetNz(1)
456 ,fCalibraMode->GetNrphi(1)));
458 // Create the 2D histo
463 fPHPlace = new Short_t[fTimeMax];
464 for (Int_t k = 0; k < fTimeMax; k++) {
467 fPHValue = new Float_t[fTimeMax];
468 for (Int_t k = 0; k < fTimeMax; k++) {
472 if (fLinearFitterOn) {
473 //fLinearFitterArray.Expand(540);
474 fLinearFitterArray.SetName("ArrayLinearFitters");
475 fEntriesLinearFitter = new Int_t[540];
476 for(Int_t k = 0; k < 540; k++){
477 fEntriesLinearFitter[k] = 0;
479 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
484 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
485 ,fCalibraMode->GetNz(2)
486 ,fCalibraMode->GetNrphi(2)));
487 // Create the 2D histo
489 CreatePRF2d(ntotal2);
496 //____________Offline tracking in the AliTRDtracker____________________________
497 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
500 // Use AliTRDtrack for the calibration
504 AliTRDcluster *cl = 0x0; // pointeur to cluster
505 Int_t index0 = 0; // index of the first cluster in the current chamber
506 Int_t index1 = 0; // index of the current cluster in the current chamber
508 //////////////////////////////
509 // loop over the clusters
510 ///////////////////////////////
511 while((cl = t->GetCluster(index1))){
513 /////////////////////////////////////////////////////////////////////////
514 // Fill the infos for the previous clusters if not the same detector
515 ////////////////////////////////////////////////////////////////////////
516 Int_t detector = cl->GetDetector();
517 if ((detector != fDetectorPreviousTrack) &&
518 (index0 != index1)) {
522 //If the same track, then look if the previous detector is in
523 //the same plane, if yes: not a good track
524 if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
528 // Fill only if the track doesn't touch a masked pad or doesn't
531 // drift velocity unables to cut bad tracklets
532 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
536 FillTheInfoOfTheTrackCH(index1-index0);
541 FillTheInfoOfTheTrackPH();
544 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
547 } // if a good tracklet
550 ResetfVariablestracklet();
553 } // Fill at the end the charge
556 /////////////////////////////////////////////////////////////
557 // Localise and take the calib gain object for the detector
558 ////////////////////////////////////////////////////////////
559 if (detector != fDetectorPreviousTrack) {
561 //Localise the detector
562 LocalisationDetectorXbins(detector);
565 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
567 AliInfo("Could not get calibDB");
574 fCalROCGain->~AliTRDCalROC();
575 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
577 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
581 // Reset the detectbjobsor
582 fDetectorPreviousTrack = detector;
584 ///////////////////////////////////////
585 // Store the info of the cluster
586 ///////////////////////////////////////
587 Int_t row = cl->GetPadRow();
588 Int_t col = cl->GetPadCol();
589 CheckGoodTrackletV1(cl);
590 Int_t group[2] = {0,0};
591 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
592 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
593 StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
597 } // while on clusters
599 ///////////////////////////
600 // Fill the last plane
601 //////////////////////////
602 if( index0 != index1 ){
608 // drift velocity unables to cut bad tracklets
609 Bool_t pass = FindP1TrackPHtracklet(t,index0,index1);
613 FillTheInfoOfTheTrackCH(index1-index0);
618 FillTheInfoOfTheTrackPH();
621 if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
623 } // if a good tracklet
628 ResetfVariablestracklet();
633 //____________Offline tracking in the AliTRDtracker____________________________
634 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
637 // Use AliTRDtrackV1 for the calibration
641 const AliTRDseedV1 *tracklet = 0x0; // tracklet per plane
642 AliTRDcluster *cl = 0x0; // cluster attached now to the tracklet
643 Bool_t newtr = kTRUE; // new track
646 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
648 AliInfo("Could not get calibDB");
652 ///////////////////////////
653 // loop over the tracklet
654 ///////////////////////////
655 for(Int_t itr = 0; itr < 6; itr++){
657 if(!(tracklet = t->GetTracklet(itr))) continue;
658 if(!tracklet->IsOK()) continue;
660 ResetfVariablestracklet();
662 //////////////////////////////////////////
663 // localisation of the tracklet and dqdl
664 //////////////////////////////////////////
665 Int_t layer = tracklet->GetPlane();
667 while(!(cl = tracklet->GetClusters(ic++))) continue;
668 Int_t detector = cl->GetDetector();
669 if (detector != fDetectorPreviousTrack) {
670 // if not a new track
672 // don't use the rest of this track if in the same plane
673 if (layer == GetLayer(fDetectorPreviousTrack)) {
674 //printf("bad tracklet, same layer for detector %d\n",detector);
678 //Localise the detector bin
679 LocalisationDetectorXbins(detector);
683 fCalROCGain->~AliTRDCalROC();
684 new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
686 }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
689 fDetectorPreviousTrack = detector;
693 ////////////////////////////
694 // loop over the clusters
695 ////////////////////////////
696 Int_t nbclusters = 0;
697 for(int jc=0; jc<AliTRDseed::knTimebins; jc++){
698 if(!(cl = tracklet->GetClusters(jc))) continue;
701 // Store the info bis of the tracklet
702 Int_t row = cl->GetPadRow();
703 Int_t col = cl->GetPadCol();
704 CheckGoodTrackletV1(cl);
705 Int_t group[2] = {0,0};
706 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,row,col);
707 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,row,col);
708 StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col);
711 ////////////////////////////////////////
712 // Fill the stuffs if a good tracklet
713 ////////////////////////////////////////
716 // drift velocity unables to cut bad tracklets
717 Bool_t pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
721 FillTheInfoOfTheTrackCH(nbclusters);
726 FillTheInfoOfTheTrackPH();
729 if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
731 } // if a good tracklet
737 ///////////////////////////////////////////////////////////////////////////////////
738 // Routine inside the update with AliTRDtrack
739 ///////////////////////////////////////////////////////////////////////////////////
740 //____________Offine tracking in the AliTRDtracker_____________________________
741 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
744 // Drift velocity calibration:
745 // Fit the clusters with a straight line
746 // From the slope find the drift velocity
749 //Number of points: if less than 3 return kFALSE
750 Int_t npoints = index1-index0;
751 if(npoints <= 2) return kFALSE;
756 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
757 // parameters of the track
758 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
759 Double_t tgl = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
760 Double_t tnp = 0.0; // tan angle in the plan xy track
761 if( TMath::Abs(snp) < 1.){
762 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
764 Float_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx now from dz/dl
765 // tilting pad and cross row
766 Int_t crossrow = 0; // if it crosses a pad row
767 Int_t rowp = -1; // if it crosses a pad row
768 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
769 Double_t tiltingangle = padplane->GetTiltingAngle(); // tiltingangle of the pad
770 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
772 fLinearFitterTracklet->ClearPoints();
773 Double_t dydt = 0.0; // dydt tracklet after straight line fit
774 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
775 Double_t pointError = 0.0; // error after straight line fit
776 Int_t nbli = 0; // number linear fitter points
778 //////////////////////////////
779 // loop over clusters
780 ////////////////////////////
781 for(Int_t k = 0; k < npoints; k++){
783 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
784 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
785 Double_t ycluster = cl->GetY();
786 Int_t time = cl->GetPadTime();
787 Double_t timeis = time/fSf;
788 //Double_t q = cl->GetQ();
789 //See if cross two pad rows
790 Int_t row = cl->GetPadRow();
792 if(row != rowp) crossrow = 1;
794 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
799 //////////////////////////////
801 //////////////////////////////
803 fLinearFitterTracklet->ClearPoints();
807 fLinearFitterTracklet->Eval();
808 fLinearFitterTracklet->GetParameters(pars);
809 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
810 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
812 fLinearFitterTracklet->ClearPoints();
814 /////////////////////////////
816 ////////////////////////////
818 if ( !fDebugStreamer ) {
820 TDirectory *backup = gDirectory;
821 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
822 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
826 (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
827 //"snpright="<<snpright<<
828 "npoints="<<npoints<<
830 "detector="<<detector<<
837 "crossrow="<<crossrow<<
838 "errorpar="<<errorpar<<
839 "pointError="<<pointError<<
843 Int_t nbclusters = index1-index0;
844 Int_t layer = GetLayer(fDetectorPreviousTrack);
846 (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
847 //"snpright="<<snpright<<
848 "nbclusters="<<nbclusters<<
849 "detector="<<fDetectorPreviousTrack<<
855 ///////////////////////////
857 ///////////////////////////
858 if(npoints < fNumberClusters) return kFALSE;
859 if(npoints > fNumberClustersf) return kFALSE;
860 if(pointError >= 0.1) return kFALSE;
861 if(crossrow == 1) return kFALSE;
863 ////////////////////////////
865 ////////////////////////////
867 //Add to the linear fitter of the detector
868 if( TMath::Abs(snp) < 1.){
869 Double_t x = tnp-dzdx*tnt;
870 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
871 if(fLinearFitterDebugOn) {
872 fLinearVdriftFit->Update(detector,x,pars[1]);
874 fEntriesLinearFitter[detector]++;
880 //____________Offine tracking in the AliTRDtracker_____________________________
881 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
884 // Drift velocity calibration:
885 // Fit the clusters with a straight line
886 // From the slope find the drift velocity
889 ////////////////////////////////////////////////
890 //Number of points: if less than 3 return kFALSE
891 /////////////////////////////////////////////////
892 if(nbclusters <= 2) return kFALSE;
897 // results of the linear fit
898 Double_t dydt = 0.0; // dydt tracklet after straight line fit
899 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
900 Double_t pointError = 0.0; // error after straight line fit
901 // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
902 Int_t crossrow = 0; // if it crosses a pad row
903 Int_t rowp = -1; // if it crosses a pad row
904 Float_t tnt = tracklet->GetTilt(); // tan tiltingangle
905 fLinearFitterTracklet->ClearPoints();
908 ///////////////////////////////////////////
909 // Take the parameters of the track
910 //////////////////////////////////////////
911 // take now the snp, tnp and tgl from the track
912 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
913 Double_t tnp = 0.0; // dy/dx at the end of the chamber
914 if( TMath::Abs(snp) < 1.){
915 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
917 Double_t tgl = tracklet->GetTgl(); // dz/dl
918 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
920 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
921 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
922 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
923 // at the end with correction due to linear fit
924 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
925 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
928 ////////////////////////////
929 // loop over the clusters
930 ////////////////////////////
932 AliTRDcluster *cl = 0x0;
933 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
934 if(!(cl = tracklet->GetClusters(ic))) continue;
935 if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
937 Double_t ycluster = cl->GetY();
938 Int_t time = cl->GetPadTime();
939 Double_t timeis = time/fSf;
940 //See if cross two pad rows
941 Int_t row = cl->GetPadRow();
942 if(rowp==-1) rowp = row;
943 if(row != rowp) crossrow = 1;
945 fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
950 ////////////////////////////////////
951 // Do the straight line fit now
952 ///////////////////////////////////
954 fLinearFitterTracklet->ClearPoints();
958 fLinearFitterTracklet->Eval();
959 fLinearFitterTracklet->GetParameters(pars);
960 pointError = TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
961 errorpar = fLinearFitterTracklet->GetParError(1)*pointError;
963 //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
964 fLinearFitterTracklet->ClearPoints();
966 ////////////////////////////////
968 ///////////////////////////////
972 if ( !fDebugStreamer ) {
974 TDirectory *backup = gDirectory;
975 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
976 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
980 Int_t layer = GetLayer(fDetectorPreviousTrack);
982 (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
983 //"snpright="<<snpright<<
985 "nbclusters="<<nbclusters<<
986 "detector="<<fDetectorPreviousTrack<<
994 "crossrow="<<crossrow<<
995 "errorpar="<<errorpar<<
996 "pointError="<<pointError<<
1001 /////////////////////////
1003 ////////////////////////
1005 if(nbclusters < fNumberClusters) return kFALSE;
1006 if(nbclusters > fNumberClustersf) return kFALSE;
1007 if(pointError >= 0.3) return kFALSE;
1008 if(crossrow == 1) return kFALSE;
1010 ///////////////////////
1012 //////////////////////
1014 if(fLinearFitterOn){
1015 //Add to the linear fitter of the detector
1016 if( TMath::Abs(snp) < 1.){
1017 Double_t x = tnp-dzdx*tnt;
1018 (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1019 if(fLinearFitterDebugOn) {
1020 fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1022 fEntriesLinearFitter[fDetectorPreviousTrack]++;
1028 //____________Offine tracking in the AliTRDtracker_____________________________
1029 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1032 // PRF width calibration
1033 // Assume a Gaussian shape: determinate the position of the three pad clusters
1034 // Fit with a straight line
1035 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1036 // Fill the PRF as function of angle of the track
1041 //////////////////////////
1043 /////////////////////////
1044 Int_t npoints = index1-index0; // number of total points
1045 Int_t nb3pc = 0; // number of three pads clusters used for fit
1046 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1047 // To see the difference due to the fit
1048 Double_t *padPositions;
1049 padPositions = new Double_t[npoints];
1050 for(Int_t k = 0; k < npoints; k++){
1051 padPositions[k] = 0.0;
1053 // Take the tgl and snp with the track t now
1054 Double_t tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1055 Double_t snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1056 Float_t dzdx = 0.0; // dzdx
1058 if(TMath::Abs(snp) < 1.0){
1059 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1060 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1063 fLinearFitterTracklet->ClearPoints();
1065 ///////////////////////////
1066 // calcul the tnp group
1067 ///////////////////////////
1068 Bool_t echec = kFALSE;
1069 Double_t shift = 0.0;
1070 //Calculate the shift in x coresponding to this tnp
1071 if(fNgroupprf != 0.0){
1072 shift = -3.0*(fNgroupprf-1)-1.5;
1073 Double_t limithigh = -0.2*(fNgroupprf-1);
1074 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1076 while(tnp > limithigh){
1083 delete [] padPositions;
1087 //////////////////////
1089 /////////////////////
1090 for(Int_t k = 0; k < npoints; k++){
1092 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1093 Short_t *signals = cl->GetSignals();
1094 Double_t time = cl->GetLocalTimeBin();
1095 //Calculate x if possible
1096 Float_t xcenter = 0.0;
1097 Bool_t echec1 = kTRUE;
1098 if((time<=7) || (time>=21)) continue;
1099 // Center 3 balanced: position with the center of the pad
1100 if ((((Float_t) signals[3]) > 0.0) &&
1101 (((Float_t) signals[2]) > 0.0) &&
1102 (((Float_t) signals[4]) > 0.0)) {
1104 // Security if the denomiateur is 0
1105 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1106 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1107 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1108 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1109 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1115 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1117 //if no echec: calculate with the position of the pad
1118 // Position of the cluster
1119 Double_t padPosition = xcenter + cl->GetPadCol();
1120 padPositions[k] = padPosition;
1122 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1126 /////////////////////////////
1128 ////////////////////////////
1130 delete [] padPositions;
1131 fLinearFitterTracklet->ClearPoints();
1134 fLinearFitterTracklet->Eval();
1136 fLinearFitterTracklet->GetParameters(line);
1137 Float_t pointError = -1.0;
1138 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1139 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1141 fLinearFitterTracklet->ClearPoints();
1143 /////////////////////////////////////////////////////
1144 // Now fill the PRF: second loop over clusters
1145 /////////////////////////////////////////////////////
1146 for(Int_t k = 0; k < npoints; k++){
1148 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
1149 Short_t *signals = cl->GetSignals(); // signal
1150 Double_t time = cl->GetLocalTimeBin(); // time bin
1151 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1152 Float_t padPos = cl->GetPadCol(); // middle pad
1153 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1154 Float_t ycenter = 0.0; // relative center charge
1155 Float_t ymin = 0.0; // relative left charge
1156 Float_t ymax = 0.0; // relative right charge
1158 //Requiere simply two pads clusters at least
1159 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1160 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1161 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1162 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1163 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1164 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1168 Int_t rowcl = cl->GetPadRow(); // row of cluster
1169 Int_t colcl = cl->GetPadCol(); // col of cluster
1170 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1171 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1172 Float_t xcl = cl->GetY(); // y cluster
1173 Float_t qcl = cl->GetQ(); // charge cluster
1174 Int_t layer = GetLayer(detector); // layer
1175 Int_t stack = GetStack(detector); // stack
1176 Double_t xdiff = dpad; // reconstructed position constant
1177 Double_t x = dpad; // reconstructed position moved
1178 Float_t ep = pointError; // error of fit
1179 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1180 Float_t signal3 = (Float_t)signals[3]; // signal
1181 Float_t signal2 = (Float_t)signals[2]; // signal
1182 Float_t signal4 = (Float_t)signals[4]; // signal
1183 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1185 //////////////////////////////
1187 /////////////////////////////
1188 if(fDebugLevel > 0){
1189 if ( !fDebugStreamer ) {
1191 TDirectory *backup = gDirectory;
1192 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1193 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1199 Float_t y = ycenter;
1200 (* fDebugStreamer) << "HandlePRFtracklet"<<
1201 "caligroup="<<caligroup<<
1202 "detector="<<detector<<
1205 "npoints="<<npoints<<
1214 "padPosition="<<padPositions[k]<<
1215 "padPosTracklet="<<padPosTracklet<<
1220 "signal1="<<signal1<<
1221 "signal2="<<signal2<<
1222 "signal3="<<signal3<<
1223 "signal4="<<signal4<<
1224 "signal5="<<signal5<<
1230 (* fDebugStreamer) << "HandlePRFtracklet"<<
1231 "caligroup="<<caligroup<<
1232 "detector="<<detector<<
1235 "npoints="<<npoints<<
1244 "padPosition="<<padPositions[k]<<
1245 "padPosTracklet="<<padPosTracklet<<
1250 "signal1="<<signal1<<
1251 "signal2="<<signal2<<
1252 "signal3="<<signal3<<
1253 "signal4="<<signal4<<
1254 "signal5="<<signal5<<
1260 (* fDebugStreamer) << "HandlePRFtracklet"<<
1261 "caligroup="<<caligroup<<
1262 "detector="<<detector<<
1265 "npoints="<<npoints<<
1274 "padPosition="<<padPositions[k]<<
1275 "padPosTracklet="<<padPosTracklet<<
1280 "signal1="<<signal1<<
1281 "signal2="<<signal2<<
1282 "signal3="<<signal3<<
1283 "signal4="<<signal4<<
1284 "signal5="<<signal5<<
1290 ////////////////////////////
1292 ///////////////////////////
1293 if(npoints < fNumberClusters) continue;
1294 if(npoints > fNumberClustersf) continue;
1295 if(nb3pc <= 5) continue;
1296 if((time >= 21) || (time < 7)) continue;
1297 if(TMath::Abs(snp) >= 1.0) continue;
1298 if(TMath::Abs(qcl) < 80) continue;
1300 ////////////////////////////
1302 ///////////////////////////
1304 if(TMath::Abs(dpad) < 1.5) {
1305 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1306 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1308 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1309 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1310 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1312 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1313 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1314 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1318 if(TMath::Abs(dpad) < 1.5) {
1319 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1320 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1322 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1323 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1324 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1326 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1327 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1328 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1332 delete [] padPositions;
1336 //____________Offine tracking in the AliTRDtracker_____________________________
1337 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1340 // PRF width calibration
1341 // Assume a Gaussian shape: determinate the position of the three pad clusters
1342 // Fit with a straight line
1343 // Take the fitted values for all the clusters (3 or 2 pad clusters)
1344 // Fill the PRF as function of angle of the track
1348 //printf("begin\n");
1349 ///////////////////////////////////////////
1350 // Take the parameters of the track
1351 //////////////////////////////////////////
1352 // take now the snp, tnp and tgl from the track
1353 Double_t snp = tracklet->GetSnp(); // sin dy/dx at the end of the chamber
1354 Double_t tnp = 0.0; // dy/dx at the end of the chamber
1355 if( TMath::Abs(snp) < 1.){
1356 tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1358 Double_t tgl = tracklet->GetTgl(); // dz/dl
1359 Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp); // dz/dx calculated from dz/dl
1361 //Double_t tnp = tracklet->GetYref(1); // dy/dx at the entrance of the chamber
1362 //Double_t tgl = tracklet->GetZref(1); // dz/dl at the entrance of the chamber
1363 //Double_t dzdx = tgl; //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1364 // at the end with correction due to linear fit
1365 //Double_t tnp = tracklet->GetYfit(1); // dy/dx at the end of the chamber after fit correction
1366 //Double_t tgl = tracklet->GetZfit(1); // dz/dl at the end of the chamber after fit correction
1368 ///////////////////////////////
1369 // Calculate tnp group shift
1370 ///////////////////////////////
1371 Bool_t echec = kFALSE;
1372 Double_t shift = 0.0;
1373 //Calculate the shift in x coresponding to this tnp
1374 if(fNgroupprf != 0.0){
1375 shift = -3.0*(fNgroupprf-1)-1.5;
1376 Double_t limithigh = -0.2*(fNgroupprf-1);
1377 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1379 while(tnp > limithigh){
1385 // do nothing if out of tnp range
1386 //printf("echec %d\n",(Int_t)echec);
1387 if(echec) return kFALSE;
1389 ///////////////////////
1391 //////////////////////
1393 Int_t nb3pc = 0; // number of three pads clusters used for fit
1394 Double_t *padPositions; // to see the difference between the fit and the 3 pad clusters position
1395 padPositions = new Double_t[AliTRDseed::knTimebins];
1396 for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1397 padPositions[k] = 0.0;
1399 fLinearFitterTracklet->ClearPoints();
1401 //printf("loop clusters \n");
1402 ////////////////////////////
1403 // loop over the clusters
1404 ////////////////////////////
1405 AliTRDcluster *cl = 0x0;
1406 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1407 if(!(cl = tracklet->GetClusters(ic))) continue;
1409 Double_t time = cl->GetLocalTimeBin();
1410 if((time<=7) || (time>=21)) continue;
1411 Short_t *signals = cl->GetSignals();
1412 Float_t xcenter = 0.0;
1413 Bool_t echec1 = kTRUE;
1415 /////////////////////////////////////////////////////////////
1416 // Center 3 balanced: position with the center of the pad
1417 /////////////////////////////////////////////////////////////
1418 if ((((Float_t) signals[3]) > 0.0) &&
1419 (((Float_t) signals[2]) > 0.0) &&
1420 (((Float_t) signals[4]) > 0.0)) {
1422 // Security if the denomiateur is 0
1423 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1424 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1425 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1426 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1427 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1433 if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1434 if(echec1) continue;
1436 ////////////////////////////////////////////////////////
1437 //if no echec1: calculate with the position of the pad
1438 // Position of the cluster
1439 // fill the linear fitter
1440 ///////////////////////////////////////////////////////
1441 Double_t padPosition = xcenter + cl->GetPadCol();
1442 padPositions[ic] = padPosition;
1444 fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1449 //printf("Fin loop clusters \n");
1450 //////////////////////////////
1451 // fit with a straight line
1452 /////////////////////////////
1454 delete [] padPositions;
1455 fLinearFitterTracklet->ClearPoints();
1456 delete [] padPositions;
1459 fLinearFitterTracklet->Eval();
1461 fLinearFitterTracklet->GetParameters(line);
1462 Float_t pointError = -1.0;
1463 if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1464 pointError = TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1466 fLinearFitterTracklet->ClearPoints();
1468 //printf("PRF second loop \n");
1469 ////////////////////////////////////////////////
1470 // Fill the PRF: Second loop over clusters
1471 //////////////////////////////////////////////
1472 for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1473 if(!(cl = tracklet->GetClusters(ic))) continue;
1475 Short_t *signals = cl->GetSignals(); // signal
1476 Double_t time = cl->GetLocalTimeBin(); // time bin
1477 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1478 Float_t padPos = cl->GetPadCol(); // middle pad
1479 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1480 Float_t ycenter = 0.0; // relative center charge
1481 Float_t ymin = 0.0; // relative left charge
1482 Float_t ymax = 0.0; // relative right charge
1484 ////////////////////////////////////////////////////////////////
1485 // Calculate ycenter, ymin and ymax even for two pad clusters
1486 ////////////////////////////////////////////////////////////////
1487 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1488 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1489 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1490 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1491 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
1492 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
1495 /////////////////////////
1496 // Calibration group
1497 ////////////////////////
1498 Int_t rowcl = cl->GetPadRow(); // row of cluster
1499 Int_t colcl = cl->GetPadCol(); // col of cluster
1500 Int_t grouplocal = CalculateCalibrationGroup(2,rowcl,colcl); // calcul the corresponding group
1501 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
1502 Float_t xcl = cl->GetY(); // y cluster
1503 Float_t qcl = cl->GetQ(); // charge cluster
1504 Int_t layer = GetLayer(fDetectorPreviousTrack); // layer
1505 Int_t stack = GetStack(fDetectorPreviousTrack); // stack
1506 Double_t xdiff = dpad; // reconstructed position constant
1507 Double_t x = dpad; // reconstructed position moved
1508 Float_t ep = pointError; // error of fit
1509 Float_t signal1 = (Float_t)signals[1]; // signal at the border
1510 Float_t signal3 = (Float_t)signals[3]; // signal
1511 Float_t signal2 = (Float_t)signals[2]; // signal
1512 Float_t signal4 = (Float_t)signals[4]; // signal
1513 Float_t signal5 = (Float_t)signals[5]; // signal at the border
1517 /////////////////////
1519 ////////////////////
1521 if(fDebugLevel > 0){
1522 if ( !fDebugStreamer ) {
1524 TDirectory *backup = gDirectory;
1525 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1526 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1531 Float_t y = ycenter;
1532 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1533 "caligroup="<<caligroup<<
1534 "detector="<<fDetectorPreviousTrack<<
1537 "npoints="<<nbclusters<<
1546 "padPosition="<<padPositions[ic]<<
1547 "padPosTracklet="<<padPosTracklet<<
1552 "signal1="<<signal1<<
1553 "signal2="<<signal2<<
1554 "signal3="<<signal3<<
1555 "signal4="<<signal4<<
1556 "signal5="<<signal5<<
1562 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1563 "caligroup="<<caligroup<<
1564 "detector="<<fDetectorPreviousTrack<<
1567 "npoints="<<nbclusters<<
1576 "padPosition="<<padPositions[ic]<<
1577 "padPosTracklet="<<padPosTracklet<<
1582 "signal1="<<signal1<<
1583 "signal2="<<signal2<<
1584 "signal3="<<signal3<<
1585 "signal4="<<signal4<<
1586 "signal5="<<signal5<<
1592 (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1593 "caligroup="<<caligroup<<
1594 "detector="<<fDetectorPreviousTrack<<
1597 "npoints="<<nbclusters<<
1606 "padPosition="<<padPositions[ic]<<
1607 "padPosTracklet="<<padPosTracklet<<
1612 "signal1="<<signal1<<
1613 "signal2="<<signal2<<
1614 "signal3="<<signal3<<
1615 "signal4="<<signal4<<
1616 "signal5="<<signal5<<
1622 /////////////////////
1624 /////////////////////
1625 if(nbclusters < fNumberClusters) continue;
1626 if(nbclusters > fNumberClustersf) continue;
1627 if(nb3pc <= 5) continue;
1628 if((time >= 21) || (time < 7)) continue;
1629 if(TMath::Abs(qcl) < 80) continue;
1630 if( TMath::Abs(snp) > 1.) continue;
1633 ////////////////////////
1635 ///////////////////////
1637 if(TMath::Abs(dpad) < 1.5) {
1638 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1639 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1640 //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1642 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1643 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1644 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1646 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1647 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1648 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1653 if(TMath::Abs(dpad) < 1.5) {
1654 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1655 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1657 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1658 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1659 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1661 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1662 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1663 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1666 } // second loop over clusters
1669 delete [] padPositions;
1673 ///////////////////////////////////////////////////////////////////////////////////////
1674 // Pad row col stuff: see if masked or not
1675 ///////////////////////////////////////////////////////////////////////////////////////
1676 //_____________________________________________________________________________
1677 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1680 // See if we are not near a masked pad
1683 if(cl->IsMasked()) fGoodTracklet = kFALSE;
1687 //_____________________________________________________________________________
1688 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1691 // See if we are not near a masked pad
1694 if (!IsPadOn(detector, col, row)) {
1695 fGoodTracklet = kFALSE;
1699 if (!IsPadOn(detector, col-1, row)) {
1700 fGoodTracklet = kFALSE;
1705 if (!IsPadOn(detector, col+1, row)) {
1706 fGoodTracklet = kFALSE;
1711 //_____________________________________________________________________________
1712 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1715 // Look in the choosen database if the pad is On.
1716 // If no the track will be "not good"
1719 // Get the parameter object
1720 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1722 AliInfo("Could not get calibDB");
1726 if (!cal->IsChamberInstalled(detector) ||
1727 cal->IsChamberMasked(detector) ||
1728 cal->IsPadMasked(detector,col,row)) {
1736 ///////////////////////////////////////////////////////////////////////////////////////
1737 // Calibration groups: calculate the number of groups, localise...
1738 ////////////////////////////////////////////////////////////////////////////////////////
1739 //_____________________________________________________________________________
1740 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1743 // Calculate the calibration group number for i
1746 // Row of the cluster and position in the pad groups
1748 if (fCalibraMode->GetNnZ(i) != 0) {
1749 posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1753 // Col of the cluster and position in the pad groups
1755 if (fCalibraMode->GetNnRphi(i) != 0) {
1756 posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1759 return posc*fCalibraMode->GetNfragZ(i)+posr;
1762 //____________________________________________________________________________________
1763 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1766 // Calculate the total number of calibration groups
1772 if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1774 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1779 if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1781 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1786 fCalibraMode->ModePadCalibration(2,i);
1787 fCalibraMode->ModePadFragmentation(0,2,0,i);
1788 fCalibraMode->SetDetChamb2(i);
1789 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1790 fCalibraMode->ModePadCalibration(0,i);
1791 fCalibraMode->ModePadFragmentation(0,0,0,i);
1792 fCalibraMode->SetDetChamb0(i);
1793 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1794 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1798 //_____________________________________________________________________________
1799 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1802 // Set the mode of calibration group in the z direction for the parameter i
1807 fCalibraMode->SetNz(i, Nz);
1810 AliInfo("You have to choose between 0 and 4");
1814 //_____________________________________________________________________________
1815 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1818 // Set the mode of calibration group in the rphi direction for the parameter i
1823 fCalibraMode->SetNrphi(i ,Nrphi);
1826 AliInfo("You have to choose between 0 and 6");
1831 //_____________________________________________________________________________
1832 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1835 // Set the mode of calibration group all together
1837 if(fVector2d == kTRUE) {
1838 AliInfo("Can not work with the vector method");
1841 fCalibraMode->SetAllTogether(i);
1845 //_____________________________________________________________________________
1846 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1849 // Set the mode of calibration group per supermodule
1851 if(fVector2d == kTRUE) {
1852 AliInfo("Can not work with the vector method");
1855 fCalibraMode->SetPerSuperModule(i);
1859 //____________Set the pad calibration variables for the detector_______________
1860 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1863 // For the detector calcul the first Xbins and set the number of row
1864 // and col pads per calibration groups, the number of calibration
1865 // groups in the detector.
1868 // first Xbins of the detector
1870 fCalibraMode->CalculXBins(detector,0);
1873 fCalibraMode->CalculXBins(detector,1);
1876 fCalibraMode->CalculXBins(detector,2);
1879 // fragmentation of idect
1880 for (Int_t i = 0; i < 3; i++) {
1881 fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1882 fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1883 , (Int_t) GetStack(detector)
1884 , (Int_t) GetSector(detector),i);
1890 //_____________________________________________________________________________
1891 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1894 // Should be between 0 and 6
1897 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1898 AliInfo("The number of groups must be between 0 and 6!");
1901 fNgroupprf = numberGroupsPRF;
1905 ///////////////////////////////////////////////////////////////////////////////////////////
1906 // Per tracklet: store or reset the info, fill the histos with the info
1907 //////////////////////////////////////////////////////////////////////////////////////////
1908 //_____________________________________________________________________________
1909 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1912 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1913 // Correct from the gain correction before
1916 // time bin of the cluster not corrected
1917 Int_t time = cl->GetPadTime();
1919 //Correct for the gain coefficient used in the database for reconstruction
1920 Float_t correctthegain = 1.0;
1921 if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1922 else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1923 Float_t correction = 1.0;
1924 Float_t normalisation = 6.67;
1925 // we divide with gain in AliTRDclusterizer::Transform...
1926 if( correctthegain > 0 ) normalisation /= correctthegain;
1929 // take dd/dl corrected from the angle of the track
1930 correction = dqdl / (normalisation);
1933 // Fill the fAmpTotal with the charge
1935 if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1938 // Fill the fPHPlace and value
1940 fPHPlace[time] = group[1];
1941 fPHValue[time] = correction;
1945 //____________Offine tracking in the AliTRDtracker_____________________________
1946 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1949 // Reset values per tracklet
1952 //Reset good tracklet
1953 fGoodTracklet = kTRUE;
1955 // Reset the fPHValue
1957 //Reset the fPHValue and fPHPlace
1958 for (Int_t k = 0; k < fTimeMax; k++) {
1964 // Reset the fAmpTotal where we put value
1966 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1971 //____________Offine tracking in the AliTRDtracker_____________________________
1972 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1975 // For the offline tracking or mcm tracklets
1976 // This function will be called in the functions UpdateHistogram...
1977 // to fill the info of a track for the relativ gain calibration
1980 Int_t nb = 0; // Nombre de zones traversees
1981 Int_t fd = -1; // Premiere zone non nulle
1982 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
1984 if(nbclusters < fNumberClusters) return;
1985 if(nbclusters > fNumberClustersf) return;
1988 // Normalize with the number of clusters
1989 Double_t normalizeCst = fRelativeScale;
1990 if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1992 // See if the track goes through different zones
1993 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1994 if (fAmpTotal[k] > 0.0) {
1995 totalcharge += fAmpTotal[k];
2008 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2010 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2011 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2014 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2018 if ((fAmpTotal[fd] > 0.0) &&
2019 (fAmpTotal[fd+1] > 0.0)) {
2020 // One of the two very big
2021 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2023 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2024 //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2027 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2030 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2032 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2034 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2035 //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2038 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2041 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2044 if (fCalibraMode->GetNfragZ(0) > 1) {
2045 if (fAmpTotal[fd] > 0.0) {
2046 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2047 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2048 // One of the two very big
2049 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
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);
2058 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2060 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2062 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2063 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2066 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2068 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2079 //____________Offine tracking in the AliTRDtracker_____________________________
2080 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2083 // For the offline tracking or mcm tracklets
2084 // This function will be called in the functions UpdateHistogram...
2085 // to fill the info of a track for the drift velocity calibration
2088 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2089 Int_t fd1 = -1; // Premiere zone non nulle
2090 Int_t fd2 = -1; // Deuxieme zone non nulle
2091 Int_t k1 = -1; // Debut de la premiere zone
2092 Int_t k2 = -1; // Debut de la seconde zone
2093 Int_t nbclusters = 0; // number of clusters
2097 // See if the track goes through different zones
2098 for (Int_t k = 0; k < fTimeMax; k++) {
2099 if (fPHValue[k] > 0.0) {
2105 if (fPHPlace[k] != fd1) {
2111 if (fPHPlace[k] != fd2) {
2118 // See if noise before and after
2119 if(fMaxCluster > 0) {
2120 if(fPHValue[0] > fMaxCluster) return;
2121 if(fTimeMax > fNbMaxCluster) {
2122 for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2123 if(fPHValue[k] > fMaxCluster) return;
2128 //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2130 if(nbclusters < fNumberClusters) return;
2131 if(nbclusters > fNumberClustersf) return;
2137 for (Int_t i = 0; i < fTimeMax; i++) {
2139 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2141 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2143 //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2146 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2148 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2154 if ((fd1 == fd2+1) ||
2156 // One of the two fast all the think
2157 if (k2 > (k1+fDifference)) {
2158 //we choose to fill the fd1 with all the values
2160 for (Int_t i = 0; i < fTimeMax; i++) {
2162 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2164 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2168 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2170 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2175 if ((k2+fDifference) < fTimeMax) {
2176 //we choose to fill the fd2 with all the values
2178 for (Int_t i = 0; i < fTimeMax; i++) {
2180 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2182 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2186 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2188 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2194 // Two zones voisines sinon rien!
2195 if (fCalibraMode->GetNfragZ(1) > 1) {
2197 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2198 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2199 // One of the two fast all the think
2200 if (k2 > (k1+fDifference)) {
2201 //we choose to fill the fd1 with all the values
2203 for (Int_t i = 0; i < fTimeMax; i++) {
2205 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2207 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2211 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2213 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2218 if ((k2+fDifference) < fTimeMax) {
2219 //we choose to fill the fd2 with all the values
2221 for (Int_t i = 0; i < fTimeMax; i++) {
2223 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2225 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2229 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2231 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2238 // Two zones voisines sinon rien!
2240 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2241 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2242 // One of the two fast all the think
2243 if (k2 > (k1 + fDifference)) {
2244 //we choose to fill the fd1 with all the values
2246 for (Int_t i = 0; i < fTimeMax; i++) {
2248 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2250 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2254 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2256 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2261 if ((k2+fDifference) < fTimeMax) {
2262 //we choose to fill the fd2 with all the values
2264 for (Int_t i = 0; i < fTimeMax; i++) {
2266 if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2268 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2272 if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2274 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2286 //////////////////////////////////////////////////////////////////////////////////////////
2287 // DAQ process functions
2288 /////////////////////////////////////////////////////////////////////////////////////////
2289 //_____________________________________________________________________
2290 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2293 // Event Processing loop - AliTRDrawStreamBase
2294 // TestBeam 2007 version
2295 // 0 timebin problem
2299 // Same algorithm as TestBeam but different reader
2302 Int_t withInput = 1;
2304 Double_t phvalue[16][144][36];
2305 for(Int_t k = 0; k < 36; k++){
2306 for(Int_t j = 0; j < 16; j++){
2307 for(Int_t c = 0; c < 144; c++){
2308 phvalue[j][c][k] = 0.0;
2313 fDetectorPreviousTrack = -1;
2316 Int_t nbtimebin = 0;
2317 Int_t baseline = 10;
2324 while (rawStream->Next()) {
2326 Int_t idetector = rawStream->GetDet(); // current detector
2327 Int_t imcm = rawStream->GetMCM(); // current MCM
2328 Int_t irob = rawStream->GetROB(); // current ROB
2330 //printf("Detector %d\n",idetector);
2332 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2335 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2339 for(Int_t k = 0; k < 36; k++){
2340 for(Int_t j = 0; j < 16; j++){
2341 for(Int_t c = 0; c < 144; c++){
2342 phvalue[j][c][k] = 0.0;
2348 fDetectorPreviousTrack = idetector;
2349 fMCMPrevious = imcm;
2350 fROBPrevious = irob;
2352 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2353 if(nbtimebin == 0) return 0;
2354 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2355 fTimeMax = nbtimebin;
2357 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2358 fNumberClustersf = fTimeMax;
2359 fNumberClusters = (Int_t)(0.6*fTimeMax);
2362 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2363 Int_t col = rawStream->GetCol();
2364 Int_t row = rawStream->GetRow();
2367 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2370 for(Int_t itime = 0; itime < nbtimebin; itime++){
2371 phvalue[row][col][itime] = signal[itime]-baseline;
2375 // fill the last one
2376 if(fDetectorPreviousTrack != -1){
2379 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2382 for(Int_t k = 0; k < 36; k++){
2383 for(Int_t j = 0; j < 16; j++){
2384 for(Int_t c = 0; c < 144; c++){
2385 phvalue[j][c][k] = 0.0;
2394 while (rawStream->Next()) {
2396 Int_t idetector = rawStream->GetDet(); // current detector
2397 Int_t imcm = rawStream->GetMCM(); // current MCM
2398 Int_t irob = rawStream->GetROB(); // current ROB
2400 //printf("Detector %d\n",idetector);
2402 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2405 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2408 for(Int_t k = 0; k < 36; k++){
2409 for(Int_t j = 0; j < 16; j++){
2410 for(Int_t c = 0; c < 144; c++){
2411 phvalue[j][c][k] = 0.0;
2417 fDetectorPreviousTrack = idetector;
2418 fMCMPrevious = imcm;
2419 fROBPrevious = irob;
2421 //baseline = rawStream->GetCommonAdditive(); // common baseline
2423 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2424 fNumberClustersf = fTimeMax;
2425 fNumberClusters = (Int_t)(0.6*fTimeMax);
2426 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2427 Int_t col = rawStream->GetCol();
2428 Int_t row = rawStream->GetRow();
2431 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2433 for(Int_t itime = 0; itime < fTimeMax; itime++){
2434 phvalue[row][col][itime] = signal[itime]-baseline;
2438 // fill the last one
2439 if(fDetectorPreviousTrack != -1){
2442 withInput = TMath::Max(FillDAQ(phvalue),withInput);
2445 for(Int_t k = 0; k < 36; k++){
2446 for(Int_t j = 0; j < 16; j++){
2447 for(Int_t c = 0; c < 144; c++){
2448 phvalue[j][c][k] = 0.0;
2458 //_____________________________________________________________________
2459 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2462 // Event Processing loop - AliTRDrawStreamBase
2463 // Use old AliTRDmcmtracklet code
2464 // 0 timebin problem
2468 // Algorithm with mcm tracklet
2471 Int_t withInput = 1;
2473 AliTRDmcm mcm = AliTRDmcm(0);
2474 AliTRDtrigParam *param = AliTRDtrigParam::Instance();
2475 rawStream->SetSharedPadReadout(kTRUE);
2477 fDetectorPreviousTrack = -1;
2481 Int_t nbtimebin = 0;
2482 Int_t baseline = 10;
2489 while (rawStream->Next()) {
2491 Int_t idetector = rawStream->GetDet(); // current detector
2492 Int_t imcm = rawStream->GetMCM(); // current MCM
2493 Int_t irob = rawStream->GetROB(); // current ROB
2494 row = rawStream->GetRow();
2497 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2500 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2506 mcm.SetChaId(idetector);
2508 mcm.SetColRange(0,21);
2511 if(fDetectorPreviousTrack == -1){
2514 mcm.SetChaId(idetector);
2516 mcm.SetColRange(0,21);
2520 fDetectorPreviousTrack = idetector;
2521 fMCMPrevious = imcm;
2522 fROBPrevious = irob;
2524 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2525 if(nbtimebin == 0) return 0;
2526 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2527 fTimeMax = nbtimebin;
2528 fNumberClustersf = fTimeMax;
2529 fNumberClusters = (Int_t)(0.6*fTimeMax);
2530 param->SetTimeRange(0,fTimeMax);
2532 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
2534 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2535 Int_t adc = rawStream->GetADC();
2538 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2540 for(Int_t itime = 0; itime < nbtimebin; itime++){
2541 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2545 // fill the last one
2546 if(fDetectorPreviousTrack != -1){
2549 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2554 mcm.SetRobId(fROBPrevious);
2555 mcm.SetChaId(fDetectorPreviousTrack);
2557 mcm.SetColRange(0,21);
2564 while (rawStream->Next()) {
2566 Int_t idetector = rawStream->GetDet(); // current detector
2567 Int_t imcm = rawStream->GetMCM(); // current MCM
2568 Int_t irob = rawStream->GetROB(); // current ROB
2569 row = rawStream->GetRow();
2571 if(((fDetectorPreviousTrack != idetector) || (fMCMPrevious != imcm) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
2574 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2580 mcm.SetChaId(idetector);
2582 mcm.SetColRange(0,21);
2586 fDetectorPreviousTrack = idetector;
2587 fMCMPrevious = imcm;
2588 fROBPrevious = irob;
2590 //baseline = rawStream->GetCommonAdditive(); // common baseline
2592 fTimeMax = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
2593 fNumberClustersf = fTimeMax;
2594 fNumberClusters = (Int_t)(0.6*fTimeMax);
2595 param->SetTimeRange(0,fTimeMax);
2596 Int_t *signal = rawStream->GetSignals(); // current ADC signal
2597 Int_t adc = rawStream->GetADC();
2600 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2602 for(Int_t itime = 0; itime < fTimeMax; itime++){
2603 mcm.SetADC(adc,itime,(signal[itime]-baseline));
2607 // fill the last one
2608 if(fDetectorPreviousTrack != -1){
2611 withInput = TMath::Max(FillDAQ(&mcm),withInput);
2615 mcm.SetRobId(fROBPrevious);
2616 mcm.SetChaId(fDetectorPreviousTrack);
2618 mcm.SetColRange(0,21);
2626 //_____________________________________________________________________
2627 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2630 // Event processing loop - AliRawReader
2631 // Testbeam 2007 version
2634 AliTRDrawStreamBase rawStream(rawReader);
2636 rawReader->Select("TRD");
2638 return ProcessEventDAQ(&rawStream, nocheck);
2641 //_________________________________________________________________________
2642 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2644 eventHeaderStruct *event,
2647 eventHeaderStruct* /*event*/,
2654 // process date event
2655 // Testbeam 2007 version
2658 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2659 Int_t result=ProcessEventDAQ(rawReader, nocheck);
2663 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2668 //_____________________________________________________________________
2669 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(AliRawReader *rawReader, Bool_t nocheck)
2672 // Event processing loop - AliRawReader
2673 // use the old mcm traklet code
2676 AliTRDrawStreamBase rawStream(rawReader);
2678 rawReader->Select("TRD");
2680 return ProcessEventDAQV1(&rawStream, nocheck);
2682 //_________________________________________________________________________
2683 Int_t AliTRDCalibraFillHisto::ProcessEventDAQV1(
2685 eventHeaderStruct *event,
2688 eventHeaderStruct* /*event*/,
2695 // process date event
2696 // use the old mcm tracklet code
2699 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2700 Int_t result=ProcessEventDAQV1(rawReader, nocheck);
2704 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2709 //////////////////////////////////////////////////////////////////////////////
2710 // Routine inside the DAQ process
2711 /////////////////////////////////////////////////////////////////////////////
2712 //_______________________________________________________________________
2713 Int_t AliTRDCalibraFillHisto::FillDAQ(AliTRDmcm *mcm){
2716 // Return 2 if some tracklets are found and used, 1 if nothing
2723 for (Int_t iSeed = 0; iSeed < 4; iSeed++) {
2725 if (mcm->GetSeedCol()[iSeed] < 0) {
2729 nbev += TestTracklet(mcm->GetChaId(),mcm->GetRow(),iSeed,mcm);
2734 if(nbev > 0) nbev = 2;
2740 //__________________________________________________________________________
2741 Int_t AliTRDCalibraFillHisto::TestTracklet( Int_t idet, Int_t row, Int_t iSeed, AliTRDmcm *mcm){
2744 // Build the tracklet and return if the tracklet if finally used or not (1/0)
2749 AliTRDmcmTracklet mcmtracklet = AliTRDmcmTracklet();
2750 //mcmtracklet.Reset();
2751 mcmtracklet.SetDetector(idet);
2752 mcmtracklet.SetRow(row);
2753 mcmtracklet.SetN(0);
2755 Int_t iCol, iCol1, iCol2, track[3];
2756 iCol = mcm->GetSeedCol()[iSeed]; // 0....20 (MCM)
2757 mcm->GetColRange(iCol1,iCol2); // range in the pad plane
2760 for (Int_t iTime = 0; iTime < fTimeMax; iTime++) {
2762 amp[0] = mcm->GetADC(iCol-1,iTime);
2763 amp[1] = mcm->GetADC(iCol ,iTime);
2764 amp[2] = mcm->GetADC(iCol+1,iTime);
2766 if(mcm->IsCluster(iCol,iTime)) {
2768 mcmtracklet.AddCluster(iCol+iCol1,iTime,amp,track);
2771 else if ((iCol+1+1) < 21) {
2773 amp[0] = mcm->GetADC(iCol-1+1,iTime);
2774 amp[1] = mcm->GetADC(iCol +1,iTime);
2775 amp[2] = mcm->GetADC(iCol+1+1,iTime);
2777 if(mcm->IsCluster(iCol+1,iTime)) {
2779 mcmtracklet.AddCluster(iCol+1+iCol1,iTime,amp,track);
2788 nbev = UpdateHistogramcm(&mcmtracklet);
2793 //____________Online trackling in AliTRDtrigger________________________________
2794 Int_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
2797 // Return if the tracklet is finally used or not (1/0) for calibration
2802 //fGoodTracklet = kTRUE;
2804 // Localisation of the Xbins involved
2805 Int_t idect = trk->GetDetector();
2806 Int_t idectrue = trk->GetDetector();
2809 Int_t nbclusters = trk->GetNclusters();
2811 // Eventuelle correction due to track angle in z direction
2812 Float_t correction = 1.0;
2813 if (fMcmCorrectAngle) {
2814 Float_t z = trk->GetRowz();
2815 Float_t r = trk->GetTime0();
2816 correction = r / TMath::Sqrt((r*r+z*z));
2820 Int_t row = trk->GetRow();
2823 // Boucle sur les clusters
2824 // Condition on number of cluster: don't come from the middle of the detector
2827 for(Int_t k =0; k < 36; k++) amph[k]=0.0;
2828 Double_t ampTotal = 0.0;
2830 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
2832 Float_t amp[3] = { 0.0, 0.0, 0.0 };
2833 Int_t time = trk->GetClusterTime(icl);
2834 Int_t col = trk->GetClusterCol(icl);
2836 //CheckGoodTrackletV0(idect,row,col);
2838 amp[0] = trk->GetClusterADC(icl)[0] * correction;
2839 amp[1] = trk->GetClusterADC(icl)[1] * correction;
2840 amp[2] = trk->GetClusterADC(icl)[2] * correction;
2842 ampTotal += (Float_t) (amp[0]+amp[1]+amp[2]);
2843 amph[time]=amp[0]+amp[1]+amp[2];
2845 if(fDebugLevel > 0){
2846 if ( !fDebugStreamer ) {
2848 TDirectory *backup = gDirectory;
2849 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2850 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2853 Double_t amp0 = amp[0];
2854 Double_t amp1 = amp[1];
2855 Double_t amp2 = amp[2];
2857 (* fDebugStreamer) << "UpdateHistogramcm0"<<
2858 "nbclusters="<<nbclusters<<
2865 "detector="<<idectrue<<
2869 } // Boucle clusters
2871 if((amph[0] > 100.0) || (!fGoodTracklet) || (trk->GetNclusters() < fNumberClusters) || (trk->GetNclusters() > fNumberClustersf)) used = 0;
2874 for(Int_t k = 0; k < fTimeMax; k++) UpdateDAQ(idect,0,0,k,amph[k],fTimeMax);
2875 //((TH2I *)GetCH2d()->Fill(ampTotal/30.0,idect));
2879 if(fDebugLevel > 0){
2880 if ( !fDebugStreamer ) {
2882 TDirectory *backup = gDirectory;
2883 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2884 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2887 Double_t amph0 = amph[0];
2888 Double_t amphlast = amph[fTimeMax-1];
2889 Double_t rms = TMath::RMS(fTimeMax,amph);
2890 Int_t goodtracklet = (Int_t) fGoodTracklet;
2892 (* fDebugStreamer) << "UpdateHistogramcm1"<<
2893 "nbclusters="<<nbclusters<<
2894 "ampTotal="<<ampTotal<<
2896 "detector="<<idectrue<<
2898 "amphlast="<<amphlast<<
2899 "goodtracklet="<<goodtracklet<<
2907 //_______________________________________________________________________
2908 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2911 // Look for the maximum by collapsing over the time
2912 // Sum over four pad col and two pad row
2918 Int_t idect = fDetectorPreviousTrack;
2919 //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2921 for(Int_t tb = 0; tb < 36; tb++){
2925 //fGoodTracklet = kTRUE;
2926 //fDetectorPreviousTrack = 0;
2929 ///////////////////////////
2931 /////////////////////////
2935 Double_t integralMax = -1;
2937 for (Int_t ir = 1; ir <= 15; ir++)
2939 for (Int_t ic = 2; ic <= 142; ic++)
2941 Double_t integral = 0;
2942 for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2944 for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2946 if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2947 ic + ishiftC >= 1 && ic + ishiftC <= 144)
2950 for(Int_t tb = 0; tb< fTimeMax; tb++){
2951 integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2957 if (integralMax < integral)
2961 integralMax = integral;
2962 } // check max integral
2966 //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2968 if((imaxRow == 0) || (imaxCol == 0)) {
2972 //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2973 //if(!fGoodTracklet) used = 1;;
2975 // /////////////////////////////////////////////////////
2976 // sum ober 2 row and 4 pad cols for each time bins
2977 // ////////////////////////////////////////////////////
2980 for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2982 for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2984 for(Int_t it = 0; it < fTimeMax; it++){
2985 sum[it] += phvalue[ir][ic][it];
2991 Double_t sumcharge = 0.0;
2992 for(Int_t it = 0; it < fTimeMax; it++){
2993 sumcharge += sum[it];
2994 if(sum[it] > 20.0) nbcl++;
2998 /////////////////////////////////////////////////////////
3000 ////////////////////////////////////////////////////////
3001 if(fDebugLevel > 0){
3002 if ( !fDebugStreamer ) {
3004 TDirectory *backup = gDirectory;
3005 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
3006 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
3009 Double_t amph0 = sum[0];
3010 Double_t amphlast = sum[fTimeMax-1];
3011 Double_t rms = TMath::RMS(fTimeMax,sum);
3012 Int_t goodtracklet = (Int_t) fGoodTracklet;
3013 for(Int_t it = 0; it < fTimeMax; it++){
3014 Double_t clustera = sum[it];
3016 (* fDebugStreamer) << "FillDAQa"<<
3017 "ampTotal="<<sumcharge<<
3020 "detector="<<idect<<
3022 "amphlast="<<amphlast<<
3023 "goodtracklet="<<goodtracklet<<
3024 "clustera="<<clustera<<
3031 ////////////////////////////////////////////////////////
3033 ///////////////////////////////////////////////////////
3034 if(sum[0] > 100.0) used = 1;
3035 if(nbcl < fNumberClusters) used = 1;
3036 if(nbcl > fNumberClustersf) used = 1;
3038 //if(fDetectorPreviousTrack == 15){
3039 // printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
3041 //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
3043 for(Int_t it = 0; it < fTimeMax; it++){
3044 if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
3046 if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
3051 //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
3053 //printf("Pass Detector %d\n",fDetectorPreviousTrack);
3060 //____________Online trackling in AliTRDtrigger________________________________
3061 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
3065 // Fill a simple average pulse height
3069 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
3075 //____________Write_____________________________________________________
3076 //_____________________________________________________________________
3077 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
3080 // Write infos to file
3084 if ( fDebugStreamer ) {
3085 delete fDebugStreamer;
3086 fDebugStreamer = 0x0;
3089 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
3094 ,fNumberUsedPh[1]));
3096 TDirectory *backup = gDirectory;
3102 option = "recreate";
3104 TFile f(filename,option.Data());
3106 TStopwatch stopwatch;
3109 f.WriteTObject(fCalibraVector);
3114 f.WriteTObject(fCH2d);
3119 f.WriteTObject(fPH2d);
3124 f.WriteTObject(fPRF2d);
3127 if(fLinearFitterOn){
3128 AnalyseLinearFitter();
3129 f.WriteTObject(fLinearVdriftFit);
3134 if ( backup ) backup->cd();
3136 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
3137 ,stopwatch.RealTime(),stopwatch.CpuTime()));
3139 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3141 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3142 //___________________________________________probe the histos__________________________________________________
3143 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
3146 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
3147 // debug mode with 2 for TH2I and 3 for TProfile2D
3148 // It gives a pointer to a Double_t[7] with the info following...
3149 // [0] : number of calibration groups with entries
3150 // [1] : minimal number of entries found
3151 // [2] : calibration group number of the min
3152 // [3] : maximal number of entries found
3153 // [4] : calibration group number of the max
3154 // [5] : mean number of entries found
3155 // [6] : mean relative error
3158 Double_t *info = new Double_t[7];
3160 // Number of Xbins (detectors or groups of pads)
3161 Int_t nbins = h->GetNbinsY(); //number of calibration groups
3162 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
3165 Double_t nbwe = 0; //number of calibration groups with entries
3166 Double_t minentries = 0; //minimal number of entries found
3167 Double_t maxentries = 0; //maximal number of entries found
3168 Double_t placemin = 0; //calibration group number of the min
3169 Double_t placemax = -1; //calibration group number of the max
3170 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
3171 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
3173 Double_t counter = 0;
3176 TH1F *nbEntries = 0x0;//distribution of the number of entries
3177 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
3178 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
3180 // Beginning of the loop over the calibration groups
3181 for (Int_t idect = 0; idect < nbins; idect++) {
3183 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
3184 projch->SetDirectory(0);
3186 // Number of entries for this calibration group
3187 Double_t nentries = 0.0;
3189 for (Int_t k = 0; k < nxbins; k++) {
3190 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
3194 for (Int_t k = 0; k < nxbins; k++) {
3195 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3196 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3197 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3205 if((!((Bool_t)nbEntries)) && (nentries > 0)){
3206 nbEntries = new TH1F("Number of entries","Number of entries"
3207 ,100,(Int_t)nentries/2,nentries*2);
3208 nbEntries->SetDirectory(0);
3209 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
3211 nbEntriesPerGroup->SetDirectory(0);
3212 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
3213 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
3214 nbEntriesPerSp->SetDirectory(0);
3217 if(nentries > 0) nbEntries->Fill(nentries);
3218 nbEntriesPerGroup->Fill(idect+0.5,nentries);
3219 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3224 if(nentries > maxentries){
3225 maxentries = nentries;
3229 minentries = nentries;
3231 if(nentries < minentries){
3232 minentries = nentries;
3238 meanstats += nentries;
3240 }//calibration groups loop
3242 if(nbwe > 0) meanstats /= nbwe;
3243 if(counter > 0) meanrelativerror /= counter;
3245 AliInfo(Form("There are %f calibration groups with entries",nbwe));
3246 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3247 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3248 AliInfo(Form("The mean number of entries is %f",meanstats));
3249 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3252 info[1] = minentries;
3254 info[3] = maxentries;
3256 info[5] = meanstats;
3257 info[6] = meanrelativerror;
3260 gStyle->SetPalette(1);
3261 gStyle->SetOptStat(1111);
3262 gStyle->SetPadBorderMode(0);
3263 gStyle->SetCanvasColor(10);
3264 gStyle->SetPadLeftMargin(0.13);
3265 gStyle->SetPadRightMargin(0.01);
3266 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3269 nbEntries->Draw("");
3271 nbEntriesPerSp->SetStats(0);
3272 nbEntriesPerSp->Draw("");
3273 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3275 nbEntriesPerGroup->SetStats(0);
3276 nbEntriesPerGroup->Draw("");
3282 //____________________________________________________________________________
3283 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3286 // Return a Int_t[4] with:
3287 // 0 Mean number of entries
3288 // 1 median of number of entries
3289 // 2 rms of number of entries
3290 // 3 number of group with entries
3293 Double_t *stat = new Double_t[4];
3296 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
3297 Double_t *weight = new Double_t[nbofgroups];
3298 Int_t *nonul = new Int_t[nbofgroups];
3300 for(Int_t k = 0; k < nbofgroups; k++){
3301 if(fEntriesCH[k] > 0) {
3303 nonul[(Int_t)stat[3]] = fEntriesCH[k];
3306 else weight[k] = 0.0;
3308 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
3309 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
3310 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3315 //____________________________________________________________________________
3316 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3319 // Return a Int_t[4] with:
3320 // 0 Mean number of entries
3321 // 1 median of number of entries
3322 // 2 rms of number of entries
3323 // 3 number of group with entries
3326 Double_t *stat = new Double_t[4];
3329 Int_t nbofgroups = 540;
3330 Double_t *weight = new Double_t[nbofgroups];
3331 Int_t *nonul = new Int_t[nbofgroups];
3333 for(Int_t k = 0; k < nbofgroups; k++){
3334 if(fEntriesLinearFitter[k] > 0) {
3336 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3339 else weight[k] = 0.0;
3341 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3342 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3343 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
3348 //////////////////////////////////////////////////////////////////////////////////////
3350 //////////////////////////////////////////////////////////////////////////////////////
3351 //_____________________________________________________________________________
3352 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3355 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3356 // If fNgroupprf is zero then no binning in tnp
3360 name += fCalibraMode->GetNz(2);
3362 name += fCalibraMode->GetNrphi(2);
3366 if(fNgroupprf != 0){
3368 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3369 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3370 fPRF2d->SetYTitle("Det/pad groups");
3371 fPRF2d->SetXTitle("Position x/W [pad width units]");
3372 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3373 fPRF2d->SetStats(0);
3376 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3377 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3378 fPRF2d->SetYTitle("Det/pad groups");
3379 fPRF2d->SetXTitle("Position x/W [pad width units]");
3380 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3381 fPRF2d->SetStats(0);
3386 //_____________________________________________________________________________
3387 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3390 // Create the 2D histos
3394 name += fCalibraMode->GetNz(1);
3396 name += fCalibraMode->GetNrphi(1);
3398 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3399 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3401 fPH2d->SetYTitle("Det/pad groups");
3402 fPH2d->SetXTitle("time [#mus]");
3403 fPH2d->SetZTitle("<PH> [a.u.]");
3407 //_____________________________________________________________________________
3408 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3411 // Create the 2D histos
3415 name += fCalibraMode->GetNz(0);
3417 name += fCalibraMode->GetNrphi(0);
3419 fCH2d = new TH2I("CH2d",(const Char_t *) name
3420 ,fNumberBinCharge,0,300,nn,0,nn);
3421 fCH2d->SetYTitle("Det/pad groups");
3422 fCH2d->SetXTitle("charge deposit [a.u]");
3423 fCH2d->SetZTitle("counts");
3428 //////////////////////////////////////////////////////////////////////////////////
3429 // Set relative scale
3430 /////////////////////////////////////////////////////////////////////////////////
3431 //_____________________________________________________________________________
3432 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3435 // Set the factor that will divide the deposited charge
3436 // to fit in the histo range [0,300]
3439 if (RelativeScale > 0.0) {
3440 fRelativeScale = RelativeScale;
3443 AliInfo("RelativeScale must be strict positif!");
3447 //////////////////////////////////////////////////////////////////////////////////
3448 // Quick way to fill a histo
3449 //////////////////////////////////////////////////////////////////////////////////
3450 //_____________________________________________________________________
3451 void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3454 // FillCH2d: Marian style
3457 //skip simply the value out of range
3458 if((y>=300.0) || (y<0.0)) return;
3460 //Calcul the y place
3461 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3462 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3465 fCH2d->GetArray()[place]++;
3469 //////////////////////////////////////////////////////////////////////////////////
3470 // Geometrical functions
3471 ///////////////////////////////////////////////////////////////////////////////////
3472 //_____________________________________________________________________________
3473 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3476 // Reconstruct the layer number from the detector number
3479 return ((Int_t) (d % 6));
3483 //_____________________________________________________________________________
3484 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3487 // Reconstruct the stack number from the detector number
3489 const Int_t kNlayer = 6;
3491 return ((Int_t) (d % 30) / kNlayer);
3495 //_____________________________________________________________________________
3496 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3499 // Reconstruct the sector number from the detector number
3503 return ((Int_t) (d / fg));
3506 ///////////////////////////////////////////////////////////////////////////////////
3507 // Getter functions for DAQ of the CH2d and the PH2d
3508 //////////////////////////////////////////////////////////////////////////////////
3509 //_____________________________________________________________________
3510 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3513 // return pointer to fPH2d TProfile2D
3514 // create a new TProfile2D if it doesn't exist allready
3520 fTimeMax = nbtimebin;
3521 fSf = samplefrequency;
3527 //_____________________________________________________________________
3528 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3531 // return pointer to fCH2d TH2I
3532 // create a new TH2I if it doesn't exist allready
3541 ////////////////////////////////////////////////////////////////////////////////////////////
3542 // Drift velocity calibration
3543 ///////////////////////////////////////////////////////////////////////////////////////////
3544 //_____________________________________________________________________
3545 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3548 // return pointer to TLinearFitter Calibration
3549 // if force is true create a new TLinearFitter if it doesn't exist allready
3552 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3553 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3556 // if we are forced and TLinearFitter doesn't yet exist create it
3558 // new TLinearFitter
3559 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3560 fLinearFitterArray.AddAt(linearfitter,detector);
3561 return linearfitter;
3564 //____________________________________________________________________________
3565 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3568 // Analyse array of linear fitter because can not be written
3569 // Store two arrays: one with the param the other one with the error param + number of entries
3572 for(Int_t k = 0; k < 540; k++){
3573 TLinearFitter *linearfitter = GetLinearFitter(k);
3574 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3575 TVectorD *par = new TVectorD(2);
3576 TVectorD pare = TVectorD(2);
3577 TVectorD *parE = new TVectorD(3);
3578 linearfitter->Eval();
3579 linearfitter->GetParameters(*par);
3580 linearfitter->GetErrors(pare);
3581 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3582 (*parE)[0] = pare[0]*ppointError;
3583 (*parE)[1] = pare[1]*ppointError;
3584 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3585 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3586 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);