]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDCalibraFillHisto.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Remove obsolete dedx root-file
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraFillHisto.cxx
CommitLineData
55a288e5 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/* $Id$ */
17
18/////////////////////////////////////////////////////////////////////////////////
170c35f1 19//
55a288e5 20// AliTRDCalibraFillHisto
21//
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.
170c35f1 25// The user has to choose with the functions SetNz and SetNrphi the precision of the calibration (see AliTRDCalibraMode).
55a288e5 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.
30//
31// Author:
32// R. Bailhache (R.Bailhache@gsi.de)
33//
34//////////////////////////////////////////////////////////////////////////////////////
35
36#include <TTree.h>
37#include <TProfile2D.h>
38#include <TProfile.h>
39#include <TFile.h>
40#include <TChain.h>
41#include <TStyle.h>
42#include <TCanvas.h>
43#include <TGraphErrors.h>
44#include <TObjArray.h>
45#include <TH1F.h>
46#include <TH2I.h>
47#include <TH2.h>
48#include <TStopwatch.h>
49#include <TMath.h>
50#include <TDirectory.h>
51#include <TROOT.h>
170c35f1 52#include <TTreeStream.h>
53#include <TVectorD.h>
55a288e5 54
55#include "AliLog.h"
56#include "AliCDBManager.h"
57
58#include "AliTRDCalibraFillHisto.h"
59#include "AliTRDCalibraMode.h"
60#include "AliTRDCalibraVector.h"
61#include "AliTRDcalibDB.h"
62#include "AliTRDCommonParam.h"
63#include "AliTRDmcmTracklet.h"
64#include "AliTRDpadPlane.h"
65#include "AliTRDcluster.h"
66#include "AliTRDtrack.h"
170c35f1 67#include "AliTRDRawStream.h"
68#include "AliRawReader.h"
69#include "AliRawReaderDate.h"
70
71#ifdef ALI_DATE
72#include "event.h"
73#endif
55a288e5 74
75
76ClassImp(AliTRDCalibraFillHisto)
77
78AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
79Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
80
81//_____________singleton implementation_________________________________________________
82AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
83{
84 //
85 // Singleton implementation
86 //
87
88 if (fgTerminated != kFALSE) {
89 return 0;
90 }
91
92 if (fgInstance == 0) {
93 fgInstance = new AliTRDCalibraFillHisto();
94 }
95
96 return fgInstance;
97
98}
99
100//______________________________________________________________________________________
101void AliTRDCalibraFillHisto::Terminate()
102{
103 //
104 // Singleton implementation
105 // Deletes the instance of this class
106 //
107
108 fgTerminated = kTRUE;
109
110 if (fgInstance != 0) {
111 delete fgInstance;
112 fgInstance = 0;
113 }
114
115}
55a288e5 116//______________________________________________________________________________________
117AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
118 :TObject()
119 ,fMITracking(kFALSE)
120 ,fMcmTracking(kFALSE)
121 ,fMcmCorrectAngle(kFALSE)
122 ,fCH2dOn(kFALSE)
123 ,fPH2dOn(kFALSE)
124 ,fPRF2dOn(kFALSE)
125 ,fHisto2d(kFALSE)
126 ,fVector2d(kFALSE)
170c35f1 127 ,fLinearFitterOn(kFALSE)
128 ,fLinearFitterDebugOn(kFALSE)
55a288e5 129 ,fRelativeScale(0)
170c35f1 130 ,fThresholdClusterPRF2(15.0)
131 ,fCalibraMode(new AliTRDCalibraMode())
132 ,fDebugStreamer(0)
133 ,fDebugLevel(0)
55a288e5 134 ,fDetectorAliTRDtrack(kFALSE)
55a288e5 135 ,fDetectorPreviousTrack(-1)
170c35f1 136 ,fNumberClusters(18)
137 ,fProcent(6.0)
138 ,fDifference(17)
139 ,fNumberTrack(0)
140 ,fTimeMax(0)
141 ,fSf(10.0)
142 ,fNumberBinCharge(100)
143 ,fNumberBinPRF(40)
144 ,fNgroupprf(0)
145 ,fListClusters(new TObjArray())
146 ,fPar0(0x0)
147 ,fPar1(0x0)
148 ,fPar2(0x0)
149 ,fPar3(0x0)
150 ,fPar4(0x0)
55a288e5 151 ,fAmpTotal(0x0)
152 ,fPHPlace(0x0)
153 ,fPHValue(0x0)
170c35f1 154 ,fGoodTracklet(kTRUE)
155 ,fGoodTrack(kTRUE)
156 ,fEntriesCH(0x0)
157 ,fEntriesLinearFitter(0x0)
158 ,fCalibraVector(0x0)
55a288e5 159 ,fPH2d(0x0)
160 ,fPRF2d(0x0)
161 ,fCH2d(0x0)
170c35f1 162 ,fLinearFitterArray(0)
163 ,fLinearFitterHistoArray(0)
55a288e5 164{
165 //
166 // Default constructor
167 //
168
170c35f1 169 //
170 // Init some default values
171 //
55a288e5 172
170c35f1 173 fNumberUsedCh[0] = 0;
174 fNumberUsedCh[1] = 0;
175 fNumberUsedPh[0] = 0;
176 fNumberUsedPh[1] = 0;
177
55a288e5 178}
55a288e5 179//______________________________________________________________________________________
180AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
170c35f1 181:TObject(c)
182,fMITracking(c.fMITracking)
183,fMcmTracking(c.fMcmTracking)
184,fMcmCorrectAngle(c.fMcmCorrectAngle)
185,fCH2dOn(c.fCH2dOn)
186,fPH2dOn(c.fPH2dOn)
187,fPRF2dOn(c.fPRF2dOn)
188,fHisto2d(c.fHisto2d)
189,fVector2d(c.fVector2d)
190,fLinearFitterOn(c.fLinearFitterOn)
191,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
192,fRelativeScale(c.fRelativeScale)
193,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
194,fCalibraMode(0x0)
195,fDebugStreamer(0)
196,fDebugLevel(c.fDebugLevel)
197,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
198,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
199,fNumberClusters(c.fNumberClusters)
200,fProcent(c.fProcent)
201,fDifference(c.fDifference)
202,fNumberTrack(c.fNumberTrack)
203,fTimeMax(c.fTimeMax)
204,fSf(c.fSf)
205,fNumberBinCharge(c.fNumberBinCharge)
206,fNumberBinPRF(c.fNumberBinPRF)
207,fNgroupprf(c.fNgroupprf)
208,fListClusters(new TObjArray())
209,fPar0(c.fPar0)
210,fPar1(c.fPar1)
211,fPar2(c.fPar2)
212,fPar3(c.fPar3)
213,fPar4(c.fPar4)
214,fAmpTotal(c.fAmpTotal)
215,fPHPlace(c.fPHPlace)
216,fPHValue(c.fPHValue)
217,fGoodTracklet(c.fGoodTracklet)
218,fGoodTrack(c.fGoodTrack)
219,fEntriesCH(c.fEntriesCH)
220,fEntriesLinearFitter(fEntriesLinearFitter)
221,fCalibraVector(0x0)
222,fPH2d(0x0)
223,fPRF2d(0x0)
224,fCH2d(0x0)
225,fLinearFitterArray(0)
226,fLinearFitterHistoArray(0)
55a288e5 227{
228 //
229 // Copy constructor
230 //
170c35f1 231 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
232 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
233 if(c.fPH2d) {
234 fPH2d = new TProfile2D(*c.fPH2d);
235 fPH2d->SetDirectory(0);
236 }
237 if(c.fPRF2d) {
238 fPRF2d = new TProfile2D(*c.fPRF2d);
239 fPRF2d->SetDirectory(0);
240 }
241 if(c.fCH2d) {
242 fCH2d = new TH2I(*c.fCH2d);
243 fCH2d->SetDirectory(0);
244 }
245 if(c.fLinearFitterOn){
246 fLinearFitterArray.Expand(540);
247 for (Int_t cb = 0; cb < 540; cb++){
248 const TLinearFitter *linearFitter = (TLinearFitter*)c.fLinearFitterArray.UncheckedAt(cb);
249 if ( linearFitter != 0x0 ) fLinearFitterArray.AddAt(new TLinearFitter(*linearFitter), cb);
250 }
251 }
252 if(c.fLinearFitterDebugOn){
253 fLinearFitterHistoArray.Expand(540);
254 for (Int_t cb = 0; cb < 540; cb++){
255 const TH2F *linearfitterhisto = (TH2F*)c.fLinearFitterHistoArray.UncheckedAt(cb);
256 if ( linearfitterhisto != 0x0 ){
257 TH2F *hNew = new TH2F(*linearfitterhisto);
258 hNew->SetDirectory(0);
259 fLinearFitterHistoArray.AddAt(hNew,cb);
260 }
261 }
262 }
55a288e5 263}
55a288e5 264//____________________________________________________________________________________
265AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
266{
267 //
268 // AliTRDCalibraFillHisto destructor
269 //
270
271 ClearHistos();
170c35f1 272 if ( fDebugStreamer ) delete fDebugStreamer;
55a288e5 273
274}
275
276//_____________________________________________________________________________
277void AliTRDCalibraFillHisto::Destroy()
278{
279 //
280 // Delete instance
281 //
282
283 if (fgInstance) {
284 delete fgInstance;
285 fgInstance = 0x0;
286 }
287
288}
289
290//_____________________________________________________________________________
291void AliTRDCalibraFillHisto::ClearHistos()
292{
293 //
294 // Delete the histos
295 //
296
297 if (fPH2d) {
298 delete fPH2d;
299 fPH2d = 0x0;
300 }
301 if (fCH2d) {
302 delete fCH2d;
303 fCH2d = 0x0;
304 }
305 if (fPRF2d) {
306 delete fPRF2d;
307 fPRF2d = 0x0;
308 }
309
310}
55a288e5 311//____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
312Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
313{
314 //
315 // For the offline tracking
316 // This function will be called in the function AliReconstruction::Run()
317 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
318 //
319
320 // DB Setting
321 // Get cal
322 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
323 if (!cal) {
324 AliInfo("Could not get calibDB");
325 return kFALSE;
326 }
327 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
328 if (!parCom) {
329 AliInfo("Could not get CommonParam");
330 return kFALSE;
331 }
332
333 // Some parameters
334 fTimeMax = cal->GetNumberOfTimeBins();
335 fSf = parCom->GetSamplingFrequency();
170c35f1 336 fRelativeScale = 20;
337
338 //Init the tracklet parameters
339 fPar0 = new Double_t[fTimeMax];
340 fPar1 = new Double_t[fTimeMax];
341 fPar2 = new Double_t[fTimeMax];
342 fPar3 = new Double_t[fTimeMax];
343 fPar4 = new Double_t[fTimeMax];
344
345 for(Int_t k = 0; k < fTimeMax; k++){
346 fPar0[k] = 0.0;
347 fPar1[k] = 0.0;
348 fPar2[k] = 0.0;
349 fPar3[k] = 0.0;
350 fPar4[k] = 0.0;
55a288e5 351 }
170c35f1 352
55a288e5 353 //If vector method On initialised all the stuff
354 if(fVector2d){
355 fCalibraVector = new AliTRDCalibraVector();
170c35f1 356 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
357 fCalibraVector->SetTimeMax(fTimeMax);
358 if(fNgroupprf != 0) {
359 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
360 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
361 }
362 else {
363 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
364 fCalibraVector->SetPRFRange(1.5);
365 }
55a288e5 366 }
170c35f1 367
55a288e5 368 // Create the 2D histos corresponding to the pad groupCalibration mode
369 if (fCH2dOn) {
370
371 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
372 ,fCalibraMode->GetNz(0)
373 ,fCalibraMode->GetNrphi(0)));
374
375 // Calcul the number of Xbins
170c35f1 376 Int_t Ntotal0 = CalculateTotalNumberOfBins(0);
55a288e5 377 // Create the 2D histo
378 if (fHisto2d) {
379 CreateCH2d(Ntotal0);
380 }
55a288e5 381 // Variable
382 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
383 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
384 fAmpTotal[k] = 0.0;
385 }
170c35f1 386 //Statistics
387 fEntriesCH = new Int_t[Ntotal0];
388 for(Int_t k = 0; k < Ntotal0; k++){
389 fEntriesCH[k] = 0;
390 }
55a288e5 391 }
55a288e5 392 if (fPH2dOn) {
393
394 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
395 ,fCalibraMode->GetNz(1)
396 ,fCalibraMode->GetNrphi(1)));
397
398 // Calcul the number of Xbins
170c35f1 399 Int_t Ntotal1 = CalculateTotalNumberOfBins(1);
55a288e5 400 // Create the 2D histo
401 if (fHisto2d) {
402 CreatePH2d(Ntotal1);
403 }
55a288e5 404 // Variable
405 fPHPlace = new Short_t[fTimeMax];
406 for (Int_t k = 0; k < fTimeMax; k++) {
407 fPHPlace[k] = -1;
408 }
409 fPHValue = new Float_t[fTimeMax];
410 for (Int_t k = 0; k < fTimeMax; k++) {
411 fPHValue[k] = 0.0;
412 }
170c35f1 413 }
414 if (fLinearFitterOn) {
415 fLinearFitterArray.Expand(540);
416 fLinearFitterArray.SetName("ArrayLinearFitters");
417 fEntriesLinearFitter = new Int_t[540];
418 for(Int_t k = 0; k < 540; k++){
419 fEntriesLinearFitter[k] = 0;
420 }
421 if(fLinearFitterDebugOn) {
422 fLinearFitterHistoArray.Expand(540);
423 fLinearFitterHistoArray.SetName("ArrayHistos");
424 }
55a288e5 425 }
426
427 if (fPRF2dOn) {
428
429 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
430 ,fCalibraMode->GetNz(2)
431 ,fCalibraMode->GetNrphi(2)));
432
433 // Calcul the number of Xbins
170c35f1 434 Int_t Ntotal2 = CalculateTotalNumberOfBins(2);
55a288e5 435 // Create the 2D histo
436 if (fHisto2d) {
437 CreatePRF2d(Ntotal2);
438 }
55a288e5 439 }
440
441 return kTRUE;
442
443}
444
445//____________Functions for filling the histos in the code_____________________
446
447//____________Offine tracking in the AliTRDtracker_____________________________
448Bool_t AliTRDCalibraFillHisto::ResetTrack()
449{
450 //
451 // For the offline tracking
452 // This function will be called in the function
453 // AliTRDtracker::FollowBackPropagation() at the beginning
454 // Reset the parameter to know we have a new TRD track
455 //
456
457 fDetectorAliTRDtrack = kFALSE;
170c35f1 458 fGoodTrack = kTRUE;
55a288e5 459 return kTRUE;
460
461}
170c35f1 462//____________Offine tracking in the AliTRDtracker_____________________________
463void AliTRDCalibraFillHisto::ResetfVariables()
464{
465 //
466 // Reset values per tracklet
467 //
468
469 // Reset the list of clusters
470 fListClusters->Clear();
471
472 //Reset the tracklet parameters
473 for(Int_t k = 0; k < fTimeMax; k++){
474 fPar0[k] = 0.0;
475 fPar1[k] = 0.0;
476 fPar2[k] = 0.0;
477 fPar3[k] = 0.0;
478 fPar4[k] = 0.0;
479 }
480
481
482 //Reset good tracklet
483 fGoodTracklet = kTRUE;
55a288e5 484
170c35f1 485 // Reset the fPHValue
486 if (fPH2dOn) {
487 //Reset the fPHValue and fPHPlace
488 for (Int_t k = 0; k < fTimeMax; k++) {
489 fPHValue[k] = 0.0;
490 fPHPlace[k] = -1;
491 }
492 }
493
494 // Reset the fAmpTotal where we put value
495 if (fCH2dOn) {
496 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
497 fAmpTotal[k] = 0.0;
498 }
499 }
500
501}
55a288e5 502//____________Offline tracking in the AliTRDtracker____________________________
503Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
504{
505 //
506 // For the offline tracking
507 // This function will be called in the function
508 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
509 // of TRD tracks
510 // Fill the 2D histos or the vectors with the info of the clusters at
511 // the end of a detectors if the track is "good"
512 //
513
55a288e5 514 // Localisation of the detector
515 Int_t detector = cl->GetDetector();
55a288e5 516
517 // Fill the infos for the previous clusters if not the same
518 // detector anymore or if not the same track
519 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
520 (fDetectorPreviousTrack != -1)) {
521
522 fNumberTrack++;
523
524 // If the same track, then look if the previous detector is in
525 // the same plane, if yes: not a good track
526 if (fDetectorAliTRDtrack &&
527 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
528 fGoodTrack = kFALSE;
529 }
530
531 // Fill only if the track doesn't touch a masked pad or doesn't
532 // appear in the middle (fGoodTrack)
170c35f1 533 if (fGoodTrack && fGoodTracklet) {
55a288e5 534
170c35f1 535 // drift velocity unables to cut bad tracklets
536 Bool_t pass = FindP1TrackPH();
537
55a288e5 538 // Gain calibration
539 if (fCH2dOn) {
540 FillTheInfoOfTheTrackCH();
541 }
542
543 // PH calibration
544 if (fPH2dOn) {
545 FillTheInfoOfTheTrackPH();
546 }
170c35f1 547
548 if(pass && fPRF2dOn) HandlePRF();
549
55a288e5 550
551 } // if a good track
552
553 ResetfVariables();
554
555 } // Fill at the end the charge
556
557 // Calcul the position of the detector
558 if (detector != fDetectorPreviousTrack) {
559 LocalisationDetectorXbins(detector);
560 }
170c35f1 561
55a288e5 562 // Reset the detector
563 fDetectorPreviousTrack = detector;
564 fDetectorAliTRDtrack = kTRUE;
55a288e5 565
170c35f1 566 // Store the infos of the tracklets
567 AliTRDcluster *kcl = new AliTRDcluster(*cl);
568 fListClusters->Add((TObject *)kcl);
569 Int_t time = cl->GetLocalTimeBin();
570 fPar0[time] = t->GetY();
571 fPar1[time] = t->GetZ();
572 fPar2[time] = t->GetSnp();
573 fPar3[time] = t->GetTgl();
574 fPar4[time] = t->Get1Pt();
575
576 // Store the info bis of the tracklet
577 Int_t *rowcol = CalculateRowCol(cl);
578 CheckGoodTracklet(detector,rowcol);
579 Int_t group[2] = {0,0};
580 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
581 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
582 StoreInfoCHPH(cl,t,group);
583
55a288e5 584 return kTRUE;
585
586}
55a288e5 587//____________Online trackling in AliTRDtrigger________________________________
588Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
589{
590 //
591 // For the tracking
592 // This function will be called in the function AliTRDtrigger::TestTracklet
593 // before applying the pt cut on the tracklets
594 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
595 //
596
597 // Localisation of the Xbins involved
598 Int_t idect = trk->GetDetector();
599 LocalisationDetectorXbins(idect);
600
601 // Get the parameter object
602 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
603 if (!cal) {
604 AliInfo("Could not get calibDB");
605 return kFALSE;
606 }
607
608 // Reset
609 ResetfVariables();
610
611 // Row of the tracklet and position in the pad groups
170c35f1 612 Int_t *rowcol = new Int_t[2];
613 rowcol[0] = trk->GetRow();
614 Int_t group[3] = {-1,-1,-1};
615
55a288e5 616 // Eventuelle correction due to track angle in z direction
617 Float_t correction = 1.0;
618 if (fMcmCorrectAngle) {
619 Float_t z = trk->GetRowz();
620 Float_t r = trk->GetTime0();
621 correction = r / TMath::Sqrt((r*r+z*z));
622 }
623
624 // Boucle sur les clusters
625 // Condition on number of cluster: don't come from the middle of the detector
626 if (trk->GetNclusters() >= fNumberClusters) {
627
628 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
629
630 Float_t amp[3] = { 0.0, 0.0, 0.0 };
631 Int_t time = trk->GetClusterTime(icl);
170c35f1 632 rowcol[1] = trk->GetClusterCol(icl);
55a288e5 633
634 amp[0] = trk->GetClusterADC(icl)[0] * correction;
635 amp[1] = trk->GetClusterADC(icl)[1] * correction;
636 amp[2] = trk->GetClusterADC(icl)[2] * correction;
637
638
639 if ((amp[0] < 0.0) ||
640 (amp[1] < 0.0) ||
641 (amp[2] < 0.0)) {
642 continue;
643 }
644
645 // Col of cluster and position in the pad groups
170c35f1 646 if(fCH2dOn) {
647 group[0] = CalculateCalibrationGroup(0,rowcol);
648 fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
55a288e5 649 }
170c35f1 650 if(fPH2dOn) {
651 group[1] = CalculateCalibrationGroup(1,rowcol);
652 fPHPlace[time] = group[1];
55a288e5 653 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
654 }
170c35f1 655 if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
656
657 // See if we are not near a masked pad fGoodTracklet
658 CheckGoodTracklet(idect,rowcol);
659
660 // Fill PRF direct without tnp bins...only for monitoring...
661 if (fPRF2dOn && fGoodTracklet) {
55a288e5 662
663 if ((amp[0] > fThresholdClusterPRF2) &&
664 (amp[1] > fThresholdClusterPRF2) &&
665 (amp[2] > fThresholdClusterPRF2) &&
666 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
667
668 // Security of the denomiateur is 0
669 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
670 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
671 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
672 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
673 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
674
675 if (TMath::Abs(xcenter) < 0.5) {
676 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
677 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
678 // Fill only if it is in the drift region!
679 if (((Float_t) time / fSf) > 0.3) {
680 if (fHisto2d) {
170c35f1 681 fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
682 fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
683 fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
55a288e5 684 }
685 if (fVector2d) {
170c35f1 686 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],xcenter,ycenter);
687 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],-(xcenter+1.0),yminus);
688 fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],(1.0-xcenter),ymax);
55a288e5 689 }
690 }//in the drift region
691 }//in the middle
692 }//denominateur security
693 }//cluster shape and thresholds
694 }//good and PRF On
695
696 } // Boucle clusters
697
698 // Fill the charge
170c35f1 699 if(fGoodTracklet){
700 if (fCH2dOn) FillTheInfoOfTheTrackCH();
701 if (fPH2dOn) FillTheInfoOfTheTrackPH();
55a288e5 702 }
703
704 fNumberTrack++;
705
706 } // Condition on number of clusters
707
708 return kTRUE;
709
710}
55a288e5 711//_____________________________________________________________________________
170c35f1 712Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
55a288e5 713{
714 //
170c35f1 715 // Calculate the row and col number of the cluster
55a288e5 716 //
717
170c35f1 718
719 Int_t *rowcol = new Int_t[2];
720 rowcol[0] = 0;
721 rowcol[1] = 0;
722
55a288e5 723 // Get the parameter object
170c35f1 724 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
725 if (!parCom) {
726 AliInfo("Could not get CommonParam");
727 return rowcol;
55a288e5 728 }
170c35f1 729
730 // Localisation of the detector
731 Int_t detector = cl->GetDetector();
732 Int_t chamber = GetChamber(detector);
733 Int_t plane = GetPlane(detector);
734
735 // Localisation of the cluster
736 Double_t pos[3] = { 0.0, 0.0, 0.0 };
737 pos[0] = ((AliCluster *)cl)->GetX();
738 pos[1] = cl->GetY();
739 pos[2] = cl->GetZ();
740
741 // Position of the cluster
742 AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
743 Int_t row = padplane->GetPadRowNumber(pos[2]);
744 //Do not take from here because it was corrected from ExB already....
745 //Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
746 //Double_t offsettilt = padplane->GetTiltOffset(offsetz);
747 //Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
748 //Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt);
749 Int_t col = cl->GetPad();
750
751 //return
752 rowcol[0] = row;
753 rowcol[1] = col;
754 return rowcol;
755
756}
757//_____________________________________________________________________________
758void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
759{
760 //
761 // See if we are not near a masked pad
762 //
763
764 Int_t row = rowcol[0];
765 Int_t col = rowcol[1];
766
767 if (!IsPadOn(detector, col, row)) {
768 fGoodTracklet = kFALSE;
769 }
770
771 if (col > 0) {
772 if (!IsPadOn(detector, col-1, row)) {
773 fGoodTracklet = kFALSE;
774 }
775 }
776
777 if (col < 143) {
778 if (!IsPadOn(detector, col+1, row)) {
779 fGoodTracklet = kFALSE;
780 }
781 }
782
783}
784//_____________________________________________________________________________
785Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
786{
787 //
788 // Look in the choosen database if the pad is On.
789 // If no the track will be "not good"
790 //
791
792 // Get the parameter object
793 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
794 if (!cal) {
795 AliInfo("Could not get calibDB");
796 return kFALSE;
797 }
798
799 if (!cal->IsChamberInstalled(detector) ||
55a288e5 800 cal->IsChamberMasked(detector) ||
801 cal->IsPadMasked(detector,col,row)) {
802 return kFALSE;
803 }
804 else {
805 return kTRUE;
806 }
807
808}
170c35f1 809//_____________________________________________________________________________
810Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
811{
812 //
813 // Calculate the calibration group number for i
814 //
815
816 // Row of the cluster and position in the pad groups
817 Int_t posr = 0;
818 if (fCalibraMode->GetNnZ(i) != 0) {
819 posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
820 }
821
822
823 // Col of the cluster and position in the pad groups
824 Int_t posc = 0;
825 if (fCalibraMode->GetNnRphi(i) != 0) {
826 posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
827 }
828
829 return posc*fCalibraMode->GetNfragZ(i)+posr;
830
831}
832//____________________________________________________________________________________
833Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
834{
835 //
836 // Calculate the total number of calibration groups
837 //
838
839 Int_t Ntotal = 0;
840 fCalibraMode->ModePadCalibration(2,i);
841 fCalibraMode->ModePadFragmentation(0,2,0,i);
842 fCalibraMode->SetDetChamb2(i);
843 Ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
844 fCalibraMode->ModePadCalibration(0,i);
845 fCalibraMode->ModePadFragmentation(0,0,0,i);
846 fCalibraMode->SetDetChamb0(i);
847 Ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
848 AliInfo(Form("Total number of Xbins: %d",Ntotal));
849 return Ntotal;
850
851}
852//____________Set the pad calibration variables for the detector_______________
853Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
854{
855 //
856 // For the detector calcul the first Xbins and set the number of row
857 // and col pads per calibration groups, the number of calibration
858 // groups in the detector.
859 //
860
861 // first Xbins of the detector
862 if (fCH2dOn) {
863 fCalibraMode->CalculXBins(detector,0);
864 }
865 if (fPH2dOn) {
866 fCalibraMode->CalculXBins(detector,1);
867 }
868 if (fPRF2dOn) {
869 fCalibraMode->CalculXBins(detector,2);
870 }
871
872 // fragmentation of idect
873 for (Int_t i = 0; i < 3; i++) {
874 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
875 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
876 , (Int_t) GetChamber(detector)
877 , (Int_t) GetSector(detector),i);
878 }
879
880 return kTRUE;
55a288e5 881
170c35f1 882}
883//_____________________________________________________________________________
884void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group)
885{
886 //
887 // Store the infos in fAmpTotal, fPHPlace and fPHValue
888 //
889
890 // Charge in the cluster
891 Float_t q = TMath::Abs(cl->GetQ());
892 Int_t time = cl->GetLocalTimeBin();
893
894 // Correction due to the track angle
895 Float_t correction = 1.0;
896 Float_t normalisation = 6.67;
897 if ((q >0) && (t->GetNdedx() > 0)) {
898 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
899 }
900
901 // Fill the fAmpTotal with the charge
902 if (fCH2dOn) {
903 fAmpTotal[(Int_t) group[0]] += correction;
904 }
905
906 // Fill the fPHPlace and value
907 if (fPH2dOn) {
908 fPHPlace[time] = group[1];
909 fPHValue[time] = correction;
910 }
911
912}
913//_____________________________________________________________________
914Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDRawStream *rawStream)
915{
916 //
917 // Event Processing loop - AliTRDRawStream
918 //
919
920 Bool_t withInput = kFALSE;
921
922 Int_t phvalue[36];
923 //Int_t row[36];
924 //Int_t col[36];
925 for(Int_t k = 0; k < 36; k++){
926 phvalue[k] = 0;
927 //row[k] = -1;
928 //col[36] = -1;
929 }
930 fDetectorPreviousTrack = -1;
931 Int_t nbtimebin = 0;
932
933 while (rawStream->Next()) {
934
935 Int_t idetector = rawStream->GetDet(); // current detector
936 if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
937 if(TMath::Mean(nbtimebin,phvalue)>0.0){
938 withInput = kTRUE;
939 for(Int_t k = 0; k < nbtimebin; k++){
940 UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
941 phvalue[k] = 0;
942 //row[k] = -1;
943 //col[k] = -1;
944 }
945 }
946 }
947 fDetectorPreviousTrack = idetector;
948 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
949 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
950 //row[iTimeBin] = rawStream->GetRow(); // current row
951 //col[iTimeBin] = rawStream->GetCol(); // current col
952 Int_t *signal = rawStream->GetSignals(); // current ADC signal
953
954 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
955 Int_t n = 0;
956 for(Int_t itime = iTimeBin; itime < fin; itime++){
957 // should extract baseline here!
958 if(signal[n]>5.0) phvalue[itime] = signal[n];
959 n++;
960 }
961 }
962
963 // fill the last one
964 if(fDetectorPreviousTrack != -1){
965 if(TMath::Mean(nbtimebin,phvalue)>0.0){
966 withInput = kTRUE;
967 for(Int_t k = 0; k < nbtimebin; k++){
968 UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
969 phvalue[k] = 0;
970 //row[k] = -1;
971 //col[k] = -1;
972 }
973 }
974 }
975
976 return withInput;
977
978}
979//_____________________________________________________________________
980Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
981{
982 //
983 // Event processing loop - AliRawReader
984 //
985
986
987 AliTRDRawStream rawStream(rawReader);
988
989 rawReader->Select("TRD");
990
991 return ProcessEventDAQ(&rawStream);
992}
993//_________________________________________________________________________
994Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(
995#ifdef ALI_DATE
996 eventHeaderStruct *event
997#else
998 eventHeaderStruct* /*event*/
999
1000#endif
1001 )
1002{
1003 //
1004 // process date event
1005 //
1006#ifdef ALI_DATE
1007 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1008 Bool_t result=ProcessEventDAQ(rawReader);
1009 delete rawReader;
1010 return result;
1011#else
1012 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1013 return 0;
1014#endif
1015
1016}
1017//____________Online trackling in AliTRDtrigger________________________________
1018Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1019{
1020 //
1021 // For the DAQ
1022 // Fill a simple average pulse height
1023 //
1024
1025 // Localisation of the Xbins involved
1026 //LocalisationDetectorXbins(det);
1027
1028 // Row and position in the pad groups
1029 //Int_t posr = 0;
1030 //if (fCalibraMode->GetNnZ(1) != 0) {
1031 // posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1032 //}
1033
1034 // Col of cluster and position in the pad groups
1035 //Int_t posc = 0;
1036 //if (fCalibraMode->GetNnRphi(1) != 0) {
1037 // posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1038 //}
1039
1040 //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1041
1042 ((TProfile2D *)GetPH2d(nbtimebins,fSf,kTRUE))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1043
1044 return kTRUE;
1045
1046}
55a288e5 1047//____________Functions for plotting the 2D____________________________________
1048
1049//_____________________________________________________________________________
1050void AliTRDCalibraFillHisto::Plot2d()
1051{
1052 //
1053 // Plot the 2D histos
1054 //
1055
1056 if (fPH2dOn) {
1057 TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
1058 cph2d->cd();
1059 fPH2d->Draw("LEGO");
1060 }
1061 if (fCH2dOn) {
1062 TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
1063 cch2d->cd();
1064 fCH2d->Draw("LEGO");
1065 }
1066 if (fPRF2dOn) {
1067 TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
1068 cPRF2d->cd();
1069 fPRF2d->Draw("LEGO");
1070 }
1071
1072}
170c35f1 1073//_____________________Reset____________________________________________________
1074//_______________________________________________________________________________
1075void AliTRDCalibraFillHisto::ResetLinearFitter()
55a288e5 1076{
170c35f1 1077 fLinearFitterArray.SetOwner();
1078 fLinearFitterArray.Clear();
1079 fLinearFitterHistoArray.SetOwner();
1080 fLinearFitterHistoArray.Clear();
1081}
1082//____________Write_____________________________________________________
1083//_____________________________________________________________________
1084void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1085{
1086 //
1087 // Write infos to file
1088 //
1089
1090 //For debugging
1091 if ( fDebugStreamer ) {
1092 delete fDebugStreamer;
1093 fDebugStreamer = 0x0;
55a288e5 1094 }
170c35f1 1095
55a288e5 1096 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
170c35f1 1097 ,fNumberTrack
1098 ,fNumberUsedCh[0]
1099 ,fNumberUsedCh[1]
1100 ,fNumberUsedPh[0]
1101 ,fNumberUsedPh[1]));
1102
1103 TDirectory *backup = gDirectory;
1104 TString option;
1105
1106 if ( append )
1107 option = "update";
1108 else
1109 option = "recreate";
1110
1111 TFile f(filename,option.Data());
55a288e5 1112
1113 TStopwatch stopwatch;
1114 stopwatch.Start();
55a288e5 1115
170c35f1 1116 if (fCH2dOn ) {
55a288e5 1117 if (fHisto2d) {
170c35f1 1118 f.WriteTObject(fCH2d);
55a288e5 1119 }
1120 if (fVector2d) {
1121 TString name("Nz");
1122 name += fCalibraMode->GetNz(0);
1123 name += "Nrphi";
1124 name += fCalibraMode->GetNrphi(0);
1125 TTree *treeCH2d = fCalibraVector->ConvertVectorCTTreeHisto(fCalibraVector->GetVectorCH(),fCalibraVector->GetPlaCH(),"treeCH2d",(const char *) name);
170c35f1 1126 f.WriteTObject(treeCH2d);
55a288e5 1127 }
1128 }
170c35f1 1129 if (fPH2dOn ) {
55a288e5 1130 if (fHisto2d) {
170c35f1 1131 f.WriteTObject(fPH2d);
55a288e5 1132 }
1133 if (fVector2d) {
1134 TString name("Nz");
1135 name += fCalibraMode->GetNz(1);
1136 name += "Nrphi";
1137 name += fCalibraMode->GetNrphi(1);
1138 TTree *treePH2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPH(),fCalibraVector->GetPlaPH(),"treePH2d",(const char *) name);
170c35f1 1139 f.WriteTObject(treePH2d);
55a288e5 1140 }
1141 }
170c35f1 1142 if (fPRF2dOn) {
55a288e5 1143 if (fHisto2d) {
170c35f1 1144 f.WriteTObject(fPRF2d);
55a288e5 1145 }
1146 if (fVector2d) {
1147 TString name("Nz");
1148 name += fCalibraMode->GetNz(2);
1149 name += "Nrphi";
1150 name += fCalibraMode->GetNrphi(2);
170c35f1 1151 name += "Ngp";
1152 name += fNgroupprf;
55a288e5 1153 TTree *treePRF2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPRF(),fCalibraVector->GetPlaPRF(),"treePRF2d",(const char *) name);
170c35f1 1154 f.WriteTObject(treePRF2d);
55a288e5 1155 }
1156 }
170c35f1 1157 if(fLinearFitterOn && fLinearFitterDebugOn){
1158 f.WriteTObject(&fLinearFitterHistoArray);
1159 }
55a288e5 1160
170c35f1 1161 f.Close();
55a288e5 1162
170c35f1 1163 if ( backup ) backup->cd();
55a288e5 1164
170c35f1 1165 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1166 ,stopwatch.RealTime(),stopwatch.CpuTime()));
55a288e5 1167}
170c35f1 1168//___________________________________________probe the histos__________________________________________________
55a288e5 1169Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1170{
1171 //
1172 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1173 // debug mode with 2 for TH2I and 3 for TProfile2D
1174 // It gives a pointer to a Double_t[7] with the info following...
1175 // [0] : number of calibration groups with entries
1176 // [1] : minimal number of entries found
1177 // [2] : calibration group number of the min
1178 // [3] : maximal number of entries found
1179 // [4] : calibration group number of the max
1180 // [5] : mean number of entries found
1181 // [6] : mean relativ error
1182 //
1183
1184 Double_t *info = new Double_t[7];
1185
1186 // Number of Xbins (detectors or groups of pads)
170c35f1 1187 Int_t nbins = h->GetNbinsY(); //number of calibration groups
1188 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
55a288e5 1189
1190 // Initialise
1191 Double_t nbwe = 0; //number of calibration groups with entries
1192 Double_t minentries = 0; //minimal number of entries found
1193 Double_t maxentries = 0; //maximal number of entries found
1194 Double_t placemin = 0; //calibration group number of the min
1195 Double_t placemax = -1; //calibration group number of the max
1196 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1197 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1198
1199 Double_t counter = 0;
1200
1201 //Debug
1202 TH1F *NbEntries = 0x0;//distribution of the number of entries
1203 TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1204 TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1205
1206 // Beginning of the loop over the calibration groups
1207 for (Int_t idect = 0; idect < nbins; idect++) {
1208
170c35f1 1209 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
55a288e5 1210 projch->SetDirectory(0);
1211
1212 // Number of entries for this calibration group
1213 Double_t nentries = 0.0;
1214 if((i%2) == 0){
170c35f1 1215 for (Int_t k = 0; k < nxbins; k++) {
1216 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
55a288e5 1217 }
1218 }
1219 else{
170c35f1 1220 for (Int_t k = 0; k < nxbins; k++) {
1221 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1222 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1223 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
55a288e5 1224 counter++;
1225 }
1226 }
1227 }
1228
1229 //Debug
1230 if(i > 1){
1231 if((!((Bool_t)NbEntries)) && (nentries > 0)){
1232 NbEntries = new TH1F("Number of entries","Number of entries"
1233 ,100,(Int_t)nentries/2,nentries*2);
1234 NbEntries->SetDirectory(0);
1235 NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1236 ,nbins,0,nbins);
1237 NbEntriesPerGroup->SetDirectory(0);
1238 NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1239 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1240 NbEntriesPerSp->SetDirectory(0);
1241 }
1242 if(NbEntries){
1243 if(nentries > 0) NbEntries->Fill(nentries);
1244 NbEntriesPerGroup->Fill(idect+0.5,nentries);
1245 NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1246 }
1247 }
1248
1249 //min amd max
1250 if(nentries > maxentries){
1251 maxentries = nentries;
1252 placemax = idect;
1253 }
1254 if(idect == 0) {
1255 minentries = nentries;
1256 }
1257 if(nentries < minentries){
1258 minentries = nentries;
1259 placemin = idect;
1260 }
1261 //nbwe
1262 if(nentries > 0) {
1263 nbwe++;
1264 meanstats += nentries;
1265 }
55a288e5 1266 }//calibration groups loop
1267
1268 if(nbwe > 0) meanstats /= nbwe;
1269 if(counter > 0) meanrelativerror /= counter;
1270
1271 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1272 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1273 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1274 AliInfo(Form("The mean number of entries is %f",meanstats));
1275 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1276
1277 info[0] = nbwe;
1278 info[1] = minentries;
1279 info[2] = placemin;
1280 info[3] = maxentries;
1281 info[4] = placemax;
1282 info[5] = meanstats;
1283 info[6] = meanrelativerror;
1284
1285 if(i > 1){
1286 gStyle->SetPalette(1);
1287 gStyle->SetOptStat(1111);
1288 gStyle->SetPadBorderMode(0);
1289 gStyle->SetCanvasColor(10);
1290 gStyle->SetPadLeftMargin(0.13);
1291 gStyle->SetPadRightMargin(0.01);
1292 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1293 stat->Divide(2,1);
1294 stat->cd(1);
1295 NbEntries->Draw("");
1296 stat->cd(2);
1297 NbEntriesPerSp->SetStats(0);
1298 NbEntriesPerSp->Draw("");
1299 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1300 stat1->cd();
1301 NbEntriesPerGroup->SetStats(0);
1302 NbEntriesPerGroup->Draw("");
1303 }
1304
1305 return info;
1306
1307}
170c35f1 1308//____________________________________________________________________________
1309Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1310{
1311 //
1312 // Return a Int_t[4] with:
1313 // 0 Mean number of entries
1314 // 1 median of number of entries
1315 // 2 rms of number of entries
1316 // 3 number of group with entries
1317 //
1318
1319 Double_t *stat = new Double_t[4];
1320 stat[3] = 0.0;
1321
1322 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
1323 Double_t *weight = new Double_t[nbofgroups];
1324 Int_t *nonul = new Int_t[nbofgroups];
1325
1326 for(Int_t k = 0; k < nbofgroups; k++){
1327 if(fEntriesCH[k] > 0) {
1328 weight[k] = 1.0;
1329 nonul[(Int_t)stat[3]] = fEntriesCH[k];
1330 stat[3]++;
1331 }
1332 else weight[k] = 0.0;
1333 }
1334 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
1335 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
1336 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1337
1338 return stat;
1339
1340}
1341//____________________________________________________________________________
1342Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1343{
1344 //
1345 // Return a Int_t[4] with:
1346 // 0 Mean number of entries
1347 // 1 median of number of entries
1348 // 2 rms of number of entries
1349 // 3 number of group with entries
1350 //
1351
1352 Double_t *stat = new Double_t[4];
1353 stat[3] = 0.0;
1354
1355 Int_t nbofgroups = 540;
1356 Double_t *weight = new Double_t[nbofgroups];
1357 Int_t *nonul = new Int_t[nbofgroups];
1358
1359 for(Int_t k = 0; k < nbofgroups; k++){
1360 if(fEntriesLinearFitter[k] > 0) {
1361 weight[k] = 1.0;
1362 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1363 stat[3]++;
1364 }
1365 else weight[k] = 0.0;
1366 }
1367 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1368 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1369 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1370
1371 return stat;
1372
1373}
1374//_____________________________________________________________________________
1375void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1376{
1377 //
1378 // Should be between 0 and 6
1379 //
1380
1381 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1382 AliInfo("The number of groups must be between 0 and 6!");
1383 }
1384 else {
1385 fNgroupprf = numberGroupsPRF;
1386 }
55a288e5 1387
170c35f1 1388}
55a288e5 1389//_____________________________________________________________________________
1390void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1391{
1392 //
1393 // Set the factor that will divide the deposited charge
1394 // to fit in the histo range [0,300]
1395 //
1396
1397 if (RelativeScale > 0.0) {
1398 fRelativeScale = RelativeScale;
1399 }
1400 else {
1401 AliInfo("RelativeScale must be strict positif!");
1402 }
1403
1404}
1405
1406//_____________________________________________________________________________
1407void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1408{
1409 //
1410 // Set the mode of calibration group in the z direction for the parameter i
1411 //
1412
1413 if ((Nz >= 0) &&
1414 (Nz < 5)) {
1415 fCalibraMode->SetNz(i, Nz);
1416 }
1417 else {
1418 AliInfo("You have to choose between 0 and 4");
1419 }
1420
1421}
1422
1423//_____________________________________________________________________________
1424void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1425{
1426 //
1427 // Set the mode of calibration group in the rphi direction for the parameter i
1428 //
1429
1430 if ((Nrphi >= 0) &&
1431 (Nrphi < 7)) {
1432 fCalibraMode->SetNrphi(i ,Nrphi);
1433 }
1434 else {
1435 AliInfo("You have to choose between 0 and 6");
1436 }
1437
1438}
55a288e5 1439//____________Protected Functions______________________________________________
1440//____________Create the 2D histo to be filled online__________________________
1441//
55a288e5 1442//_____________________________________________________________________________
1443void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1444{
1445 //
170c35f1 1446 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1447 // If fNgroupprf is zero then no binning in tnp
55a288e5 1448 //
1449
1450 TString name("Nz");
1451 name += fCalibraMode->GetNz(2);
1452 name += "Nrphi";
1453 name += fCalibraMode->GetNrphi(2);
170c35f1 1454 name += "Ngp";
1455 name += fNgroupprf;
55a288e5 1456
170c35f1 1457 if(fNgroupprf != 0){
1458
1459 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1460 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1461 fPRF2d->SetYTitle("Det/pad groups");
1462 fPRF2d->SetXTitle("Position x/W [pad width units]");
1463 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1464 fPRF2d->SetStats(0);
1465 }
1466 else{
1467 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1468 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1469 fPRF2d->SetYTitle("Det/pad groups");
1470 fPRF2d->SetXTitle("Position x/W [pad width units]");
1471 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1472 fPRF2d->SetStats(0);
1473 }
55a288e5 1474
1475}
1476
1477//_____________________________________________________________________________
1478void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1479{
1480 //
1481 // Create the 2D histos
1482 //
1483
1484 TString name("Nz");
1485 name += fCalibraMode->GetNz(1);
1486 name += "Nrphi";
1487 name += fCalibraMode->GetNrphi(1);
170c35f1 1488
55a288e5 1489 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
170c35f1 1490 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1491 ,nn,0,nn);
1492 fPH2d->SetYTitle("Det/pad groups");
1493 fPH2d->SetXTitle("time [#mus]");
55a288e5 1494 fPH2d->SetZTitle("<PH> [a.u.]");
1495 fPH2d->SetStats(0);
1496
1497}
1498
1499//_____________________________________________________________________________
1500void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1501{
1502 //
1503 // Create the 2D histos
1504 //
1505
1506 TString name("Nz");
1507 name += fCalibraMode->GetNz(0);
1508 name += "Nrphi";
1509 name += fCalibraMode->GetNrphi(0);
1510
1511 fCH2d = new TH2I("CH2d",(const Char_t *) name
170c35f1 1512 ,fNumberBinCharge,0,300,nn,0,nn);
1513 fCH2d->SetYTitle("Det/pad groups");
1514 fCH2d->SetXTitle("charge deposit [a.u]");
55a288e5 1515 fCH2d->SetZTitle("counts");
1516 fCH2d->SetStats(0);
1517 fCH2d->Sumw2();
1518
1519}
55a288e5 1520//____________Offine tracking in the AliTRDtracker_____________________________
1521void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1522{
1523 //
1524 // For the offline tracking or mcm tracklets
1525 // This function will be called in the functions UpdateHistogram...
1526 // to fill the info of a track for the relativ gain calibration
1527 //
1528
170c35f1 1529 Int_t nb = 0; // Nombre de zones traversees
1530 Int_t fd = -1; // Premiere zone non nulle
1531 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
55a288e5 1532
1533
1534 // See if the track goes through different zones
1535 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1536 if (fAmpTotal[k] > 0.0) {
170c35f1 1537 totalcharge += fAmpTotal[k];
55a288e5 1538 nb++;
1539 if (nb == 1) {
1540 fd = k;
1541 }
1542 }
1543 }
55a288e5 1544
170c35f1 1545
1546 switch (nb)
1547 {
1548 case 1:
55a288e5 1549 fNumberUsedCh[0]++;
170c35f1 1550 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 1551 if (fHisto2d) {
170c35f1 1552 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1553 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 1554 }
1555 if (fVector2d) {
170c35f1 1556 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 1557 }
170c35f1 1558 break;
1559 case 2:
55a288e5 1560 if ((fAmpTotal[fd] > 0.0) &&
170c35f1 1561 (fAmpTotal[fd+1] > 0.0)) {
55a288e5 1562 // One of the two very big
1563 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1564 if (fHisto2d) {
170c35f1 1565 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1566 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 1567 }
1568 if (fVector2d) {
170c35f1 1569 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 1570 }
1571 fNumberUsedCh[1]++;
170c35f1 1572 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 1573 }
1574 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1575 if (fHisto2d) {
170c35f1 1576 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1577 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
55a288e5 1578 }
1579 if (fVector2d) {
170c35f1 1580 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
55a288e5 1581 }
1582 fNumberUsedCh[1]++;
170c35f1 1583 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
55a288e5 1584 }
1585 }
55a288e5 1586 if (fCalibraMode->GetNfragZ(0) > 1) {
1587 if (fAmpTotal[fd] > 0.0) {
1588 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1589 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1590 // One of the two very big
1591 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1592 if (fHisto2d) {
170c35f1 1593 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1594 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 1595 }
1596 if (fVector2d) {
170c35f1 1597 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 1598 }
1599 fNumberUsedCh[1]++;
170c35f1 1600 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 1601 }
1602 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1603 if (fHisto2d) {
170c35f1 1604 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1605 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
55a288e5 1606 }
1607 fNumberUsedCh[1]++;
170c35f1 1608 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
55a288e5 1609 if (fVector2d) {
170c35f1 1610 fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
55a288e5 1611 }
1612 }
1613 }
1614 }
1615 }
1616 }
170c35f1 1617 break;
1618 default: break;
55a288e5 1619 }
55a288e5 1620}
55a288e5 1621//____________Offine tracking in the AliTRDtracker_____________________________
1622void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1623{
1624 //
1625 // For the offline tracking or mcm tracklets
1626 // This function will be called in the functions UpdateHistogram...
1627 // to fill the info of a track for the drift velocity calibration
1628 //
1629
1630 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
1631 Int_t fd1 = -1; // Premiere zone non nulle
1632 Int_t fd2 = -1; // Deuxieme zone non nulle
1633 Int_t k1 = -1; // Debut de la premiere zone
1634 Int_t k2 = -1; // Debut de la seconde zone
1635
1636 // See if the track goes through different zones
1637 for (Int_t k = 0; k < fTimeMax; k++) {
1638 if (fPHValue[k] > 0.0) {
1639 if (fd1 == -1) {
1640 fd1 = fPHPlace[k];
1641 k1 = k;
1642 }
1643 if (fPHPlace[k] != fd1) {
1644 if (fd2 == -1) {
1645 k2 = k;
1646 fd2 = fPHPlace[k];
1647 nb = 2;
1648 }
1649 if (fPHPlace[k] != fd2) {
1650 nb = 3;
1651 }
1652 }
1653 }
1654 }
170c35f1 1655
1656 switch(nb)
1657 {
1658 case 1:
1659 fNumberUsedPh[0]++;
1660 for (Int_t i = 0; i < fTimeMax; i++) {
1661 if (fHisto2d) {
1662 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1663 }
1664 if (fVector2d) {
1665 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
55a288e5 1666 }
1667 }
170c35f1 1668 break;
1669 case 2:
1670 if ((fd1 == fd2+1) ||
1671 (fd2 == fd1+1)) {
1672 // One of the two fast all the think
1673 if (k2 > (k1+fDifference)) {
1674 //we choose to fill the fd1 with all the values
1675 fNumberUsedPh[1]++;
1676 for (Int_t i = 0; i < fTimeMax; i++) {
1677 if (fHisto2d) {
1678 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1679 }
1680 if (fVector2d) {
1681 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1682 }
55a288e5 1683 }
170c35f1 1684 }
1685 if ((k2+fDifference) < fTimeMax) {
1686 //we choose to fill the fd2 with all the values
1687 fNumberUsedPh[1]++;
1688 for (Int_t i = 0; i < fTimeMax; i++) {
1689 if (fHisto2d) {
1690 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1691 }
55a288e5 1692 if (fVector2d) {
1693 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1694 }
170c35f1 1695 }
55a288e5 1696 }
1697 }
170c35f1 1698 // Two zones voisines sinon rien!
1699 if (fCalibraMode->GetNfragZ(1) > 1) {
1700 // Case 2
1701 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1702 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1703 // One of the two fast all the think
1704 if (k2 > (k1+fDifference)) {
1705 //we choose to fill the fd1 with all the values
1706 fNumberUsedPh[1]++;
1707 for (Int_t i = 0; i < fTimeMax; i++) {
1708 if (fHisto2d) {
1709 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1710 }
1711 if (fVector2d) {
1712 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1713 }
55a288e5 1714 }
1715 }
170c35f1 1716 if ((k2+fDifference) < fTimeMax) {
1717 //we choose to fill the fd2 with all the values
1718 fNumberUsedPh[1]++;
1719 for (Int_t i = 0; i < fTimeMax; i++) {
1720 if (fHisto2d) {
1721 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1722 }
1723 if (fVector2d) {
1724 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1725 }
55a288e5 1726 }
1727 }
1728 }
1729 }
170c35f1 1730 // Two zones voisines sinon rien!
1731 // Case 3
1732 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1733 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1734 // One of the two fast all the think
1735 if (k2 > (k1 + fDifference)) {
1736 //we choose to fill the fd1 with all the values
1737 fNumberUsedPh[1]++;
1738 for (Int_t i = 0; i < fTimeMax; i++) {
1739 if (fHisto2d) {
1740 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1741 }
1742 if (fVector2d) {
1743 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1744 }
55a288e5 1745 }
1746 }
170c35f1 1747 if ((k2+fDifference) < fTimeMax) {
1748 //we choose to fill the fd2 with all the values
1749 fNumberUsedPh[1]++;
1750 for (Int_t i = 0; i < fTimeMax; i++) {
1751 if (fHisto2d) {
1752 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1753 }
1754 if (fVector2d) {
1755 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1756 }
55a288e5 1757 }
1758 }
1759 }
1760 }
1761 }
170c35f1 1762 break;
1763 default: break;
1764 }
55a288e5 1765}
170c35f1 1766//____________Offine tracking in the AliTRDtracker_____________________________
1767Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
55a288e5 1768{
1769 //
170c35f1 1770 // For the offline tracking
1771 // This function will be called in the functions UpdateHistogram...
1772 // to fill the find the parameter P1 of a track for the drift velocity calibration
55a288e5 1773 //
170c35f1 1774
1775 // Get the parameter object
1776 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1777 if (!parCom) {
1778 AliInfo("Could not get CommonParam");
1779 return kFALSE;
55a288e5 1780 }
170c35f1 1781
1782 //Number of points: if less than 3 return kFALSE
1783 Int_t Npoints = fListClusters->GetEntriesFast();
1784 if(Npoints <= 2) return kFALSE;
1785
1786 //Variables
1787 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
1788 Double_t snp = 0.0; // sin angle in the plan yx track
1789 Double_t y = 0.0; // y clusters in the middle of the chamber
1790 Double_t z = 0.0; // z cluster in the middle of the chamber
1791 Double_t dydt = 0.0; // dydt tracklet after straight line fit
1792 Double_t tnp = 0.0; // tan angle in the plan xy track
1793 Double_t tgl = 0.0; // dz/dl and not dz/dx!
1794 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
1795 Double_t pointError = 0.0; // error after straight line fit
1796 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
1797 Int_t snpright = 1; // if we took in the middle snp
1798 Int_t crossrow = 0; // if it crosses a pad row
1799 Double_t tiltingangle = 0; // tiltingangle of the pad
1800 Float_t dzdx = 0; // dz/dx now from dz/dl
1801 Int_t nbli = 0; // number linear fitter points
1802 AliTRDpadPlane *padplane = parCom->GetPadPlane(GetPlane(detector),GetChamber(detector));
1803
1804 linearFitterTracklet.StoreData(kFALSE);
1805 linearFitterTracklet.ClearPoints();
1806
1807 //if more than one row
1808 Int_t rowp = -1; // if it crosses a pad row
1809
1810 //tiltingangle
1811 tiltingangle = padplane->GetTiltingAngle();
1812 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
1813
1814 //Fill with points
1815 for(Int_t k = 0; k < Npoints; k++){
1816
1817 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1818 Double_t ycluster = cl->GetY();
1819 Int_t time = cl->GetLocalTimeBin();
1820 Double_t timeis = time/fSf;
1821 //See if cross two pad rows
1822 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
1823 if(k==0) rowp = row;
1824 if(row != rowp) crossrow = 1;
1825 //Take in the middle of the chamber
1826 if(time > (Int_t) 10) {
1827 z = cl->GetZ();
1828 y = cl->GetY();
1829 snp = fPar2[time];
1830 tgl = fPar3[time];
1831 }
1832 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1833 nbli++;
55a288e5 1834 }
170c35f1 1835 if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() <= 10) snpright = 0;
1836 if(nbli <= 2) return kFALSE;
1837
1838 // Do the straight line fit now
1839 TVectorD pars;
1840 linearFitterTracklet.Eval();
1841 linearFitterTracklet.GetParameters(pars);
1842 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
1843 errorpar = linearFitterTracklet.GetParError(1)*pointError;
1844 dydt = pars[1];
1845
1846 if( TMath::Abs(snp) < 1.){
1847 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1848 }
1849 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1850
1851 if(fDebugLevel > 0){
1852 if ( !fDebugStreamer ) {
1853 //debug stream
1854 TDirectory *backup = gDirectory;
1855 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
1856 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1857 }
1858
1859 (* fDebugStreamer) << "VDRIFT0"<<
1860 "Npoints="<<Npoints<<
1861 "\n";
1862
1863
1864 (* fDebugStreamer) << "VDRIFT"<<
1865 "snpright="<<snpright<<
1866 "Npoints="<<Npoints<<
1867 "nbli="<<nbli<<
1868 "detector="<<detector<<
1869 "snp="<<snp<<
1870 "tnp="<<tnp<<
1871 "tgl="<<tgl<<
1872 "tnt="<<tnt<<
1873 "y="<<y<<
1874 "z="<<z<<
1875 "dydt="<<dydt<<
1876 "dzdx="<<dzdx<<
1877 "crossrow="<<crossrow<<
1878 "errorpar="<<errorpar<<
1879 "pointError="<<pointError<<
1880 "\n";
55a288e5 1881
55a288e5 1882 }
1883
170c35f1 1884 if(Npoints < fNumberClusters) return kFALSE;
1885 if(snpright == 0) return kFALSE;
1886 if(pointError >= 0.1) return kFALSE;
1887 if(crossrow == 1) return kFALSE;
1888
1889 if(fLinearFitterOn){
1890 //Add to the linear fitter of the detector
1891 if( TMath::Abs(snp) < 1.){
1892 Double_t x = tnp-dzdx*tnt;
1893 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
1894 if(fLinearFitterDebugOn) {
1895 ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars[1]);
1896 }
1897 fEntriesLinearFitter[detector]++;
1898 }
1899 }
1900 //AliInfo("End of FindP1TrackPH with success!")
55a288e5 1901 return kTRUE;
1902
1903}
170c35f1 1904//____________Offine tracking in the AliTRDtracker_____________________________
1905Bool_t AliTRDCalibraFillHisto::HandlePRF()
1906{
1907 //
1908 // For the offline tracking
1909 // Fit the tracklet with a line and take the position as reference for the PRF
1910 //
1911
1912 // Get the parameter object
1913 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
1914 if (!parCom) {
1915 AliInfo("Could not get CommonParam");
1916 return kFALSE;
1917 }
1918
1919 //Number of points
1920 Int_t Npoints = fListClusters->GetEntriesFast(); // number of total points
1921 Int_t Nb3pc = 0; // number of three pads clusters used for fit
1922 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
1923
1924
1925 // To see the difference due to the fit
1926 Double_t *padPositions;
1927 padPositions = new Double_t[Npoints];
1928 for(Int_t k = 0; k < Npoints; k++){
1929 padPositions[k] = 0.0;
1930 }
1931
1932
1933 //Find the position by a fit
1934 TLinearFitter fitter(2,"pol1");
1935 fitter.StoreData(kFALSE);
1936 fitter.ClearPoints();
1937 for(Int_t k = 0; k < Npoints; k++){
1938 //Take the cluster
1939 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1940 Short_t *signals = cl->GetSignals();
1941 Double_t time = cl->GetLocalTimeBin();
1942 //Calculate x if possible
1943 Float_t xcenter = 0.0;
1944 Bool_t echec = kTRUE;
1945 if((time<=7) || (time>=21)) continue;
1946 // Center 3 balanced: position with the center of the pad
1947 if ((((Float_t) signals[3]) > 0.0) &&
1948 (((Float_t) signals[2]) > 0.0) &&
1949 (((Float_t) signals[4]) > 0.0)) {
1950 echec = kFALSE;
1951 // Security if the denomiateur is 0
1952 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1953 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1954 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1955 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1956 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1957 }
1958 else {
1959 echec = kTRUE;
1960 }
1961 }
1962 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1963 if(echec) continue;
1964 //if no echec: calculate with the position of the pad
1965 // Position of the cluster
1966 Double_t padPosition = xcenter + cl->GetPad();
1967 padPositions[k] = padPosition;
1968 Nb3pc++;
1969 fitter.AddPoint(&time, padPosition,1);
1970 }//clusters loop
1971
1972 //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
1973 if(Nb3pc < 3) return kFALSE;
1974 fitter.Eval();
1975 TVectorD line(2);
1976 fitter.GetParameters(line);
1977 Float_t pointError = -1.0;
1978 pointError = TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
1979
1980
1981 // Now fill the PRF
1982 for(Int_t k = 0; k < Npoints; k++){
1983 //Take the cluster
1984 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
1985 Short_t *signals = cl->GetSignals(); // signal
1986 Double_t time = cl->GetLocalTimeBin(); // time bin
1987 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
1988 Float_t padPos = cl->GetPad(); // middle pad
1989 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
1990 Float_t ycenter = 0.0; // relative center charge
1991 Float_t ymin = 0.0; // relative left charge
1992 Float_t ymax = 0.0; // relative right charge
1993 Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx
1994 Double_t pt = fPar4[(Int_t)time]; // pt
1995 Float_t dzdx = 0.0; // dzdx
1996
1997
1998 //Requiere simply two pads clusters at least
1999 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2000 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2001 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2002 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2003 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2004 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2005 }
2006
2007 //calibration group
2008 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2009 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2010 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponidng group
2011 Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan
2012 Float_t xcl = cl->GetY(); // y cluster
2013 Float_t qcl = cl->GetQ(); // charge cluster
2014 Int_t plane = GetPlane(detector); // plane
2015 Int_t chamber = GetChamber(detector); // chamber
2016 Double_t xdiff = dpad; // reconstructed position constant
2017 Double_t x = dpad; // reconstructed position moved
2018 Float_t Ep = pointError; // error of fit
2019 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2020 Float_t signal3 = (Float_t)signals[3]; // signal
2021 Float_t signal2 = (Float_t)signals[2]; // signal
2022 Float_t signal4 = (Float_t)signals[4]; // signal
2023 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2024 Float_t tnp = 0.0;
2025 if(TMath::Abs(snp) < 1.0){
2026 tnp = snp / (TMath::Sqrt(1-snp*snp));
2027 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2028 }
2029
55a288e5 2030
170c35f1 2031 if(fDebugLevel > 0){
2032 if ( !fDebugStreamer ) {
2033 //debug stream
2034 TDirectory *backup = gDirectory;
2035 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
2036 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2037 }
2038
2039 (* fDebugStreamer) << "PRF0"<<
2040 "caligroup="<<caligroup<<
2041 "detector="<<detector<<
2042 "plane="<<plane<<
2043 "chamber="<<chamber<<
2044 "Npoints="<<Npoints<<
2045 "Np="<<Nb3pc<<
2046 "Ep="<<Ep<<
2047 "snp="<<snp<<
2048 "tnp="<<tnp<<
2049 "tgl="<<tgl<<
2050 "dzdx="<<dzdx<<
2051 "pt="<<pt<<
2052 "padPos="<<padPos<<
2053 "padPosition="<<padPositions[k]<<
2054 "padPosTracklet="<<padPosTracklet<<
2055 "x="<<x<<
2056 "ycenter="<<ycenter<<
2057 "ymin="<<ymin<<
2058 "ymax="<<ymax<<
2059 "xcl="<<xcl<<
2060 "qcl="<<qcl<<
2061 "signal1="<<signal1<<
2062 "signal2="<<signal2<<
2063 "signal3="<<signal3<<
2064 "signal4="<<signal4<<
2065 "signal5="<<signal5<<
2066 "time="<<time<<
2067 "\n";
2068 x = xdiff;
2069 Int_t type=0;
2070 Float_t y = ycenter;
2071 (* fDebugStreamer) << "PRFALL"<<
2072 "caligroup="<<caligroup<<
2073 "detector="<<detector<<
2074 "plane="<<plane<<
2075 "chamber="<<chamber<<
2076 "Npoints="<<Npoints<<
2077 "Np="<<Nb3pc<<
2078 "Ep="<<Ep<<
2079 "type="<<type<<
2080 "snp="<<snp<<
2081 "tnp="<<tnp<<
2082 "tgl="<<tgl<<
2083 "dzdx="<<dzdx<<
2084 "pt="<<pt<<
2085 "padPos="<<padPos<<
2086 "padPosition="<<padPositions[k]<<
2087 "padPosTracklet="<<padPosTracklet<<
2088 "x="<<x<<
2089 "y="<<y<<
2090 "xcl="<<xcl<<
2091 "qcl="<<qcl<<
2092 "signal1="<<signal1<<
2093 "signal2="<<signal2<<
2094 "signal3="<<signal3<<
2095 "signal4="<<signal4<<
2096 "signal5="<<signal5<<
2097 "time="<<time<<
2098 "\n";
2099 x=-(xdiff+1);
2100 y = ymin;
2101 type=-1;
2102 (* fDebugStreamer) << "PRFALL"<<
2103 "caligroup="<<caligroup<<
2104 "detector="<<detector<<
2105 "plane="<<plane<<
2106 "chamber="<<chamber<<
2107 "Npoints="<<Npoints<<
2108 "Np="<<Nb3pc<<
2109 "Ep="<<Ep<<
2110 "type="<<type<<
2111 "snp="<<snp<<
2112 "tnp="<<tnp<<
2113 "tgl="<<tgl<<
2114 "dzdx="<<dzdx<<
2115 "pt="<<pt<<
2116 "padPos="<<padPos<<
2117 "padPosition="<<padPositions[k]<<
2118 "padPosTracklet="<<padPosTracklet<<
2119 "x="<<x<<
2120 "y="<<y<<
2121 "xcl="<<xcl<<
2122 "qcl="<<qcl<<
2123 "signal1="<<signal1<<
2124 "signal2="<<signal2<<
2125 "signal3="<<signal3<<
2126 "signal4="<<signal4<<
2127 "signal5="<<signal5<<
2128 "time="<<time<<
2129 "\n";
2130 x=1-xdiff;
2131 y = ymax;
2132 type=1;
2133
2134 (* fDebugStreamer) << "PRFALL"<<
2135 "caligroup="<<caligroup<<
2136 "detector="<<detector<<
2137 "plane="<<plane<<
2138 "chamber="<<chamber<<
2139 "Npoints="<<Npoints<<
2140 "Np="<<Nb3pc<<
2141 "Ep="<<Ep<<
2142 "type="<<type<<
2143 "snp="<<snp<<
2144 "tnp="<<tnp<<
2145 "tgl="<<tgl<<
2146 "dzdx="<<dzdx<<
2147 "pt="<<pt<<
2148 "padPos="<<padPos<<
2149 "padPosition="<<padPositions[k]<<
2150 "padPosTracklet="<<padPosTracklet<<
2151 "x="<<x<<
2152 "y="<<y<<
2153 "xcl="<<xcl<<
2154 "qcl="<<qcl<<
2155 "signal1="<<signal1<<
2156 "signal2="<<signal2<<
2157 "signal3="<<signal3<<
2158 "signal4="<<signal4<<
2159 "signal5="<<signal5<<
2160 "time="<<time<<
2161 "\n";
2162
2163 }
2164
2165 // some cuts
2166 if(Npoints < fNumberClusters) continue;
2167 if(Nb3pc <= 5) continue;
2168 if((time >= 21) || (time < 7)) continue;
2169 if(TMath::Abs(snp) >= 1.0) continue;
2170 if(qcl < 80) continue;
2171
2172 Bool_t echec = kFALSE;
2173 Double_t shift = 0.0;
2174 //Calculate the shift in x coresponding to this tnp
2175 if(fNgroupprf != 0.0){
2176 shift = -3.0*(fNgroupprf-1)-1.5;
2177 Double_t limithigh = -0.2*(fNgroupprf-1);
2178 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2179 else{
2180 while(tnp > limithigh){
2181 limithigh += 0.2;
2182 shift += 3.0;
2183 }
2184 }
2185 }
2186 if (fHisto2d && !echec) {
2187 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2188 if(ymin > 0.0) {
2189 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2190 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2191 }
2192 if(ymax > 0.0) {
2193 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2194 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2195 }
2196 }
2197 //Not equivalent anymore here!
2198 if (fVector2d && !echec) {
2199 fCalibraVector->UpdateVectorPRF(caligroup,shift+dpad,ycenter);
2200 if(ymin > 0.0) {
2201 fCalibraVector->UpdateVectorPRF(caligroup,shift-(dpad+1.0),ymin);
2202 fCalibraVector->UpdateVectorPRF(caligroup,shift+(dpad+1.0),ymin);
2203 }
2204 if(ymax > 0.0) {
2205 fCalibraVector->UpdateVectorPRF(caligroup,shift+1.0-dpad,ymax);
2206 fCalibraVector->UpdateVectorPRF(caligroup,shift-1.0+dpad,ymax);
2207 }
2208 }
2209 }
2210 return kTRUE;
2211
2212}
55a288e5 2213//
2214//____________Some basic geometry function_____________________________________
2215//
55a288e5 2216//_____________________________________________________________________________
2217Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2218{
2219 //
2220 // Reconstruct the plane number from the detector number
2221 //
2222
2223 return ((Int_t) (d % 6));
2224
2225}
2226
2227//_____________________________________________________________________________
2228Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2229{
2230 //
2231 // Reconstruct the chamber number from the detector number
2232 //
2233 Int_t fgkNplan = 6;
2234
2235 return ((Int_t) (d % 30) / fgkNplan);
2236
2237}
2238
2239//_____________________________________________________________________________
2240Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2241{
2242 //
2243 // Reconstruct the sector number from the detector number
2244 //
2245 Int_t fg = 30;
2246
2247 return ((Int_t) (d / fg));
2248
2249}
170c35f1 2250//_____________________________________________________________________
2251TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency, Bool_t force)
2252{
2253 //
2254 // return pointer to fPH2d TProfile2D
2255 // if force is true create a new TProfile2D if it doesn't exist allready
2256 //
2257 if ( !force || fPH2d )
2258 return fPH2d;
2259
2260 // Some parameters
2261 fTimeMax = nbtimebin;
2262 fSf = samplefrequency;
2263
2264 /*
2265 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
2266 ,fCalibraMode->GetNz(1)
2267 ,fCalibraMode->GetNrphi(1)));
2268
2269 // Calcul the number of Xbins
2270 Int_t Ntotal1 = 0;
2271 fCalibraMode->ModePadCalibration(2,1);
2272 fCalibraMode->ModePadFragmentation(0,2,0,1);
2273 fCalibraMode->SetDetChamb2(1);
2274 Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
2275 fCalibraMode->ModePadCalibration(0,1);
2276 fCalibraMode->ModePadFragmentation(0,0,0,1);
2277 fCalibraMode->SetDetChamb0(1);
2278 Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
2279 AliInfo(Form("Total number of Xbins: %d",Ntotal1));
2280
2281 CreatePH2d(Ntotal1);
2282 */
2283
2284 CreatePH2d(540);
2285
2286 return fPH2d;
2287}
2288//_____________________________________________________________________
2289TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2290{
2291 //
2292 // return pointer to TLinearFitter Calibration
2293 // if force is true create a new TLinearFitter if it doesn't exist allready
2294 //
2295 if ( !force || fLinearFitterArray.UncheckedAt(detector) )
2296 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2297
2298 // if we are forced and TLinearFitter doesn't yet exist create it
2299
2300 // new TLinearFitter
2301 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2302 fLinearFitterArray.AddAt(linearfitter,detector);
2303 return linearfitter;
2304}
2305//_____________________________________________________________________
2306TH2F* AliTRDCalibraFillHisto::GetLinearFitterHisto(Int_t detector, Bool_t force)
2307{
2308 //
2309 // return pointer to TH2F histo
2310 // if force is true create a new histo if it doesn't exist allready
2311 //
2312 if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
2313 return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
2314
2315 // if we are forced and TLinearFitter doesn't yes exist create it
2316
2317 // new TH2F
2318 TString name("LFDV");
2319 name += detector;
2320
2321 TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
2322 ,100,-1.0,1.0,100
2323 ,-2.0,2.0);
2324 lfdv->SetXTitle("tan(phi_{track})");
2325 lfdv->SetYTitle("dy/dt");
2326 lfdv->SetZTitle("Number of clusters");
2327 lfdv->SetStats(0);
2328 lfdv->SetDirectory(0);
2329
2330 fLinearFitterHistoArray.AddAt(lfdv,detector);
2331 return lfdv;
2332}
2333//_____________________________________________________________________
2334void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2335{
2336 //
2337 // FillCH2d: Marian style
2338 //
2339
2340 //skip simply the value out of range
2341 if((y>=300.0) || (y<0.0)) return;
2342
2343 //Calcul the y place
2344 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2345 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2346
2347 //Fill
2348 fCH2d->GetArray()[place]++;
2349
2350}