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