]> 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
reading RAW without data
[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
55a288e5 36#include <TProfile2D.h>
37#include <TProfile.h>
38#include <TFile.h>
55a288e5 39#include <TStyle.h>
40#include <TCanvas.h>
55a288e5 41#include <TObjArray.h>
d0569428 42#include <TObject.h>
55a288e5 43#include <TH1F.h>
44#include <TH2I.h>
45#include <TH2.h>
46#include <TStopwatch.h>
47#include <TMath.h>
48#include <TDirectory.h>
170c35f1 49#include <TTreeStream.h>
50#include <TVectorD.h>
55a288e5 51
52#include "AliLog.h"
55a288e5 53
54#include "AliTRDCalibraFillHisto.h"
55#include "AliTRDCalibraMode.h"
56#include "AliTRDCalibraVector.h"
3a0f6479 57#include "AliTRDCalibraVdriftLinearFit.h"
55a288e5 58#include "AliTRDcalibDB.h"
59#include "AliTRDCommonParam.h"
60#include "AliTRDmcmTracklet.h"
61#include "AliTRDpadPlane.h"
62#include "AliTRDcluster.h"
63#include "AliTRDtrack.h"
d4232bb6 64#include "AliTRDrawStreamTB.h"
3a0f6479 65#include "AliRawReader.h"
66#include "AliRawReaderDate.h"
f162af62 67#include "AliTRDgeometry.h"
d0569428 68#include "./Cal/AliTRDCalROC.h"
69#include "./Cal/AliTRDCalDet.h"
170c35f1 70
71#ifdef ALI_DATE
72#include "event.h"
73#endif
55a288e5 74
3a0f6479 75
55a288e5 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}
3a0f6479 116
55a288e5 117//______________________________________________________________________________________
118AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
119 :TObject()
f162af62 120 ,fGeo(0)
55a288e5 121 ,fMITracking(kFALSE)
122 ,fMcmTracking(kFALSE)
123 ,fMcmCorrectAngle(kFALSE)
124 ,fCH2dOn(kFALSE)
125 ,fPH2dOn(kFALSE)
126 ,fPRF2dOn(kFALSE)
127 ,fHisto2d(kFALSE)
128 ,fVector2d(kFALSE)
170c35f1 129 ,fLinearFitterOn(kFALSE)
130 ,fLinearFitterDebugOn(kFALSE)
55a288e5 131 ,fRelativeScale(0)
170c35f1 132 ,fThresholdClusterPRF2(15.0)
133 ,fCalibraMode(new AliTRDCalibraMode())
134 ,fDebugStreamer(0)
135 ,fDebugLevel(0)
55a288e5 136 ,fDetectorAliTRDtrack(kFALSE)
55a288e5 137 ,fDetectorPreviousTrack(-1)
e4db522f 138 ,fMCMPrevious(-1)
139 ,fROBPrevious(-1)
170c35f1 140 ,fNumberClusters(18)
141 ,fProcent(6.0)
142 ,fDifference(17)
143 ,fNumberTrack(0)
144 ,fTimeMax(0)
145 ,fSf(10.0)
146 ,fNumberBinCharge(100)
147 ,fNumberBinPRF(40)
148 ,fNgroupprf(0)
149 ,fListClusters(new TObjArray())
150 ,fPar0(0x0)
151 ,fPar1(0x0)
152 ,fPar2(0x0)
153 ,fPar3(0x0)
154 ,fPar4(0x0)
55a288e5 155 ,fAmpTotal(0x0)
156 ,fPHPlace(0x0)
157 ,fPHValue(0x0)
170c35f1 158 ,fGoodTracklet(kTRUE)
159 ,fGoodTrack(kTRUE)
160 ,fEntriesCH(0x0)
161 ,fEntriesLinearFitter(0x0)
162 ,fCalibraVector(0x0)
55a288e5 163 ,fPH2d(0x0)
164 ,fPRF2d(0x0)
165 ,fCH2d(0x0)
d0569428 166 ,fLinearFitterArray(540)
3a0f6479 167 ,fLinearVdriftFit(0x0)
d0569428 168 ,fCalDetGain(0x0)
169 ,fCalROCGain(0x0)
170 ,fCalDetT0(0x0)
171 ,fCalROCT0(0x0)
55a288e5 172{
173 //
174 // Default constructor
175 //
176
170c35f1 177 //
178 // Init some default values
179 //
55a288e5 180
170c35f1 181 fNumberUsedCh[0] = 0;
182 fNumberUsedCh[1] = 0;
183 fNumberUsedPh[0] = 0;
184 fNumberUsedPh[1] = 0;
3a0f6479 185
f162af62 186 fGeo = new AliTRDgeometry();
3a0f6479 187
55a288e5 188}
f162af62 189
55a288e5 190//______________________________________________________________________________________
191AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
f162af62 192 :TObject(c)
193 ,fGeo(0)
194 ,fMITracking(c.fMITracking)
195 ,fMcmTracking(c.fMcmTracking)
196 ,fMcmCorrectAngle(c.fMcmCorrectAngle)
197 ,fCH2dOn(c.fCH2dOn)
198 ,fPH2dOn(c.fPH2dOn)
199 ,fPRF2dOn(c.fPRF2dOn)
200 ,fHisto2d(c.fHisto2d)
201 ,fVector2d(c.fVector2d)
202 ,fLinearFitterOn(c.fLinearFitterOn)
203 ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
204 ,fRelativeScale(c.fRelativeScale)
205 ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
206 ,fCalibraMode(0x0)
207 ,fDebugStreamer(0)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
210 ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
e4db522f 211 ,fMCMPrevious(c.fMCMPrevious)
212 ,fROBPrevious(c.fROBPrevious)
f162af62 213 ,fNumberClusters(c.fNumberClusters)
214 ,fProcent(c.fProcent)
215 ,fDifference(c.fDifference)
216 ,fNumberTrack(c.fNumberTrack)
217 ,fTimeMax(c.fTimeMax)
218 ,fSf(c.fSf)
219 ,fNumberBinCharge(c.fNumberBinCharge)
220 ,fNumberBinPRF(c.fNumberBinPRF)
221 ,fNgroupprf(c.fNgroupprf)
222 ,fListClusters(new TObjArray())
3a0f6479 223 ,fPar0(0x0)
224 ,fPar1(0x0)
225 ,fPar2(0x0)
226 ,fPar3(0x0)
227 ,fPar4(0x0)
228 ,fAmpTotal(0x0)
229 ,fPHPlace(0x0)
230 ,fPHValue(0x0)
f162af62 231 ,fGoodTracklet(c.fGoodTracklet)
232 ,fGoodTrack(c.fGoodTrack)
3a0f6479 233 ,fEntriesCH(0x0)
234 ,fEntriesLinearFitter(0x0)
f162af62 235 ,fCalibraVector(0x0)
236 ,fPH2d(0x0)
237 ,fPRF2d(0x0)
238 ,fCH2d(0x0)
d0569428 239 ,fLinearFitterArray(540)
3a0f6479 240 ,fLinearVdriftFit(0x0)
d0569428 241 ,fCalDetGain(0x0)
242 ,fCalROCGain(0x0)
243 ,fCalDetT0(0x0)
244 ,fCalROCT0(0x0)
55a288e5 245{
246 //
247 // Copy constructor
248 //
170c35f1 249 if(c.fCalibraMode) fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
250 if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
251 if(c.fPH2d) {
252 fPH2d = new TProfile2D(*c.fPH2d);
253 fPH2d->SetDirectory(0);
254 }
255 if(c.fPRF2d) {
256 fPRF2d = new TProfile2D(*c.fPRF2d);
257 fPRF2d->SetDirectory(0);
258 }
259 if(c.fCH2d) {
260 fCH2d = new TH2I(*c.fCH2d);
261 fCH2d->SetDirectory(0);
262 }
3a0f6479 263 if(c.fLinearVdriftFit){
264 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
170c35f1 265 }
3a0f6479 266
d0569428 267 if(c.fCalDetGain) fCalDetGain = new AliTRDCalDet(*c.fCalDetGain);
268 if(c.fCalDetT0) fCalDetT0 = new AliTRDCalDet(*c.fCalDetT0);
269 if(c.fCalROCGain) fCalROCGain = new AliTRDCalROC(*c.fCalROCGain);
270 if(c.fCalROCT0) fCalROCT0 = new AliTRDCalROC(*c.fCalROCT0);
271
f162af62 272 if (fGeo) {
273 delete fGeo;
274 }
275 fGeo = new AliTRDgeometry();
55a288e5 276}
f162af62 277
55a288e5 278//____________________________________________________________________________________
279AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
280{
281 //
282 // AliTRDCalibraFillHisto destructor
283 //
284
285 ClearHistos();
170c35f1 286 if ( fDebugStreamer ) delete fDebugStreamer;
f162af62 287
d0569428 288 if ( fCalDetGain ) delete fCalDetGain;
289 if ( fCalDetT0 ) delete fCalDetT0;
290 if ( fCalROCGain ) delete fCalROCGain;
291 if ( fCalROCT0 ) delete fCalROCT0;
292
f162af62 293 if (fGeo) {
294 delete fGeo;
295 }
55a288e5 296
297}
298
299//_____________________________________________________________________________
300void AliTRDCalibraFillHisto::Destroy()
301{
302 //
303 // Delete instance
304 //
305
306 if (fgInstance) {
307 delete fgInstance;
308 fgInstance = 0x0;
309 }
310
311}
312
313//_____________________________________________________________________________
314void AliTRDCalibraFillHisto::ClearHistos()
315{
316 //
317 // Delete the histos
318 //
319
320 if (fPH2d) {
321 delete fPH2d;
322 fPH2d = 0x0;
323 }
324 if (fCH2d) {
325 delete fCH2d;
326 fCH2d = 0x0;
327 }
328 if (fPRF2d) {
329 delete fPRF2d;
330 fPRF2d = 0x0;
331 }
3a0f6479 332
55a288e5 333}
3a0f6479 334
55a288e5 335//____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
336Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
337{
338 //
339 // For the offline tracking
340 // This function will be called in the function AliReconstruction::Run()
341 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
342 //
343
d0569428 344 Init2Dhistostrack();
345
346 //Init the tracklet parameters
347 fPar0 = new Double_t[fTimeMax];
348 fPar1 = new Double_t[fTimeMax];
349 fPar2 = new Double_t[fTimeMax];
350 fPar3 = new Double_t[fTimeMax];
351 fPar4 = new Double_t[fTimeMax];
352
353 for(Int_t k = 0; k < fTimeMax; k++){
354 fPar0[k] = 0.0;
355 fPar1[k] = 0.0;
356 fPar2[k] = 0.0;
357 fPar3[k] = 0.0;
358 fPar4[k] = 0.0;
359 }
360 return kTRUE;
361}
362
363//____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
364Bool_t AliTRDCalibraFillHisto::Init2Dhistostrack()
365{
366 //
367 // For the offline tracking
368 // This function will be called in the function AliReconstruction::Run()
369 // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
370 //
371
55a288e5 372 // DB Setting
373 // Get cal
374 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
375 if (!cal) {
376 AliInfo("Could not get calibDB");
377 return kFALSE;
378 }
379 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
380 if (!parCom) {
381 AliInfo("Could not get CommonParam");
382 return kFALSE;
383 }
384
385 // Some parameters
386 fTimeMax = cal->GetNumberOfTimeBins();
387 fSf = parCom->GetSamplingFrequency();
170c35f1 388 fRelativeScale = 20;
d0569428 389
390 //calib object from database used for reconstruction
391 if(fCalDetGain) delete fCalDetGain;
392 fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
393 if(fCalDetT0) delete fCalDetT0;
394 fCalDetT0 = new AliTRDCalDet(*(cal->GetT0Det()));
170c35f1 395
3a0f6479 396 // Calcul Xbins Chambd0, Chamb2
0c349049 397 Int_t ntotal0 = CalculateTotalNumberOfBins(0);
398 Int_t ntotal1 = CalculateTotalNumberOfBins(1);
399 Int_t ntotal2 = CalculateTotalNumberOfBins(2);
3a0f6479 400
401 // If vector method On initialised all the stuff
402 if(fVector2d){
55a288e5 403 fCalibraVector = new AliTRDCalibraVector();
170c35f1 404 fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
405 fCalibraVector->SetTimeMax(fTimeMax);
406 if(fNgroupprf != 0) {
407 fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
408 fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
409 }
410 else {
411 fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
412 fCalibraVector->SetPRFRange(1.5);
413 }
3a0f6479 414 for(Int_t k = 0; k < 3; k++){
415 fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
416 fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
417 }
418 TString namech("Nz");
419 namech += fCalibraMode->GetNz(0);
420 namech += "Nrphi";
421 namech += fCalibraMode->GetNrphi(0);
422 fCalibraVector->SetNameCH((const char* ) namech);
423 TString nameph("Nz");
424 nameph += fCalibraMode->GetNz(1);
425 nameph += "Nrphi";
426 nameph += fCalibraMode->GetNrphi(1);
427 fCalibraVector->SetNamePH((const char* ) nameph);
428 TString nameprf("Nz");
429 nameprf += fCalibraMode->GetNz(2);
430 nameprf += "Nrphi";
431 nameprf += fCalibraMode->GetNrphi(2);
432 nameprf += "Ngp";
433 nameprf += fNgroupprf;
434 fCalibraVector->SetNamePRF((const char* ) nameprf);
55a288e5 435 }
170c35f1 436
55a288e5 437 // Create the 2D histos corresponding to the pad groupCalibration mode
438 if (fCH2dOn) {
439
440 AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
441 ,fCalibraMode->GetNz(0)
442 ,fCalibraMode->GetNrphi(0)));
443
55a288e5 444 // Create the 2D histo
445 if (fHisto2d) {
0c349049 446 CreateCH2d(ntotal0);
55a288e5 447 }
55a288e5 448 // Variable
449 fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
450 for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
451 fAmpTotal[k] = 0.0;
452 }
170c35f1 453 //Statistics
0c349049 454 fEntriesCH = new Int_t[ntotal0];
455 for(Int_t k = 0; k < ntotal0; k++){
170c35f1 456 fEntriesCH[k] = 0;
457 }
3a0f6479 458
55a288e5 459 }
55a288e5 460 if (fPH2dOn) {
461
462 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
463 ,fCalibraMode->GetNz(1)
464 ,fCalibraMode->GetNrphi(1)));
465
55a288e5 466 // Create the 2D histo
467 if (fHisto2d) {
0c349049 468 CreatePH2d(ntotal1);
55a288e5 469 }
55a288e5 470 // Variable
471 fPHPlace = new Short_t[fTimeMax];
472 for (Int_t k = 0; k < fTimeMax; k++) {
473 fPHPlace[k] = -1;
474 }
475 fPHValue = new Float_t[fTimeMax];
476 for (Int_t k = 0; k < fTimeMax; k++) {
477 fPHValue[k] = 0.0;
478 }
170c35f1 479 }
480 if (fLinearFitterOn) {
d0569428 481 //fLinearFitterArray.Expand(540);
170c35f1 482 fLinearFitterArray.SetName("ArrayLinearFitters");
483 fEntriesLinearFitter = new Int_t[540];
484 for(Int_t k = 0; k < 540; k++){
485 fEntriesLinearFitter[k] = 0;
486 }
3a0f6479 487 fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
55a288e5 488 }
489
490 if (fPRF2dOn) {
491
492 AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
493 ,fCalibraMode->GetNz(2)
494 ,fCalibraMode->GetNrphi(2)));
55a288e5 495 // Create the 2D histo
496 if (fHisto2d) {
0c349049 497 CreatePRF2d(ntotal2);
55a288e5 498 }
55a288e5 499 }
500
501 return kTRUE;
502
503}
504
505//____________Functions for filling the histos in the code_____________________
506
507//____________Offine tracking in the AliTRDtracker_____________________________
508Bool_t AliTRDCalibraFillHisto::ResetTrack()
509{
510 //
511 // For the offline tracking
512 // This function will be called in the function
513 // AliTRDtracker::FollowBackPropagation() at the beginning
514 // Reset the parameter to know we have a new TRD track
515 //
516
517 fDetectorAliTRDtrack = kFALSE;
d0569428 518 //fGoodTrack = kTRUE;
55a288e5 519 return kTRUE;
520
521}
170c35f1 522//____________Offine tracking in the AliTRDtracker_____________________________
523void AliTRDCalibraFillHisto::ResetfVariables()
524{
525 //
526 // Reset values per tracklet
527 //
528
529 // Reset the list of clusters
530 fListClusters->Clear();
531
532 //Reset the tracklet parameters
533 for(Int_t k = 0; k < fTimeMax; k++){
534 fPar0[k] = 0.0;
535 fPar1[k] = 0.0;
536 fPar2[k] = 0.0;
537 fPar3[k] = 0.0;
538 fPar4[k] = 0.0;
539 }
540
d0569428 541 ResetfVariablestrack();
542
543}
544//____________Offine tracking in the AliTRDtracker_____________________________
545void AliTRDCalibraFillHisto::ResetfVariablestrack()
546{
547 //
548 // Reset values per tracklet
549 //
550
170c35f1 551 //Reset good tracklet
552 fGoodTracklet = kTRUE;
55a288e5 553
170c35f1 554 // Reset the fPHValue
555 if (fPH2dOn) {
556 //Reset the fPHValue and fPHPlace
557 for (Int_t k = 0; k < fTimeMax; k++) {
558 fPHValue[k] = 0.0;
559 fPHPlace[k] = -1;
560 }
561 }
562
563 // Reset the fAmpTotal where we put value
564 if (fCH2dOn) {
565 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
566 fAmpTotal[k] = 0.0;
567 }
568 }
d0569428 569}
570//____________Offline tracking in the AliTRDtracker____________________________
571Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
572{
573 //
574 // For the offline tracking
575 // This function will be called in the function
576 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
577 // of TRD tracks
578 // Fill the 2D histos or the vectors with the info of the clusters at
579 // the end of a detectors if the track is "good"
580 //
581
582
583
584 AliTRDcluster *cl = 0x0;
585 Int_t index0 = 0;
586 Int_t index1 = 0;
587
588 // reset if good track
589 fGoodTrack = kTRUE;
590
591
592 // loop over the clusters
593 while((cl = t->GetCluster(index1))){
594
595 // Localisation of the detector
596 Int_t detector = cl->GetDetector();
597
598
599 // Fill the infos for the previous clusters if not the same
600 // detector anymore but this time it should be the same track
601 if ((detector != fDetectorPreviousTrack) &&
602 (index0 != index1)) {
603
604 fNumberTrack++;
605
606 //printf("detector %d, fPreviousdetector %d, plane %d, planeprevious %d, index0 %d, index1 %d la\n",detector,fDetectorPreviousTrack,GetPlane(detector),GetPlane(fDetectorPreviousTrack),index0,index1);
607
608 //If the same track, then look if the previous detector is in
609 //the same plane, if yes: not a good track
610 //FollowBack
611 //if (fDetectorAliTRDtrack &&
612 // (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
613 //Follow
614 if ((GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) {
615 fGoodTrack = kFALSE;
616 }
617
618 // Fill only if the track doesn't touch a masked pad or doesn't
619 // appear in the middle (fGoodTrack)
620 if (fGoodTrack && fGoodTracklet) {
621
622 // drift velocity unables to cut bad tracklets
623 Bool_t pass = FindP1TrackPHtrack(t,index0,index1);
624
625 // Gain calibration
626 if (fCH2dOn) {
627 FillTheInfoOfTheTrackCH();
628 }
629
630 // PH calibration
631 if (fPH2dOn) {
632 FillTheInfoOfTheTrackPH();
633 }
634
635 if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1);
636
637
638 } // if a good track
639
640 // reset stuff
641 ResetfVariablestrack();
642 index0 = index1;
643
644 } // Fill at the end the charge
645
646 // Calcul the position of the detector and take the calib objects
647 if (detector != fDetectorPreviousTrack) {
648
649 //Localise the detector
650 LocalisationDetectorXbins(detector);
651
652 // Get cal
653 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
654 if (!cal) {
655 AliInfo("Could not get calibDB");
656 return kFALSE;
657 }
658
659 // Get calib objects
660 if( fCalROCGain ) delete fCalROCGain;
661 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
662 if( fCalROCT0 ) delete fCalROCT0;
663 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector)));
664
665 }
666
667 // Reset the detectbjobsor
668 fDetectorPreviousTrack = detector;
669
670 // Store the info bis of the tracklet
671 Int_t *rowcol = CalculateRowCol(cl);
672 CheckGoodTracklet(detector,rowcol);
673 Int_t group[2] = {0,0};
674 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
675 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
676 StoreInfoCHPHtrack(cl,t,index1,group,rowcol);
677
678 index1++;
679
680 } // while on clusters
681
682 // Fill the last plane
683 if( index0 != index1 ){
684
685 //printf("fPreviousdetector %d, planeprevious %d, index0 %d, index1 %d li\n",fDetectorPreviousTrack,GetPlane(fDetectorPreviousTrack),index0,index1);
686
687 fNumberTrack++;
688
689 if (fGoodTrack && fGoodTracklet) {
690
691 // drift velocity unables to cut bad tracklets
692 Bool_t pass = FindP1TrackPHtrack(t,index0,index1);
693
694 // Gain calibration
695 if (fCH2dOn) {
696 FillTheInfoOfTheTrackCH();
697 }
698
699 // PH calibration
700 if (fPH2dOn) {
701 FillTheInfoOfTheTrackPH();
702 }
703
704 if(pass && fPRF2dOn) HandlePRFtrack(t,index0,index1);
705
706 } // if a good track
707
708 }
170c35f1 709
d0569428 710 // reset stuff
711 ResetfVariablestrack();
712
713 return kTRUE;
714
170c35f1 715}
55a288e5 716//____________Offline tracking in the AliTRDtracker____________________________
717Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
718{
719 //
720 // For the offline tracking
721 // This function will be called in the function
722 // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
723 // of TRD tracks
724 // Fill the 2D histos or the vectors with the info of the clusters at
725 // the end of a detectors if the track is "good"
726 //
727
55a288e5 728 // Localisation of the detector
729 Int_t detector = cl->GetDetector();
3a0f6479 730
55a288e5 731 // Fill the infos for the previous clusters if not the same
732 // detector anymore or if not the same track
733 if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
734 (fDetectorPreviousTrack != -1)) {
735
736 fNumberTrack++;
737
738 // If the same track, then look if the previous detector is in
739 // the same plane, if yes: not a good track
3a0f6479 740 //FollowBack
55a288e5 741 if (fDetectorAliTRDtrack &&
3a0f6479 742 (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
743 //Follow
744 //if (fDetectorAliTRDtrack &&
745 // (GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) {
55a288e5 746 fGoodTrack = kFALSE;
747 }
748
749 // Fill only if the track doesn't touch a masked pad or doesn't
750 // appear in the middle (fGoodTrack)
170c35f1 751 if (fGoodTrack && fGoodTracklet) {
55a288e5 752
170c35f1 753 // drift velocity unables to cut bad tracklets
754 Bool_t pass = FindP1TrackPH();
755
55a288e5 756 // Gain calibration
757 if (fCH2dOn) {
758 FillTheInfoOfTheTrackCH();
759 }
760
761 // PH calibration
762 if (fPH2dOn) {
763 FillTheInfoOfTheTrackPH();
764 }
170c35f1 765
766 if(pass && fPRF2dOn) HandlePRF();
767
55a288e5 768
769 } // if a good track
770
771 ResetfVariables();
d0569428 772 if(!fDetectorAliTRDtrack) fGoodTrack = kTRUE;
55a288e5 773
774 } // Fill at the end the charge
775
d0569428 776 // Calcul the position of the detector and take the calib objects
55a288e5 777 if (detector != fDetectorPreviousTrack) {
d0569428 778 //Localise the detector
55a288e5 779 LocalisationDetectorXbins(detector);
d0569428 780
781 // Get cal
782 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
783 if (!cal) {
784 AliInfo("Could not get calibDB");
785 return kFALSE;
786 }
787
788 // Get calib objects
789 if( fCalROCGain ) delete fCalROCGain;
790 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
791 if( fCalROCT0 ) delete fCalROCT0;
792 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(detector)));
55a288e5 793 }
170c35f1 794
55a288e5 795 // Reset the detector
796 fDetectorPreviousTrack = detector;
797 fDetectorAliTRDtrack = kTRUE;
55a288e5 798
170c35f1 799 // Store the infos of the tracklets
800 AliTRDcluster *kcl = new AliTRDcluster(*cl);
801 fListClusters->Add((TObject *)kcl);
802 Int_t time = cl->GetLocalTimeBin();
803 fPar0[time] = t->GetY();
804 fPar1[time] = t->GetZ();
805 fPar2[time] = t->GetSnp();
806 fPar3[time] = t->GetTgl();
6c23ffed 807 fPar4[time] = t->GetSigned1Pt();
170c35f1 808
809 // Store the info bis of the tracklet
810 Int_t *rowcol = CalculateRowCol(cl);
811 CheckGoodTracklet(detector,rowcol);
812 Int_t group[2] = {0,0};
813 if(fCH2dOn) group[0] = CalculateCalibrationGroup(0,rowcol);
814 if(fPH2dOn) group[1] = CalculateCalibrationGroup(1,rowcol);
d0569428 815 StoreInfoCHPH(cl,t,group,rowcol);
170c35f1 816
55a288e5 817 return kTRUE;
818
819}
55a288e5 820//____________Online trackling in AliTRDtrigger________________________________
821Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
822{
823 //
824 // For the tracking
825 // This function will be called in the function AliTRDtrigger::TestTracklet
826 // before applying the pt cut on the tracklets
827 // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
828 //
829
830 // Localisation of the Xbins involved
831 Int_t idect = trk->GetDetector();
3a0f6479 832 fDetectorPreviousTrack = idect;
55a288e5 833 LocalisationDetectorXbins(idect);
834
835 // Get the parameter object
836 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
837 if (!cal) {
838 AliInfo("Could not get calibDB");
839 return kFALSE;
840 }
841
842 // Reset
843 ResetfVariables();
d0569428 844
845 // Get calib objects
846 if( fCalROCGain ) delete fCalROCGain;
847 fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(idect)));
848 if( fCalROCT0 ) delete fCalROCT0;
849 fCalROCT0 = new AliTRDCalROC(*(cal->GetT0ROC(idect)));
850
55a288e5 851 // Row of the tracklet and position in the pad groups
170c35f1 852 Int_t *rowcol = new Int_t[2];
853 rowcol[0] = trk->GetRow();
854 Int_t group[3] = {-1,-1,-1};
855
55a288e5 856 // Eventuelle correction due to track angle in z direction
857 Float_t correction = 1.0;
858 if (fMcmCorrectAngle) {
859 Float_t z = trk->GetRowz();
860 Float_t r = trk->GetTime0();
861 correction = r / TMath::Sqrt((r*r+z*z));
862 }
863
864 // Boucle sur les clusters
865 // Condition on number of cluster: don't come from the middle of the detector
866 if (trk->GetNclusters() >= fNumberClusters) {
867
868 for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
869
870 Float_t amp[3] = { 0.0, 0.0, 0.0 };
871 Int_t time = trk->GetClusterTime(icl);
170c35f1 872 rowcol[1] = trk->GetClusterCol(icl);
55a288e5 873
874 amp[0] = trk->GetClusterADC(icl)[0] * correction;
875 amp[1] = trk->GetClusterADC(icl)[1] * correction;
876 amp[2] = trk->GetClusterADC(icl)[2] * correction;
877
878
879 if ((amp[0] < 0.0) ||
880 (amp[1] < 0.0) ||
881 (amp[2] < 0.0)) {
882 continue;
883 }
884
885 // Col of cluster and position in the pad groups
170c35f1 886 if(fCH2dOn) {
887 group[0] = CalculateCalibrationGroup(0,rowcol);
888 fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
55a288e5 889 }
170c35f1 890 if(fPH2dOn) {
891 group[1] = CalculateCalibrationGroup(1,rowcol);
892 fPHPlace[time] = group[1];
55a288e5 893 fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
894 }
170c35f1 895 if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
896
897 // See if we are not near a masked pad fGoodTracklet
898 CheckGoodTracklet(idect,rowcol);
899
900 // Fill PRF direct without tnp bins...only for monitoring...
901 if (fPRF2dOn && fGoodTracklet) {
55a288e5 902
903 if ((amp[0] > fThresholdClusterPRF2) &&
904 (amp[1] > fThresholdClusterPRF2) &&
905 (amp[2] > fThresholdClusterPRF2) &&
906 ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
907
908 // Security of the denomiateur is 0
909 if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
910 / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
911 Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
912 / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
913 Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
914
915 if (TMath::Abs(xcenter) < 0.5) {
916 Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
917 Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
918 // Fill only if it is in the drift region!
919 if (((Float_t) time / fSf) > 0.3) {
920 if (fHisto2d) {
170c35f1 921 fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
922 fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
923 fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
55a288e5 924 }
925 if (fVector2d) {
3a0f6479 926 fCalibraVector->UpdateVectorPRF(idect,group[2],xcenter,ycenter);
927 fCalibraVector->UpdateVectorPRF(idect,group[2],-(xcenter+1.0),yminus);
928 fCalibraVector->UpdateVectorPRF(idect,group[2],(1.0-xcenter),ymax);
55a288e5 929 }
930 }//in the drift region
931 }//in the middle
932 }//denominateur security
933 }//cluster shape and thresholds
934 }//good and PRF On
935
936 } // Boucle clusters
937
938 // Fill the charge
170c35f1 939 if(fGoodTracklet){
940 if (fCH2dOn) FillTheInfoOfTheTrackCH();
941 if (fPH2dOn) FillTheInfoOfTheTrackPH();
55a288e5 942 }
943
944 fNumberTrack++;
945
946 } // Condition on number of clusters
947
948 return kTRUE;
949
950}
55a288e5 951//_____________________________________________________________________________
170c35f1 952Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
55a288e5 953{
954 //
170c35f1 955 // Calculate the row and col number of the cluster
55a288e5 956 //
957
170c35f1 958
959 Int_t *rowcol = new Int_t[2];
960 rowcol[0] = 0;
961 rowcol[1] = 0;
962
170c35f1 963 // Localisation of the detector
964 Int_t detector = cl->GetDetector();
965 Int_t chamber = GetChamber(detector);
966 Int_t plane = GetPlane(detector);
967
968 // Localisation of the cluster
969 Double_t pos[3] = { 0.0, 0.0, 0.0 };
970 pos[0] = ((AliCluster *)cl)->GetX();
971 pos[1] = cl->GetY();
972 pos[2] = cl->GetZ();
973
974 // Position of the cluster
f162af62 975 AliTRDpadPlane *padplane = fGeo->GetPadPlane(plane,chamber);
170c35f1 976 Int_t row = padplane->GetPadRowNumber(pos[2]);
977 //Do not take from here because it was corrected from ExB already....
978 //Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
979 //Double_t offsettilt = padplane->GetTiltOffset(offsetz);
980 //Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
981 //Int_t col = padplane->GetPadColNumber(pos[1]+offsettilt);
af26ce80 982 Int_t col = cl->GetPadCol();
170c35f1 983
984 //return
985 rowcol[0] = row;
986 rowcol[1] = col;
987 return rowcol;
988
989}
990//_____________________________________________________________________________
991void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
992{
993 //
994 // See if we are not near a masked pad
995 //
996
997 Int_t row = rowcol[0];
998 Int_t col = rowcol[1];
999
1000 if (!IsPadOn(detector, col, row)) {
1001 fGoodTracklet = kFALSE;
1002 }
1003
1004 if (col > 0) {
1005 if (!IsPadOn(detector, col-1, row)) {
1006 fGoodTracklet = kFALSE;
1007 }
1008 }
1009
1010 if (col < 143) {
1011 if (!IsPadOn(detector, col+1, row)) {
1012 fGoodTracklet = kFALSE;
1013 }
1014 }
1015
1016}
1017//_____________________________________________________________________________
1018Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
1019{
1020 //
1021 // Look in the choosen database if the pad is On.
1022 // If no the track will be "not good"
1023 //
1024
1025 // Get the parameter object
1026 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1027 if (!cal) {
1028 AliInfo("Could not get calibDB");
1029 return kFALSE;
1030 }
1031
1032 if (!cal->IsChamberInstalled(detector) ||
55a288e5 1033 cal->IsChamberMasked(detector) ||
1034 cal->IsPadMasked(detector,col,row)) {
1035 return kFALSE;
1036 }
1037 else {
1038 return kTRUE;
1039 }
1040
1041}
170c35f1 1042//_____________________________________________________________________________
1043Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
1044{
1045 //
1046 // Calculate the calibration group number for i
1047 //
1048
1049 // Row of the cluster and position in the pad groups
1050 Int_t posr = 0;
1051 if (fCalibraMode->GetNnZ(i) != 0) {
1052 posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
1053 }
1054
1055
1056 // Col of the cluster and position in the pad groups
1057 Int_t posc = 0;
1058 if (fCalibraMode->GetNnRphi(i) != 0) {
1059 posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
1060 }
1061
1062 return posc*fCalibraMode->GetNfragZ(i)+posr;
1063
1064}
1065//____________________________________________________________________________________
1066Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1067{
1068 //
1069 // Calculate the total number of calibration groups
1070 //
1071
0c349049 1072 Int_t ntotal = 0;
170c35f1 1073 fCalibraMode->ModePadCalibration(2,i);
1074 fCalibraMode->ModePadFragmentation(0,2,0,i);
1075 fCalibraMode->SetDetChamb2(i);
0c349049 1076 ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
170c35f1 1077 fCalibraMode->ModePadCalibration(0,i);
1078 fCalibraMode->ModePadFragmentation(0,0,0,i);
1079 fCalibraMode->SetDetChamb0(i);
0c349049 1080 ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1081 AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1082 return ntotal;
170c35f1 1083
1084}
1085//____________Set the pad calibration variables for the detector_______________
1086Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1087{
1088 //
1089 // For the detector calcul the first Xbins and set the number of row
1090 // and col pads per calibration groups, the number of calibration
1091 // groups in the detector.
1092 //
1093
1094 // first Xbins of the detector
1095 if (fCH2dOn) {
1096 fCalibraMode->CalculXBins(detector,0);
1097 }
1098 if (fPH2dOn) {
1099 fCalibraMode->CalculXBins(detector,1);
1100 }
1101 if (fPRF2dOn) {
1102 fCalibraMode->CalculXBins(detector,2);
1103 }
1104
1105 // fragmentation of idect
1106 for (Int_t i = 0; i < 3; i++) {
1107 fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1108 fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1109 , (Int_t) GetChamber(detector)
1110 , (Int_t) GetSector(detector),i);
1111 }
1112
1113 return kTRUE;
55a288e5 1114
170c35f1 1115}
1116//_____________________________________________________________________________
d0569428 1117void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group, Int_t *rowcol)
170c35f1 1118{
1119 //
1120 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1121 //
1122
1123 // Charge in the cluster
1124 Float_t q = TMath::Abs(cl->GetQ());
1125 Int_t time = cl->GetLocalTimeBin();
d0569428 1126
1127 //Correct for the gain coefficient used in the database for reconstruction
1128 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]);
1129 Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]);
1130
1131 // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458)
1132 Int_t timec = Arrondi((Double_t)(time+correcttheT0));
1133 if((correcttheT0+0.5)==(int(correcttheT0+0.5))) {
1134 timec++;
1135 }
1136 if( timec < 0 ) return;
1137
1138
170c35f1 1139 // Correction due to the track angle
1140 Float_t correction = 1.0;
1141 Float_t normalisation = 6.67;
d0569428 1142 // we divide with gain in AliTRDclusterizerV1::Transform...
1143 if( correctthegain > 0 ) normalisation /= correctthegain;
170c35f1 1144 if ((q >0) && (t->GetNdedx() > 0)) {
1145 correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
1146 }
1147
1148 // Fill the fAmpTotal with the charge
1149 if (fCH2dOn) {
1150 fAmpTotal[(Int_t) group[0]] += correction;
1151 }
1152
1153 // Fill the fPHPlace and value
1154 if (fPH2dOn) {
d0569428 1155 fPHPlace[timec] = group[1];
1156 fPHValue[timec] = correction;
1157 }
1158
1159}
1160//_____________________________________________________________________________
1161void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, AliTRDtrack *t, Int_t index, Int_t *group, Int_t *rowcol)
1162{
1163 //
1164 // Store the infos in fAmpTotal, fPHPlace and fPHValue
1165 //
1166
1167 // Charge in the cluster
1168 Float_t q = TMath::Abs(cl->GetQ());
1169 Int_t time = cl->GetLocalTimeBin();
1170
1171 //Correct for the gain coefficient used in the database for reconstruction
1172 Float_t correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(rowcol[1],rowcol[0]);
1173 Float_t correcttheT0 = fCalDetT0->GetValue(fDetectorPreviousTrack)+fCalROCT0->GetValue(rowcol[1],rowcol[0]);
1174
1175 // we substract correcttheT0 in AliTRDclusterizerV1::MakeClusters (line 458)
1176 Int_t timec = Arrondi((Double_t)(time+correcttheT0));
1177 if((correcttheT0+0.5)==(int(correcttheT0+0.5))) {
1178 timec++;
1179 }
1180 if( timec < 0 ) return;
1181
1182 // Correction due to the track angle
1183 Float_t correction = 1.0;
1184 Float_t normalisation = 6.67;
1185 // we divide with gain in AliTRDclusterizerV1::Transform...
1186 if( correctthegain > 0 ) normalisation /= correctthegain;
1187 if (q >0) {
1188 correction = t->GetClusterdQdl(index) / (normalisation);
1189 }
1190
1191 // Fill the fAmpTotal with the charge
1192 if (fCH2dOn) {
1193 fAmpTotal[(Int_t) group[0]] += correction;
1194 }
1195
1196 // Fill the fPHPlace and value
1197 if (fPH2dOn) {
1198 fPHPlace[timec] = group[1];
1199 fPHValue[timec] = correction;
170c35f1 1200 }
1201
1202}
1203//_____________________________________________________________________
d4232bb6 1204Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamTB *rawStream, Bool_t nocheck)
170c35f1 1205{
1206 //
d4232bb6 1207 // Event Processing loop - AliTRDrawStreamTB
3a0f6479 1208 // 0 timebin problem
1209 // 1 no input
1210 // 2 input
170c35f1 1211 //
e4db522f 1212
3a0f6479 1213 Int_t withInput = 1;
e4db522f 1214
1215 Int_t phvalue[21][36];
170c35f1 1216 for(Int_t k = 0; k < 36; k++){
e4db522f 1217 for(Int_t j = 0; j < 21; j++){
1218 phvalue[j][k] = 10;
1219 }
170c35f1 1220 }
e4db522f 1221
170c35f1 1222 fDetectorPreviousTrack = -1;
e4db522f 1223 fMCMPrevious = -1;
1224 fROBPrevious = -1;
1225 Int_t nbtimebin = 0;
1226 Int_t baseline = 10;
1227
1228 // For selecting the signal
1229 Double_t mean[21] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1230 Int_t first[21] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1231 Int_t select[21] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
170c35f1 1232
3a0f6479 1233 if(!nocheck){
1234
1235 fTimeMax = 0;
e4db522f 1236
3a0f6479 1237 while (rawStream->Next()) {
1238
1239 Int_t idetector = rawStream->GetDet(); // current detector
e4db522f 1240 Int_t imcm = rawStream->GetMCM(); // current MCM
1241 Int_t irob = rawStream->GetROB(); // current ROB
1242
1243 if(((fMCMPrevious != imcm) || (fDetectorPreviousTrack != idetector) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
1244
1245 // take the mean values and check the first time bin
162a3e73 1246 for(Int_t j = 0; j < 21; j++){
e4db522f 1247 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1248 else mean[j] = 0.0;
1249 if(phvalue[j][0] > 200.0) first[j] = 1;
1250 else first[j] = 0;
1251 }
1252
1253 // select
1254 for(Int_t j = 1; j < 20; j++){
1255 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1256 select[j] = 1;
1257 }
1258 else select[j] = 0;
1259 }
1260
1261 // fill
1262 for(Int_t j = 1; j < 20; j++){
1263 if(select[j] == 1){
1264 withInput = 2;
1265 for(Int_t k = 0; k < fTimeMax; k++){
1266 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1267 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1268 }
1269 else{
1270 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1271 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1272 }
1273 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,(3*baseline),fTimeMax);
1274 }
1275 }
1276 }
1277 }
1278
1279 // reset
162a3e73 1280 for(Int_t j = 0; j < 21; j++){
1281 mean[j] = 0.0;
1282 first[j] = 0;
1283 select[j] = 0;
1284 for(Int_t k = 0; k < 36; k++){
e4db522f 1285 phvalue[j][k] = baseline;
3a0f6479 1286 }
1287 }
1288 }
e4db522f 1289
3a0f6479 1290 fDetectorPreviousTrack = idetector;
e4db522f 1291 fMCMPrevious = imcm;
1292 fROBPrevious = irob;
1293
3a0f6479 1294 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1295 if(nbtimebin == 0) return 0;
1296 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
1297 fTimeMax = nbtimebin;
e4db522f 1298
162a3e73 1299 baseline = rawStream->GetCommonAdditive(); // common additive baseline
e4db522f 1300
3a0f6479 1301 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
162a3e73 1302 Int_t *signal = rawStream->GetSignals(); // current ADC signal
e4db522f 1303 Int_t col = (rawStream->GetCol())%18; // current COL MCM
162a3e73 1304
1305 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
1306
e4db522f 1307 if((col < 0) || (col >= 21)) return 0;
1308 if((imcm>=16) || (imcm < 0)) return 0;
1309
3a0f6479 1310 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
1311 Int_t n = 0;
1312 for(Int_t itime = iTimeBin; itime < fin; itime++){
162a3e73 1313 phvalue[col][itime] = signal[n];
3a0f6479 1314 n++;
1315 }
1316 }
e4db522f 1317
3a0f6479 1318 // fill the last one
1319 if(fDetectorPreviousTrack != -1){
e4db522f 1320
1321 // take the mean values and check the first time bin
1322 for(Int_t j = 0; j < 21; j++){
1323 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1324 else mean[j] = 0.0;
1325 if(phvalue[j][0] > 200.0) first[j] = 1;
1326 else first[j] = 0;
1327 }
1328
1329 // select
1330 for(Int_t j = 1; j < 20; j++){
1331 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1332 select[j] = 1;
1333 }
1334 else select[j] = 0;
1335 }
1336
1337 // fill
1338 for(Int_t j = 1; j < 20; j++){
1339 if(select[j] == 1){
1340 withInput = 2;
1341 for(Int_t k = 0; k < fTimeMax; k++){
1342 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1343 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1344 }
1345 else{
1346 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1347 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1348 }
1349 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,(3*baseline),fTimeMax);
1350 }
1351 }
1352 }
1353 }
1354
1355 // reset
162a3e73 1356 for(Int_t j = 0; j < 21; j++){
1357 mean[j] = 0.0;
1358 first[j] = 0;
1359 select[j] = 0;
1360 for(Int_t k = 0; k < 36; k++){
e4db522f 1361 phvalue[j][k] = baseline;
170c35f1 1362 }
1363 }
1364 }
3a0f6479 1365
170c35f1 1366 }
3a0f6479 1367 else{
1368
1369 while (rawStream->Next()) {
1370
1371 Int_t idetector = rawStream->GetDet(); // current detector
e4db522f 1372 Int_t imcm = rawStream->GetMCM(); // current MCM
1373 Int_t irob = rawStream->GetROB(); // current ROB
1374
1375 if(((fMCMPrevious != imcm) || (fDetectorPreviousTrack != idetector) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
1376
1377 // take the mean values and check the first time bin
1378 for(Int_t j = 0; j < 21; j++){
1379 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1380 else mean[j] = 0.0;
1381 if(phvalue[j][0] > 200.0) first[j] = 1;
1382 else first[j] = 0;
162a3e73 1383 //printf("meanvalue %f, first %d\n",mean[j],first[j]);
e4db522f 1384 }
1385
1386 // select
1387 for(Int_t j = 1; j < 20; j++){
1388 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1389 select[j] = 1;
1390 }
1391 else select[j] = 0;
1392 }
1393
1394 // fill
1395 for(Int_t j = 1; j < 20; j++){
162a3e73 1396 //printf("select %d\n",select[j]);
1397 if(select[j] == 1){
e4db522f 1398 withInput = 2;
1399 for(Int_t k = 0; k < fTimeMax; k++){
1400 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
162a3e73 1401 printf("fill\n");
e4db522f 1402 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1403 }
1404 else{
1405 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
162a3e73 1406 printf("fill\n");
e4db522f 1407 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1408 }
1409 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1410 }
1411 }
1412 }
1413 }
1414
1415 // reset
162a3e73 1416 for(Int_t j = 0; j < 21; j++){
1417 mean[j] = 0.0;
1418 first[j] = 0;
1419 select[j] = 0;
1420 for(Int_t k = 0; k < 36; k++){
e4db522f 1421 phvalue[j][k] = baseline;
3a0f6479 1422 }
1423 }
1424 }
e4db522f 1425
3a0f6479 1426 fDetectorPreviousTrack = idetector;
e4db522f 1427 fMCMPrevious = imcm;
1428 fROBPrevious = irob;
1429
1430
1431
162a3e73 1432 baseline = rawStream->GetCommonAdditive(); // common baseline
e4db522f 1433
3a0f6479 1434 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1435 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
162a3e73 1436 Int_t *signal = rawStream->GetSignals(); // current ADC signal
e4db522f 1437 Int_t col = (rawStream->GetCol())%18; // current COL MCM
1438
3a0f6479 1439 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
1440 Int_t n = 0;
162a3e73 1441
1442 //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
1443
1444
e4db522f 1445 if((col < 0) || (col >= 21)) return 0;
1446 if((imcm>=16) || (imcm < 0)) return 0;
1447
3a0f6479 1448 for(Int_t itime = iTimeBin; itime < fin; itime++){
162a3e73 1449 phvalue[col][itime] = signal[n];
3a0f6479 1450 n++;
1451 }
1452 }
1453
1454 // fill the last one
1455 if(fDetectorPreviousTrack != -1){
e4db522f 1456
1457 // take the mean values and check the first time bin
1458 for(Int_t j = 0; j < 21; j++){
1459 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1460 else mean[j] = 0.0;
1461 if(phvalue[j][0] > 200.0) first[j] = 1;
1462 else first[j] = 0;
1463 }
1464
1465 // select
1466 for(Int_t j = 1; j < 20; j++){
1467 if((first[j-1] == 0) && (first[j] ==0) && (first[j+1] == 0) && (mean[j-1] > (baseline+5.0)) && (mean[j] > (baseline+10.0)) && (mean[j+1] > (baseline+5.0)) && (mean[j] >= mean[j-1]) && (mean[j] >= mean[j+1])){
1468 select[j] = 1;
1469 }
1470 else select[j] = 0;
1471 }
1472
1473 // fill
1474 for(Int_t j = 1; j < 20; j++){
1475 if(select[j] == 1){
1476 withInput = 2;
1477 for(Int_t k = 0; k < fTimeMax; k++){
1478 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1479 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1480 }
1481 else{
1482 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1483 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1484 }
1485 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1486 }
1487 }
1488 }
1489 }
1490
1491 // reset
162a3e73 1492 for(Int_t j = 0; j < 21; j++){
1493 mean[j] = 0.0;
1494 first[j] = 0;
1495 select[j] = 0;
1496 for(Int_t k = 0; k < 36; k++){
e4db522f 1497 phvalue[j][k] = baseline;
170c35f1 1498 }
1499 }
3a0f6479 1500 }
170c35f1 1501 }
1502
1503 return withInput;
e4db522f 1504
170c35f1 1505}
1506//_____________________________________________________________________
3a0f6479 1507Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
170c35f1 1508{
1509 //
1510 // Event processing loop - AliRawReader
1511 //
1512
1513
d4232bb6 1514 AliTRDrawStreamTB rawStream(rawReader);
170c35f1 1515
1516 rawReader->Select("TRD");
1517
3a0f6479 1518 return ProcessEventDAQ(&rawStream, nocheck);
170c35f1 1519}
1520//_________________________________________________________________________
3a0f6479 1521Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
170c35f1 1522#ifdef ALI_DATE
3a0f6479 1523 eventHeaderStruct *event,
1524 Bool_t nocheck
170c35f1 1525#else
3a0f6479 1526 eventHeaderStruct* /*event*/,
1527 Bool_t /*nocheck*/
170c35f1 1528
1529#endif
1530 )
1531{
1532 //
1533 // process date event
1534 //
1535#ifdef ALI_DATE
1536 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
3a0f6479 1537 Int_t result=ProcessEventDAQ(rawReader, nocheck);
170c35f1 1538 delete rawReader;
1539 return result;
1540#else
1541 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1542 return 0;
1543#endif
1544
1545}
1546//____________Online trackling in AliTRDtrigger________________________________
1547Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1548{
1549 //
1550 // For the DAQ
1551 // Fill a simple average pulse height
1552 //
1553
1554 // Localisation of the Xbins involved
1555 //LocalisationDetectorXbins(det);
1556
1557 // Row and position in the pad groups
1558 //Int_t posr = 0;
1559 //if (fCalibraMode->GetNnZ(1) != 0) {
1560 // posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1561 //}
1562
1563 // Col of cluster and position in the pad groups
1564 //Int_t posc = 0;
1565 //if (fCalibraMode->GetNnRphi(1) != 0) {
1566 // posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1567 //}
1568
1569 //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1570
3a0f6479 1571 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
170c35f1 1572
1573 return kTRUE;
1574
1575}
170c35f1 1576//____________Write_____________________________________________________
1577//_____________________________________________________________________
1578void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1579{
3a0f6479 1580 //
1581 // Write infos to file
1582 //
1583
170c35f1 1584 //For debugging
1585 if ( fDebugStreamer ) {
1586 delete fDebugStreamer;
1587 fDebugStreamer = 0x0;
55a288e5 1588 }
170c35f1 1589
55a288e5 1590 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
170c35f1 1591 ,fNumberTrack
1592 ,fNumberUsedCh[0]
1593 ,fNumberUsedCh[1]
1594 ,fNumberUsedPh[0]
1595 ,fNumberUsedPh[1]));
1596
1597 TDirectory *backup = gDirectory;
1598 TString option;
1599
1600 if ( append )
1601 option = "update";
1602 else
1603 option = "recreate";
1604
1605 TFile f(filename,option.Data());
55a288e5 1606
1607 TStopwatch stopwatch;
1608 stopwatch.Start();
3a0f6479 1609 if(fVector2d) {
1610 f.WriteTObject(fCalibraVector);
1611 }
55a288e5 1612
170c35f1 1613 if (fCH2dOn ) {
55a288e5 1614 if (fHisto2d) {
170c35f1 1615 f.WriteTObject(fCH2d);
55a288e5 1616 }
55a288e5 1617 }
170c35f1 1618 if (fPH2dOn ) {
55a288e5 1619 if (fHisto2d) {
170c35f1 1620 f.WriteTObject(fPH2d);
55a288e5 1621 }
55a288e5 1622 }
170c35f1 1623 if (fPRF2dOn) {
55a288e5 1624 if (fHisto2d) {
170c35f1 1625 f.WriteTObject(fPRF2d);
55a288e5 1626 }
55a288e5 1627 }
3a0f6479 1628 if(fLinearFitterOn){
1629 AnalyseLinearFitter();
1630 f.WriteTObject(fLinearVdriftFit);
170c35f1 1631 }
3a0f6479 1632
170c35f1 1633 f.Close();
55a288e5 1634
170c35f1 1635 if ( backup ) backup->cd();
55a288e5 1636
170c35f1 1637 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1638 ,stopwatch.RealTime(),stopwatch.CpuTime()));
55a288e5 1639}
170c35f1 1640//___________________________________________probe the histos__________________________________________________
55a288e5 1641Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1642{
1643 //
1644 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1645 // debug mode with 2 for TH2I and 3 for TProfile2D
1646 // It gives a pointer to a Double_t[7] with the info following...
1647 // [0] : number of calibration groups with entries
1648 // [1] : minimal number of entries found
1649 // [2] : calibration group number of the min
1650 // [3] : maximal number of entries found
1651 // [4] : calibration group number of the max
1652 // [5] : mean number of entries found
3a0f6479 1653 // [6] : mean relative error
55a288e5 1654 //
1655
1656 Double_t *info = new Double_t[7];
1657
1658 // Number of Xbins (detectors or groups of pads)
170c35f1 1659 Int_t nbins = h->GetNbinsY(); //number of calibration groups
1660 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
55a288e5 1661
1662 // Initialise
1663 Double_t nbwe = 0; //number of calibration groups with entries
1664 Double_t minentries = 0; //minimal number of entries found
1665 Double_t maxentries = 0; //maximal number of entries found
1666 Double_t placemin = 0; //calibration group number of the min
1667 Double_t placemax = -1; //calibration group number of the max
1668 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1669 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1670
1671 Double_t counter = 0;
1672
1673 //Debug
0c349049 1674 TH1F *nbEntries = 0x0;//distribution of the number of entries
1675 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
1676 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
55a288e5 1677
1678 // Beginning of the loop over the calibration groups
1679 for (Int_t idect = 0; idect < nbins; idect++) {
1680
170c35f1 1681 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
55a288e5 1682 projch->SetDirectory(0);
1683
1684 // Number of entries for this calibration group
1685 Double_t nentries = 0.0;
1686 if((i%2) == 0){
170c35f1 1687 for (Int_t k = 0; k < nxbins; k++) {
1688 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
55a288e5 1689 }
1690 }
1691 else{
170c35f1 1692 for (Int_t k = 0; k < nxbins; k++) {
1693 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1694 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1695 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
55a288e5 1696 counter++;
1697 }
1698 }
1699 }
1700
1701 //Debug
1702 if(i > 1){
0c349049 1703 if((!((Bool_t)nbEntries)) && (nentries > 0)){
1704 nbEntries = new TH1F("Number of entries","Number of entries"
55a288e5 1705 ,100,(Int_t)nentries/2,nentries*2);
0c349049 1706 nbEntries->SetDirectory(0);
1707 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
55a288e5 1708 ,nbins,0,nbins);
0c349049 1709 nbEntriesPerGroup->SetDirectory(0);
1710 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
55a288e5 1711 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
0c349049 1712 nbEntriesPerSp->SetDirectory(0);
55a288e5 1713 }
0c349049 1714 if(nbEntries){
1715 if(nentries > 0) nbEntries->Fill(nentries);
1716 nbEntriesPerGroup->Fill(idect+0.5,nentries);
1717 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
55a288e5 1718 }
1719 }
1720
1721 //min amd max
1722 if(nentries > maxentries){
1723 maxentries = nentries;
1724 placemax = idect;
1725 }
1726 if(idect == 0) {
1727 minentries = nentries;
1728 }
1729 if(nentries < minentries){
1730 minentries = nentries;
1731 placemin = idect;
1732 }
1733 //nbwe
1734 if(nentries > 0) {
1735 nbwe++;
1736 meanstats += nentries;
1737 }
55a288e5 1738 }//calibration groups loop
1739
1740 if(nbwe > 0) meanstats /= nbwe;
1741 if(counter > 0) meanrelativerror /= counter;
1742
1743 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1744 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1745 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1746 AliInfo(Form("The mean number of entries is %f",meanstats));
1747 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1748
1749 info[0] = nbwe;
1750 info[1] = minentries;
1751 info[2] = placemin;
1752 info[3] = maxentries;
1753 info[4] = placemax;
1754 info[5] = meanstats;
1755 info[6] = meanrelativerror;
1756
1757 if(i > 1){
1758 gStyle->SetPalette(1);
1759 gStyle->SetOptStat(1111);
1760 gStyle->SetPadBorderMode(0);
1761 gStyle->SetCanvasColor(10);
1762 gStyle->SetPadLeftMargin(0.13);
1763 gStyle->SetPadRightMargin(0.01);
1764 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1765 stat->Divide(2,1);
1766 stat->cd(1);
0c349049 1767 nbEntries->Draw("");
55a288e5 1768 stat->cd(2);
0c349049 1769 nbEntriesPerSp->SetStats(0);
1770 nbEntriesPerSp->Draw("");
55a288e5 1771 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1772 stat1->cd();
0c349049 1773 nbEntriesPerGroup->SetStats(0);
1774 nbEntriesPerGroup->Draw("");
55a288e5 1775 }
1776
1777 return info;
1778
1779}
170c35f1 1780//____________________________________________________________________________
1781Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1782{
1783 //
1784 // Return a Int_t[4] with:
1785 // 0 Mean number of entries
1786 // 1 median of number of entries
1787 // 2 rms of number of entries
1788 // 3 number of group with entries
1789 //
1790
1791 Double_t *stat = new Double_t[4];
1792 stat[3] = 0.0;
1793
1794 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
1795 Double_t *weight = new Double_t[nbofgroups];
1796 Int_t *nonul = new Int_t[nbofgroups];
1797
1798 for(Int_t k = 0; k < nbofgroups; k++){
1799 if(fEntriesCH[k] > 0) {
1800 weight[k] = 1.0;
1801 nonul[(Int_t)stat[3]] = fEntriesCH[k];
1802 stat[3]++;
1803 }
1804 else weight[k] = 0.0;
1805 }
1806 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
1807 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
1808 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1809
1810 return stat;
1811
1812}
1813//____________________________________________________________________________
1814Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1815{
1816 //
1817 // Return a Int_t[4] with:
1818 // 0 Mean number of entries
1819 // 1 median of number of entries
1820 // 2 rms of number of entries
1821 // 3 number of group with entries
1822 //
1823
1824 Double_t *stat = new Double_t[4];
1825 stat[3] = 0.0;
1826
1827 Int_t nbofgroups = 540;
1828 Double_t *weight = new Double_t[nbofgroups];
1829 Int_t *nonul = new Int_t[nbofgroups];
1830
1831 for(Int_t k = 0; k < nbofgroups; k++){
1832 if(fEntriesLinearFitter[k] > 0) {
1833 weight[k] = 1.0;
1834 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1835 stat[3]++;
1836 }
1837 else weight[k] = 0.0;
1838 }
1839 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1840 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1841 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1842
1843 return stat;
1844
1845}
1846//_____________________________________________________________________________
1847void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1848{
1849 //
1850 // Should be between 0 and 6
1851 //
1852
1853 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1854 AliInfo("The number of groups must be between 0 and 6!");
1855 }
1856 else {
1857 fNgroupprf = numberGroupsPRF;
1858 }
55a288e5 1859
170c35f1 1860}
55a288e5 1861//_____________________________________________________________________________
1862void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1863{
1864 //
1865 // Set the factor that will divide the deposited charge
1866 // to fit in the histo range [0,300]
1867 //
1868
1869 if (RelativeScale > 0.0) {
1870 fRelativeScale = RelativeScale;
1871 }
1872 else {
1873 AliInfo("RelativeScale must be strict positif!");
1874 }
1875
1876}
1877
1878//_____________________________________________________________________________
1879void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1880{
1881 //
1882 // Set the mode of calibration group in the z direction for the parameter i
1883 //
1884
1885 if ((Nz >= 0) &&
1886 (Nz < 5)) {
1887 fCalibraMode->SetNz(i, Nz);
1888 }
1889 else {
1890 AliInfo("You have to choose between 0 and 4");
1891 }
1892
1893}
1894
1895//_____________________________________________________________________________
1896void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1897{
1898 //
1899 // Set the mode of calibration group in the rphi direction for the parameter i
1900 //
1901
1902 if ((Nrphi >= 0) &&
1903 (Nrphi < 7)) {
1904 fCalibraMode->SetNrphi(i ,Nrphi);
1905 }
1906 else {
1907 AliInfo("You have to choose between 0 and 6");
1908 }
1909
1910}
55a288e5 1911//____________Protected Functions______________________________________________
1912//____________Create the 2D histo to be filled online__________________________
1913//
55a288e5 1914//_____________________________________________________________________________
1915void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1916{
1917 //
170c35f1 1918 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1919 // If fNgroupprf is zero then no binning in tnp
55a288e5 1920 //
1921
1922 TString name("Nz");
1923 name += fCalibraMode->GetNz(2);
1924 name += "Nrphi";
1925 name += fCalibraMode->GetNrphi(2);
170c35f1 1926 name += "Ngp";
1927 name += fNgroupprf;
55a288e5 1928
170c35f1 1929 if(fNgroupprf != 0){
1930
1931 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1932 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1933 fPRF2d->SetYTitle("Det/pad groups");
1934 fPRF2d->SetXTitle("Position x/W [pad width units]");
1935 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1936 fPRF2d->SetStats(0);
1937 }
1938 else{
1939 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1940 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1941 fPRF2d->SetYTitle("Det/pad groups");
1942 fPRF2d->SetXTitle("Position x/W [pad width units]");
1943 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1944 fPRF2d->SetStats(0);
1945 }
55a288e5 1946
1947}
1948
1949//_____________________________________________________________________________
1950void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1951{
1952 //
1953 // Create the 2D histos
1954 //
1955
1956 TString name("Nz");
1957 name += fCalibraMode->GetNz(1);
1958 name += "Nrphi";
1959 name += fCalibraMode->GetNrphi(1);
170c35f1 1960
55a288e5 1961 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
170c35f1 1962 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1963 ,nn,0,nn);
1964 fPH2d->SetYTitle("Det/pad groups");
1965 fPH2d->SetXTitle("time [#mus]");
55a288e5 1966 fPH2d->SetZTitle("<PH> [a.u.]");
1967 fPH2d->SetStats(0);
1968
1969}
55a288e5 1970//_____________________________________________________________________________
1971void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1972{
1973 //
1974 // Create the 2D histos
1975 //
1976
1977 TString name("Nz");
1978 name += fCalibraMode->GetNz(0);
1979 name += "Nrphi";
1980 name += fCalibraMode->GetNrphi(0);
1981
1982 fCH2d = new TH2I("CH2d",(const Char_t *) name
170c35f1 1983 ,fNumberBinCharge,0,300,nn,0,nn);
1984 fCH2d->SetYTitle("Det/pad groups");
1985 fCH2d->SetXTitle("charge deposit [a.u]");
55a288e5 1986 fCH2d->SetZTitle("counts");
1987 fCH2d->SetStats(0);
1988 fCH2d->Sumw2();
1989
1990}
3a0f6479 1991
55a288e5 1992//____________Offine tracking in the AliTRDtracker_____________________________
1993void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1994{
1995 //
1996 // For the offline tracking or mcm tracklets
1997 // This function will be called in the functions UpdateHistogram...
1998 // to fill the info of a track for the relativ gain calibration
1999 //
2000
170c35f1 2001 Int_t nb = 0; // Nombre de zones traversees
2002 Int_t fd = -1; // Premiere zone non nulle
2003 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
3a0f6479 2004
2005
55a288e5 2006
2007
2008 // See if the track goes through different zones
2009 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
2010 if (fAmpTotal[k] > 0.0) {
170c35f1 2011 totalcharge += fAmpTotal[k];
55a288e5 2012 nb++;
2013 if (nb == 1) {
2014 fd = k;
2015 }
2016 }
2017 }
55a288e5 2018
3a0f6479 2019
170c35f1 2020 switch (nb)
2021 {
2022 case 1:
55a288e5 2023 fNumberUsedCh[0]++;
170c35f1 2024 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2025 if (fHisto2d) {
170c35f1 2026 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2027 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2028 }
2029 if (fVector2d) {
3a0f6479 2030 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2031 }
170c35f1 2032 break;
2033 case 2:
55a288e5 2034 if ((fAmpTotal[fd] > 0.0) &&
170c35f1 2035 (fAmpTotal[fd+1] > 0.0)) {
55a288e5 2036 // One of the two very big
2037 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2038 if (fHisto2d) {
170c35f1 2039 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2040 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2041 }
2042 if (fVector2d) {
3a0f6479 2043 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2044 }
2045 fNumberUsedCh[1]++;
170c35f1 2046 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2047 }
2048 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2049 if (fHisto2d) {
170c35f1 2050 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
2051 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
55a288e5 2052 }
2053 if (fVector2d) {
3a0f6479 2054 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
55a288e5 2055 }
2056 fNumberUsedCh[1]++;
170c35f1 2057 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
55a288e5 2058 }
2059 }
55a288e5 2060 if (fCalibraMode->GetNfragZ(0) > 1) {
2061 if (fAmpTotal[fd] > 0.0) {
2062 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2063 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2064 // One of the two very big
2065 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2066 if (fHisto2d) {
170c35f1 2067 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2068 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2069 }
2070 if (fVector2d) {
3a0f6479 2071 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2072 }
2073 fNumberUsedCh[1]++;
170c35f1 2074 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2075 }
2076 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2077 if (fHisto2d) {
170c35f1 2078 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
2079 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
55a288e5 2080 }
2081 fNumberUsedCh[1]++;
170c35f1 2082 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
55a288e5 2083 if (fVector2d) {
3a0f6479 2084 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
55a288e5 2085 }
2086 }
2087 }
2088 }
2089 }
2090 }
170c35f1 2091 break;
2092 default: break;
55a288e5 2093 }
55a288e5 2094}
55a288e5 2095//____________Offine tracking in the AliTRDtracker_____________________________
2096void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2097{
2098 //
2099 // For the offline tracking or mcm tracklets
2100 // This function will be called in the functions UpdateHistogram...
2101 // to fill the info of a track for the drift velocity calibration
2102 //
2103
2104 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2105 Int_t fd1 = -1; // Premiere zone non nulle
2106 Int_t fd2 = -1; // Deuxieme zone non nulle
2107 Int_t k1 = -1; // Debut de la premiere zone
2108 Int_t k2 = -1; // Debut de la seconde zone
2109
2110 // See if the track goes through different zones
2111 for (Int_t k = 0; k < fTimeMax; k++) {
2112 if (fPHValue[k] > 0.0) {
2113 if (fd1 == -1) {
2114 fd1 = fPHPlace[k];
2115 k1 = k;
2116 }
2117 if (fPHPlace[k] != fd1) {
2118 if (fd2 == -1) {
2119 k2 = k;
2120 fd2 = fPHPlace[k];
2121 nb = 2;
2122 }
2123 if (fPHPlace[k] != fd2) {
2124 nb = 3;
2125 }
2126 }
2127 }
2128 }
170c35f1 2129
3a0f6479 2130
170c35f1 2131 switch(nb)
2132 {
2133 case 1:
2134 fNumberUsedPh[0]++;
2135 for (Int_t i = 0; i < fTimeMax; i++) {
2136 if (fHisto2d) {
2137 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2138 }
2139 if (fVector2d) {
3a0f6479 2140 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
55a288e5 2141 }
2142 }
170c35f1 2143 break;
2144 case 2:
2145 if ((fd1 == fd2+1) ||
2146 (fd2 == fd1+1)) {
2147 // One of the two fast all the think
2148 if (k2 > (k1+fDifference)) {
2149 //we choose to fill the fd1 with all the values
2150 fNumberUsedPh[1]++;
2151 for (Int_t i = 0; i < fTimeMax; i++) {
2152 if (fHisto2d) {
2153 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2154 }
2155 if (fVector2d) {
3a0f6479 2156 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2157 }
55a288e5 2158 }
170c35f1 2159 }
2160 if ((k2+fDifference) < fTimeMax) {
2161 //we choose to fill the fd2 with all the values
2162 fNumberUsedPh[1]++;
2163 for (Int_t i = 0; i < fTimeMax; i++) {
2164 if (fHisto2d) {
2165 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2166 }
55a288e5 2167 if (fVector2d) {
3a0f6479 2168 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
55a288e5 2169 }
170c35f1 2170 }
55a288e5 2171 }
2172 }
170c35f1 2173 // Two zones voisines sinon rien!
2174 if (fCalibraMode->GetNfragZ(1) > 1) {
2175 // Case 2
2176 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2177 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2178 // One of the two fast all the think
2179 if (k2 > (k1+fDifference)) {
2180 //we choose to fill the fd1 with all the values
2181 fNumberUsedPh[1]++;
2182 for (Int_t i = 0; i < fTimeMax; i++) {
2183 if (fHisto2d) {
2184 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2185 }
2186 if (fVector2d) {
3a0f6479 2187 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2188 }
55a288e5 2189 }
2190 }
170c35f1 2191 if ((k2+fDifference) < fTimeMax) {
2192 //we choose to fill the fd2 with all the values
2193 fNumberUsedPh[1]++;
2194 for (Int_t i = 0; i < fTimeMax; i++) {
2195 if (fHisto2d) {
2196 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2197 }
2198 if (fVector2d) {
3a0f6479 2199 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
170c35f1 2200 }
55a288e5 2201 }
2202 }
2203 }
2204 }
170c35f1 2205 // Two zones voisines sinon rien!
2206 // Case 3
2207 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2208 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2209 // One of the two fast all the think
2210 if (k2 > (k1 + fDifference)) {
2211 //we choose to fill the fd1 with all the values
2212 fNumberUsedPh[1]++;
2213 for (Int_t i = 0; i < fTimeMax; i++) {
2214 if (fHisto2d) {
2215 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2216 }
2217 if (fVector2d) {
3a0f6479 2218 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2219 }
55a288e5 2220 }
2221 }
170c35f1 2222 if ((k2+fDifference) < fTimeMax) {
2223 //we choose to fill the fd2 with all the values
2224 fNumberUsedPh[1]++;
2225 for (Int_t i = 0; i < fTimeMax; i++) {
2226 if (fHisto2d) {
2227 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2228 }
2229 if (fVector2d) {
3a0f6479 2230 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
170c35f1 2231 }
55a288e5 2232 }
2233 }
2234 }
2235 }
2236 }
170c35f1 2237 break;
2238 default: break;
2239 }
55a288e5 2240}
170c35f1 2241//____________Offine tracking in the AliTRDtracker_____________________________
2242Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
55a288e5 2243{
2244 //
170c35f1 2245 // For the offline tracking
2246 // This function will be called in the functions UpdateHistogram...
2247 // to fill the find the parameter P1 of a track for the drift velocity calibration
55a288e5 2248 //
3a0f6479 2249
2250
170c35f1 2251 //Number of points: if less than 3 return kFALSE
0c349049 2252 Int_t npoints = fListClusters->GetEntriesFast();
2253 if(npoints <= 2) return kFALSE;
170c35f1 2254
2255 //Variables
2256 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2257 Double_t snp = 0.0; // sin angle in the plan yx track
2258 Double_t y = 0.0; // y clusters in the middle of the chamber
2259 Double_t z = 0.0; // z cluster in the middle of the chamber
2260 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2261 Double_t tnp = 0.0; // tan angle in the plan xy track
2262 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2263 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2264 Double_t pointError = 0.0; // error after straight line fit
2265 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
2266 Int_t snpright = 1; // if we took in the middle snp
2267 Int_t crossrow = 0; // if it crosses a pad row
2268 Double_t tiltingangle = 0; // tiltingangle of the pad
2269 Float_t dzdx = 0; // dz/dx now from dz/dl
2270 Int_t nbli = 0; // number linear fitter points
f162af62 2271 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
170c35f1 2272
2273 linearFitterTracklet.StoreData(kFALSE);
2274 linearFitterTracklet.ClearPoints();
2275
2276 //if more than one row
2277 Int_t rowp = -1; // if it crosses a pad row
2278
2279 //tiltingangle
2280 tiltingangle = padplane->GetTiltingAngle();
2281 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2282
2283 //Fill with points
0c349049 2284 for(Int_t k = 0; k < npoints; k++){
170c35f1 2285
2286 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2287 Double_t ycluster = cl->GetY();
2288 Int_t time = cl->GetLocalTimeBin();
2289 Double_t timeis = time/fSf;
2290 //See if cross two pad rows
2291 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2292 if(k==0) rowp = row;
2293 if(row != rowp) crossrow = 1;
2294 //Take in the middle of the chamber
3a0f6479 2295 //FollowBack
170c35f1 2296 if(time > (Int_t) 10) {
3a0f6479 2297 //Follow
2298 //if(time < (Int_t) 11) {
170c35f1 2299 z = cl->GetZ();
2300 y = cl->GetY();
2301 snp = fPar2[time];
2302 tgl = fPar3[time];
2303 }
2304 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2305 nbli++;
55a288e5 2306 }
3a0f6479 2307 //FollowBack
2308 if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
2309 //Follow
2310 //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
170c35f1 2311 if(nbli <= 2) return kFALSE;
2312
2313 // Do the straight line fit now
2314 TVectorD pars;
2315 linearFitterTracklet.Eval();
2316 linearFitterTracklet.GetParameters(pars);
2317 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2318 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2319 dydt = pars[1];
2320
2321 if( TMath::Abs(snp) < 1.){
2322 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2323 }
2324 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2325
2326 if(fDebugLevel > 0){
2327 if ( !fDebugStreamer ) {
2328 //debug stream
2329 TDirectory *backup = gDirectory;
3a0f6479 2330 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
170c35f1 2331 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2332 }
2333
2334 (* fDebugStreamer) << "VDRIFT0"<<
0c349049 2335 "npoints="<<npoints<<
170c35f1 2336 "\n";
2337
2338
2339 (* fDebugStreamer) << "VDRIFT"<<
2340 "snpright="<<snpright<<
0c349049 2341 "npoints="<<npoints<<
170c35f1 2342 "nbli="<<nbli<<
2343 "detector="<<detector<<
2344 "snp="<<snp<<
2345 "tnp="<<tnp<<
2346 "tgl="<<tgl<<
2347 "tnt="<<tnt<<
2348 "y="<<y<<
2349 "z="<<z<<
2350 "dydt="<<dydt<<
2351 "dzdx="<<dzdx<<
2352 "crossrow="<<crossrow<<
2353 "errorpar="<<errorpar<<
2354 "pointError="<<pointError<<
2355 "\n";
55a288e5 2356
55a288e5 2357 }
2358
0c349049 2359 if(npoints < fNumberClusters) return kFALSE;
170c35f1 2360 if(snpright == 0) return kFALSE;
2361 if(pointError >= 0.1) return kFALSE;
2362 if(crossrow == 1) return kFALSE;
2363
2364 if(fLinearFitterOn){
2365 //Add to the linear fitter of the detector
2366 if( TMath::Abs(snp) < 1.){
2367 Double_t x = tnp-dzdx*tnt;
2368 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2369 if(fLinearFitterDebugOn) {
3a0f6479 2370 fLinearVdriftFit->Update(detector,x,pars[1]);
170c35f1 2371 }
2372 fEntriesLinearFitter[detector]++;
2373 }
2374 }
2375 //AliInfo("End of FindP1TrackPH with success!")
55a288e5 2376 return kTRUE;
2377
2378}
170c35f1 2379//____________Offine tracking in the AliTRDtracker_____________________________
2380Bool_t AliTRDCalibraFillHisto::HandlePRF()
2381{
2382 //
2383 // For the offline tracking
2384 // Fit the tracklet with a line and take the position as reference for the PRF
2385 //
3a0f6479 2386
170c35f1 2387 //Number of points
0c349049 2388 Int_t npoints = fListClusters->GetEntriesFast(); // number of total points
2389 Int_t nb3pc = 0; // number of three pads clusters used for fit
170c35f1 2390 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
2391
2392
2393 // To see the difference due to the fit
2394 Double_t *padPositions;
0c349049 2395 padPositions = new Double_t[npoints];
2396 for(Int_t k = 0; k < npoints; k++){
170c35f1 2397 padPositions[k] = 0.0;
2398 }
2399
2400
2401 //Find the position by a fit
2402 TLinearFitter fitter(2,"pol1");
2403 fitter.StoreData(kFALSE);
2404 fitter.ClearPoints();
0c349049 2405 for(Int_t k = 0; k < npoints; k++){
170c35f1 2406 //Take the cluster
2407 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2408 Short_t *signals = cl->GetSignals();
2409 Double_t time = cl->GetLocalTimeBin();
2410 //Calculate x if possible
2411 Float_t xcenter = 0.0;
2412 Bool_t echec = kTRUE;
2413 if((time<=7) || (time>=21)) continue;
2414 // Center 3 balanced: position with the center of the pad
2415 if ((((Float_t) signals[3]) > 0.0) &&
2416 (((Float_t) signals[2]) > 0.0) &&
2417 (((Float_t) signals[4]) > 0.0)) {
2418 echec = kFALSE;
2419 // Security if the denomiateur is 0
2420 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2421 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2422 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2423 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2424 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2425 }
2426 else {
2427 echec = kTRUE;
2428 }
2429 }
2430 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2431 if(echec) continue;
2432 //if no echec: calculate with the position of the pad
2433 // Position of the cluster
af26ce80 2434 Double_t padPosition = xcenter + cl->GetPadCol();
170c35f1 2435 padPositions[k] = padPosition;
0c349049 2436 nb3pc++;
170c35f1 2437 fitter.AddPoint(&time, padPosition,1);
2438 }//clusters loop
2439
0c349049 2440 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2441 if(nb3pc < 3) return kFALSE;
170c35f1 2442 fitter.Eval();
2443 TVectorD line(2);
2444 fitter.GetParameters(line);
2445 Float_t pointError = -1.0;
0c349049 2446 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
170c35f1 2447
2448
2449 // Now fill the PRF
0c349049 2450 for(Int_t k = 0; k < npoints; k++){
170c35f1 2451 //Take the cluster
2452 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2453 Short_t *signals = cl->GetSignals(); // signal
2454 Double_t time = cl->GetLocalTimeBin(); // time bin
2455 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
af26ce80 2456 Float_t padPos = cl->GetPadCol(); // middle pad
170c35f1 2457 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2458 Float_t ycenter = 0.0; // relative center charge
2459 Float_t ymin = 0.0; // relative left charge
2460 Float_t ymax = 0.0; // relative right charge
2461 Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx
2462 Double_t pt = fPar4[(Int_t)time]; // pt
2463 Float_t dzdx = 0.0; // dzdx
2464
2465
2466 //Requiere simply two pads clusters at least
2467 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2468 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2469 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2470 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2471 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2472 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2473 }
2474
2475 //calibration group
2476 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2477 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
3a0f6479 2478 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
170c35f1 2479 Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan
2480 Float_t xcl = cl->GetY(); // y cluster
2481 Float_t qcl = cl->GetQ(); // charge cluster
2482 Int_t plane = GetPlane(detector); // plane
2483 Int_t chamber = GetChamber(detector); // chamber
2484 Double_t xdiff = dpad; // reconstructed position constant
2485 Double_t x = dpad; // reconstructed position moved
0c349049 2486 Float_t ep = pointError; // error of fit
170c35f1 2487 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2488 Float_t signal3 = (Float_t)signals[3]; // signal
2489 Float_t signal2 = (Float_t)signals[2]; // signal
2490 Float_t signal4 = (Float_t)signals[4]; // signal
2491 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2492 Float_t tnp = 0.0;
2493 if(TMath::Abs(snp) < 1.0){
2494 tnp = snp / (TMath::Sqrt(1-snp*snp));
2495 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2496 }
2497
55a288e5 2498
170c35f1 2499 if(fDebugLevel > 0){
2500 if ( !fDebugStreamer ) {
2501 //debug stream
2502 TDirectory *backup = gDirectory;
3a0f6479 2503 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
170c35f1 2504 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2505 }
2506
2507 (* fDebugStreamer) << "PRF0"<<
2508 "caligroup="<<caligroup<<
2509 "detector="<<detector<<
2510 "plane="<<plane<<
2511 "chamber="<<chamber<<
0c349049 2512 "npoints="<<npoints<<
2513 "Np="<<nb3pc<<
2514 "ep="<<ep<<
170c35f1 2515 "snp="<<snp<<
2516 "tnp="<<tnp<<
2517 "tgl="<<tgl<<
2518 "dzdx="<<dzdx<<
2519 "pt="<<pt<<
2520 "padPos="<<padPos<<
2521 "padPosition="<<padPositions[k]<<
2522 "padPosTracklet="<<padPosTracklet<<
2523 "x="<<x<<
2524 "ycenter="<<ycenter<<
2525 "ymin="<<ymin<<
2526 "ymax="<<ymax<<
2527 "xcl="<<xcl<<
2528 "qcl="<<qcl<<
2529 "signal1="<<signal1<<
2530 "signal2="<<signal2<<
2531 "signal3="<<signal3<<
2532 "signal4="<<signal4<<
2533 "signal5="<<signal5<<
2534 "time="<<time<<
2535 "\n";
2536 x = xdiff;
2537 Int_t type=0;
2538 Float_t y = ycenter;
2539 (* fDebugStreamer) << "PRFALL"<<
2540 "caligroup="<<caligroup<<
2541 "detector="<<detector<<
2542 "plane="<<plane<<
2543 "chamber="<<chamber<<
0c349049 2544 "npoints="<<npoints<<
2545 "Np="<<nb3pc<<
2546 "ep="<<ep<<
170c35f1 2547 "type="<<type<<
2548 "snp="<<snp<<
2549 "tnp="<<tnp<<
2550 "tgl="<<tgl<<
2551 "dzdx="<<dzdx<<
2552 "pt="<<pt<<
2553 "padPos="<<padPos<<
2554 "padPosition="<<padPositions[k]<<
2555 "padPosTracklet="<<padPosTracklet<<
2556 "x="<<x<<
2557 "y="<<y<<
2558 "xcl="<<xcl<<
2559 "qcl="<<qcl<<
2560 "signal1="<<signal1<<
2561 "signal2="<<signal2<<
2562 "signal3="<<signal3<<
2563 "signal4="<<signal4<<
2564 "signal5="<<signal5<<
2565 "time="<<time<<
2566 "\n";
2567 x=-(xdiff+1);
2568 y = ymin;
2569 type=-1;
2570 (* fDebugStreamer) << "PRFALL"<<
2571 "caligroup="<<caligroup<<
2572 "detector="<<detector<<
2573 "plane="<<plane<<
2574 "chamber="<<chamber<<
0c349049 2575 "npoints="<<npoints<<
2576 "Np="<<nb3pc<<
2577 "ep="<<ep<<
170c35f1 2578 "type="<<type<<
2579 "snp="<<snp<<
2580 "tnp="<<tnp<<
2581 "tgl="<<tgl<<
2582 "dzdx="<<dzdx<<
2583 "pt="<<pt<<
2584 "padPos="<<padPos<<
2585 "padPosition="<<padPositions[k]<<
2586 "padPosTracklet="<<padPosTracklet<<
2587 "x="<<x<<
2588 "y="<<y<<
2589 "xcl="<<xcl<<
2590 "qcl="<<qcl<<
2591 "signal1="<<signal1<<
2592 "signal2="<<signal2<<
2593 "signal3="<<signal3<<
2594 "signal4="<<signal4<<
2595 "signal5="<<signal5<<
2596 "time="<<time<<
2597 "\n";
2598 x=1-xdiff;
2599 y = ymax;
2600 type=1;
2601
2602 (* fDebugStreamer) << "PRFALL"<<
2603 "caligroup="<<caligroup<<
2604 "detector="<<detector<<
2605 "plane="<<plane<<
2606 "chamber="<<chamber<<
0c349049 2607 "npoints="<<npoints<<
2608 "Np="<<nb3pc<<
2609 "ep="<<ep<<
170c35f1 2610 "type="<<type<<
2611 "snp="<<snp<<
2612 "tnp="<<tnp<<
2613 "tgl="<<tgl<<
2614 "dzdx="<<dzdx<<
2615 "pt="<<pt<<
2616 "padPos="<<padPos<<
2617 "padPosition="<<padPositions[k]<<
2618 "padPosTracklet="<<padPosTracklet<<
2619 "x="<<x<<
2620 "y="<<y<<
2621 "xcl="<<xcl<<
2622 "qcl="<<qcl<<
2623 "signal1="<<signal1<<
2624 "signal2="<<signal2<<
2625 "signal3="<<signal3<<
2626 "signal4="<<signal4<<
2627 "signal5="<<signal5<<
2628 "time="<<time<<
2629 "\n";
2630
2631 }
2632
2633 // some cuts
0c349049 2634 if(npoints < fNumberClusters) continue;
2635 if(nb3pc <= 5) continue;
170c35f1 2636 if((time >= 21) || (time < 7)) continue;
2637 if(TMath::Abs(snp) >= 1.0) continue;
2638 if(qcl < 80) continue;
2639
2640 Bool_t echec = kFALSE;
2641 Double_t shift = 0.0;
2642 //Calculate the shift in x coresponding to this tnp
2643 if(fNgroupprf != 0.0){
2644 shift = -3.0*(fNgroupprf-1)-1.5;
2645 Double_t limithigh = -0.2*(fNgroupprf-1);
2646 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2647 else{
2648 while(tnp > limithigh){
2649 limithigh += 0.2;
2650 shift += 3.0;
2651 }
2652 }
2653 }
2654 if (fHisto2d && !echec) {
3a0f6479 2655 if(TMath::Abs(dpad) < 1.5) {
2656 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2657 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2658 }
2659 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
170c35f1 2660 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2661 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2662 }
3a0f6479 2663 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
170c35f1 2664 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2665 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2666 }
2667 }
2668 //Not equivalent anymore here!
2669 if (fVector2d && !echec) {
3a0f6479 2670 if(TMath::Abs(dpad) < 1.5) {
2671 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2672 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
170c35f1 2673 }
3a0f6479 2674 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2675 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2676 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2677 }
2678 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2679 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2680 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
170c35f1 2681 }
2682 }
2683 }
2684 return kTRUE;
2685
2686}
d0569428 2687//____________Offine tracking in the AliTRDtracker_____________________________
2688Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
55a288e5 2689{
2690 //
d0569428 2691 // For the offline tracking
2692 // This function will be called in the functions UpdateHistogram...
2693 // to fill the find the parameter P1 of a track for the drift velocity calibration
55a288e5 2694 //
2695
d0569428 2696
2697 //Number of points: if less than 3 return kFALSE
0c349049 2698 Int_t npoints = index1-index0;
2699 if(npoints <= 2) return kFALSE;
55a288e5 2700
d0569428 2701 //Variables
2702 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2703 Double_t snp = 0.0; // sin angle in the plan yx track
2704 Double_t y = 0.0; // y clusters in the middle of the chamber
2705 Double_t z = 0.0; // z cluster in the middle of the chamber
2706 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2707 Double_t tnp = 0.0; // tan angle in the plan xy track
2708 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2709 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2710 Double_t pointError = 0.0; // error after straight line fit
2711 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
2712 //Int_t snpright = 1; // if we took in the middle snp
2713 Int_t crossrow = 0; // if it crosses a pad row
2714 Double_t tiltingangle = 0; // tiltingangle of the pad
2715 Float_t dzdx = 0; // dz/dx now from dz/dl
2716 Int_t nbli = 0; // number linear fitter points
2717 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
55a288e5 2718
d0569428 2719 linearFitterTracklet.StoreData(kFALSE);
2720 linearFitterTracklet.ClearPoints();
2721
2722 //if more than one row
2723 Int_t rowp = -1; // if it crosses a pad row
55a288e5 2724
d0569428 2725 //tiltingangle
2726 tiltingangle = padplane->GetTiltingAngle();
2727 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
55a288e5 2728
d0569428 2729 //Fill with points
0c349049 2730 for(Int_t k = 0; k < npoints; k++){
d0569428 2731
2732 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2733 Double_t ycluster = cl->GetY();
2734 Int_t time = cl->GetLocalTimeBin();
2735 Double_t timeis = time/fSf;
2736 //See if cross two pad rows
2737 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2738 if(k==0) rowp = row;
2739 if(row != rowp) crossrow = 1;
2740 //Take in the middle of the chamber
2741 //FollowBack
2742 //if(time > (Int_t) 10) {
2743 //Follow
2744 if(time < (Int_t) 11) {
2745 z = cl->GetZ();
2746 y = cl->GetY();
2747 }
2748 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2749 nbli++;
2750 }
55a288e5 2751
d0569428 2752 // take now the snp, tnp and tgl from the track
2753 snp = t->GetSnpPlane(GetPlane(detector));
2754 tgl = t->GetTglPlane(GetPlane(detector));
170c35f1 2755
d0569428 2756 //FollowBack
2757 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0;
2758 //Follow
2759 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0;
2760 if(nbli <= 2) return kFALSE;
170c35f1 2761
d0569428 2762 // Do the straight line fit now
2763 TVectorD pars;
2764 linearFitterTracklet.Eval();
2765 linearFitterTracklet.GetParameters(pars);
2766 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2767 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2768 dydt = pars[1];
2769
2770 if( TMath::Abs(snp) < 1.){
2771 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2772 }
2773 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2774
2775 if(fDebugLevel > 0){
2776 if ( !fDebugStreamer ) {
2777 //debug stream
2778 TDirectory *backup = gDirectory;
2779 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2780 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2781 }
2782
2783
2784 (* fDebugStreamer) << "VDRIFT"<<
2785 //"snpright="<<snpright<<
0c349049 2786 "npoints="<<npoints<<
d0569428 2787 "nbli="<<nbli<<
2788 "detector="<<detector<<
2789 "snp="<<snp<<
2790 "tnp="<<tnp<<
2791 "tgl="<<tgl<<
2792 "tnt="<<tnt<<
2793 "y="<<y<<
2794 "z="<<z<<
2795 "dydt="<<dydt<<
2796 "dzdx="<<dzdx<<
2797 "crossrow="<<crossrow<<
2798 "errorpar="<<errorpar<<
2799 "pointError="<<pointError<<
2800 "\n";
2801
2802 }
2803
0c349049 2804 if(npoints < fNumberClusters) return kFALSE;
d0569428 2805 //if(snpright == 0) return kFALSE;
2806 if(pointError >= 0.1) return kFALSE;
2807 if(crossrow == 1) return kFALSE;
2808
2809 if(fLinearFitterOn){
2810 //Add to the linear fitter of the detector
2811 if( TMath::Abs(snp) < 1.){
2812 Double_t x = tnp-dzdx*tnt;
2813 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2814 if(fLinearFitterDebugOn) {
2815 fLinearVdriftFit->Update(detector,x,pars[1]);
2816 }
2817 fEntriesLinearFitter[detector]++;
2818 }
2819 }
2820 //AliInfo("End of FindP1TrackPH with success!")
2821 return kTRUE;
2822
2823}
2824//____________Offine tracking in the AliTRDtracker_____________________________
2825Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2826{
2827 //
2828 // For the offline tracking
2829 // Fit the tracklet with a line and take the position as reference for the PRF
2830 //
2831
2832 //Number of points
0c349049 2833 Int_t npoints = index1-index0; // number of total points
2834 Int_t nb3pc = 0; // number of three pads clusters used for fit
d0569428 2835 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
2836
2837
2838 // To see the difference due to the fit
2839 Double_t *padPositions;
0c349049 2840 padPositions = new Double_t[npoints];
2841 for(Int_t k = 0; k < npoints; k++){
d0569428 2842 padPositions[k] = 0.0;
2843 }
2844
2845
2846 //Find the position by a fit
2847 TLinearFitter fitter(2,"pol1");
2848 fitter.StoreData(kFALSE);
2849 fitter.ClearPoints();
0c349049 2850 for(Int_t k = 0; k < npoints; k++){
d0569428 2851 //Take the cluster
2852 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2853 Short_t *signals = cl->GetSignals();
2854 Double_t time = cl->GetLocalTimeBin();
2855 //Calculate x if possible
2856 Float_t xcenter = 0.0;
2857 Bool_t echec = kTRUE;
2858 if((time<=7) || (time>=21)) continue;
2859 // Center 3 balanced: position with the center of the pad
2860 if ((((Float_t) signals[3]) > 0.0) &&
2861 (((Float_t) signals[2]) > 0.0) &&
2862 (((Float_t) signals[4]) > 0.0)) {
2863 echec = kFALSE;
2864 // Security if the denomiateur is 0
2865 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2866 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2867 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2868 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2869 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2870 }
2871 else {
2872 echec = kTRUE;
2873 }
2874 }
2875 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2876 if(echec) continue;
2877 //if no echec: calculate with the position of the pad
2878 // Position of the cluster
af26ce80 2879 Double_t padPosition = xcenter + cl->GetPadCol();
d0569428 2880 padPositions[k] = padPosition;
0c349049 2881 nb3pc++;
d0569428 2882 fitter.AddPoint(&time, padPosition,1);
2883 }//clusters loop
2884
0c349049 2885 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2886 if(nb3pc < 3) return kFALSE;
d0569428 2887 fitter.Eval();
2888 TVectorD line(2);
2889 fitter.GetParameters(line);
2890 Float_t pointError = -1.0;
0c349049 2891 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
d0569428 2892
2893 // Take the tgl and snp with the track t now
2894 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
2895 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
2896 Float_t dzdx = 0.0; // dzdx
2897 Float_t tnp = 0.0;
2898 if(TMath::Abs(snp) < 1.0){
2899 tnp = snp / (TMath::Sqrt(1-snp*snp));
2900 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2901 }
2902
2903
2904 // Now fill the PRF
0c349049 2905 for(Int_t k = 0; k < npoints; k++){
d0569428 2906 //Take the cluster
2907 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2908 Short_t *signals = cl->GetSignals(); // signal
2909 Double_t time = cl->GetLocalTimeBin(); // time bin
2910 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
af26ce80 2911 Float_t padPos = cl->GetPadCol(); // middle pad
d0569428 2912 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2913 Float_t ycenter = 0.0; // relative center charge
2914 Float_t ymin = 0.0; // relative left charge
2915 Float_t ymax = 0.0; // relative right charge
2916
2917
2918
2919 //Requiere simply two pads clusters at least
2920 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2921 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2922 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2923 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2924 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2925 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2926 }
2927
2928 //calibration group
2929 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2930 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2931 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
2932 Float_t xcl = cl->GetY(); // y cluster
2933 Float_t qcl = cl->GetQ(); // charge cluster
2934 Int_t plane = GetPlane(detector); // plane
2935 Int_t chamber = GetChamber(detector); // chamber
2936 Double_t xdiff = dpad; // reconstructed position constant
2937 Double_t x = dpad; // reconstructed position moved
0c349049 2938 Float_t ep = pointError; // error of fit
d0569428 2939 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2940 Float_t signal3 = (Float_t)signals[3]; // signal
2941 Float_t signal2 = (Float_t)signals[2]; // signal
2942 Float_t signal4 = (Float_t)signals[4]; // signal
2943 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2944
2945
2946
2947 if(fDebugLevel > 0){
2948 if ( !fDebugStreamer ) {
2949 //debug stream
2950 TDirectory *backup = gDirectory;
2951 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2952 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2953 }
2954
2955 (* fDebugStreamer) << "PRF0"<<
2956 "caligroup="<<caligroup<<
2957 "detector="<<detector<<
2958 "plane="<<plane<<
2959 "chamber="<<chamber<<
0c349049 2960 "npoints="<<npoints<<
2961 "Np="<<nb3pc<<
2962 "ep="<<ep<<
d0569428 2963 "snp="<<snp<<
2964 "tnp="<<tnp<<
2965 "tgl="<<tgl<<
2966 "dzdx="<<dzdx<<
2967 "padPos="<<padPos<<
2968 "padPosition="<<padPositions[k]<<
2969 "padPosTracklet="<<padPosTracklet<<
2970 "x="<<x<<
2971 "ycenter="<<ycenter<<
2972 "ymin="<<ymin<<
2973 "ymax="<<ymax<<
2974 "xcl="<<xcl<<
2975 "qcl="<<qcl<<
2976 "signal1="<<signal1<<
2977 "signal2="<<signal2<<
2978 "signal3="<<signal3<<
2979 "signal4="<<signal4<<
2980 "signal5="<<signal5<<
2981 "time="<<time<<
2982 "\n";
2983 x = xdiff;
2984 Int_t type=0;
2985 Float_t y = ycenter;
2986 (* fDebugStreamer) << "PRFALL"<<
2987 "caligroup="<<caligroup<<
2988 "detector="<<detector<<
2989 "plane="<<plane<<
2990 "chamber="<<chamber<<
0c349049 2991 "npoints="<<npoints<<
2992 "Np="<<nb3pc<<
2993 "ep="<<ep<<
d0569428 2994 "type="<<type<<
2995 "snp="<<snp<<
2996 "tnp="<<tnp<<
2997 "tgl="<<tgl<<
2998 "dzdx="<<dzdx<<
2999 "padPos="<<padPos<<
3000 "padPosition="<<padPositions[k]<<
3001 "padPosTracklet="<<padPosTracklet<<
3002 "x="<<x<<
3003 "y="<<y<<
3004 "xcl="<<xcl<<
3005 "qcl="<<qcl<<
3006 "signal1="<<signal1<<
3007 "signal2="<<signal2<<
3008 "signal3="<<signal3<<
3009 "signal4="<<signal4<<
3010 "signal5="<<signal5<<
3011 "time="<<time<<
3012 "\n";
3013 x=-(xdiff+1);
3014 y = ymin;
3015 type=-1;
3016 (* fDebugStreamer) << "PRFALL"<<
3017 "caligroup="<<caligroup<<
3018 "detector="<<detector<<
3019 "plane="<<plane<<
3020 "chamber="<<chamber<<
0c349049 3021 "npoints="<<npoints<<
3022 "Np="<<nb3pc<<
3023 "ep="<<ep<<
d0569428 3024 "type="<<type<<
3025 "snp="<<snp<<
3026 "tnp="<<tnp<<
3027 "tgl="<<tgl<<
3028 "dzdx="<<dzdx<<
3029 "padPos="<<padPos<<
3030 "padPosition="<<padPositions[k]<<
3031 "padPosTracklet="<<padPosTracklet<<
3032 "x="<<x<<
3033 "y="<<y<<
3034 "xcl="<<xcl<<
3035 "qcl="<<qcl<<
3036 "signal1="<<signal1<<
3037 "signal2="<<signal2<<
3038 "signal3="<<signal3<<
3039 "signal4="<<signal4<<
3040 "signal5="<<signal5<<
3041 "time="<<time<<
3042 "\n";
3043 x=1-xdiff;
3044 y = ymax;
3045 type=1;
3046
3047 (* fDebugStreamer) << "PRFALL"<<
3048 "caligroup="<<caligroup<<
3049 "detector="<<detector<<
3050 "plane="<<plane<<
3051 "chamber="<<chamber<<
0c349049 3052 "npoints="<<npoints<<
3053 "Np="<<nb3pc<<
3054 "ep="<<ep<<
d0569428 3055 "type="<<type<<
3056 "snp="<<snp<<
3057 "tnp="<<tnp<<
3058 "tgl="<<tgl<<
3059 "dzdx="<<dzdx<<
3060 "padPos="<<padPos<<
3061 "padPosition="<<padPositions[k]<<
3062 "padPosTracklet="<<padPosTracklet<<
3063 "x="<<x<<
3064 "y="<<y<<
3065 "xcl="<<xcl<<
3066 "qcl="<<qcl<<
3067 "signal1="<<signal1<<
3068 "signal2="<<signal2<<
3069 "signal3="<<signal3<<
3070 "signal4="<<signal4<<
3071 "signal5="<<signal5<<
3072 "time="<<time<<
3073 "\n";
3074
3075 }
3076
3077 // some cuts
0c349049 3078 if(npoints < fNumberClusters) continue;
3079 if(nb3pc <= 5) continue;
d0569428 3080 if((time >= 21) || (time < 7)) continue;
3081 if(TMath::Abs(snp) >= 1.0) continue;
3082 if(qcl < 80) continue;
3083
3084 Bool_t echec = kFALSE;
3085 Double_t shift = 0.0;
3086 //Calculate the shift in x coresponding to this tnp
3087 if(fNgroupprf != 0.0){
3088 shift = -3.0*(fNgroupprf-1)-1.5;
3089 Double_t limithigh = -0.2*(fNgroupprf-1);
3090 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
3091 else{
3092 while(tnp > limithigh){
3093 limithigh += 0.2;
3094 shift += 3.0;
3095 }
3096 }
3097 }
3098 if (fHisto2d && !echec) {
3099 if(TMath::Abs(dpad) < 1.5) {
3100 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
3101 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
3102 }
3103 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3104 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
3105 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
3106 }
3107 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3108 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
3109 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
3110 }
3111 }
3112 //Not equivalent anymore here!
3113 if (fVector2d && !echec) {
3114 if(TMath::Abs(dpad) < 1.5) {
3115 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
3116 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
3117 }
3118 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3119 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
3120 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
3121 }
3122 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3123 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
3124 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
3125 }
3126 }
3127 }
3128 return kTRUE;
3129
3130}
3131//
3132//____________Some basic geometry function_____________________________________
3133//
3134//_____________________________________________________________________________
3135Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3136{
3137 //
3138 // Reconstruct the plane number from the detector number
3139 //
3140
3141 return ((Int_t) (d % 6));
3142
3143}
3144
3145//_____________________________________________________________________________
3146Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3147{
3148 //
3149 // Reconstruct the chamber number from the detector number
3150 //
3151 Int_t fgkNplan = 6;
3152
3153 return ((Int_t) (d % 30) / fgkNplan);
3154
3155}
3156
3157//_____________________________________________________________________________
3158Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3159{
3160 //
3161 // Reconstruct the sector number from the detector number
3162 //
3163 Int_t fg = 30;
3164
3165 return ((Int_t) (d / fg));
3166
3167}
3168//_____________________________________________________________________
3169TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3170{
3171 //
3172 // return pointer to fPH2d TProfile2D
3173 // create a new TProfile2D if it doesn't exist allready
3174 //
3175 if ( fPH2d )
3176 return fPH2d;
3177
3178 // Some parameters
3179 fTimeMax = nbtimebin;
3180 fSf = samplefrequency;
3181
3182 /*
3183 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
3184 ,fCalibraMode->GetNz(1)
3185 ,fCalibraMode->GetNrphi(1)));
170c35f1 3186
3187 // Calcul the number of Xbins
0c349049 3188 Int_t ntotal1 = 0;
170c35f1 3189 fCalibraMode->ModePadCalibration(2,1);
3190 fCalibraMode->ModePadFragmentation(0,2,0,1);
3191 fCalibraMode->SetDetChamb2(1);
0c349049 3192 ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
170c35f1 3193 fCalibraMode->ModePadCalibration(0,1);
3194 fCalibraMode->ModePadFragmentation(0,0,0,1);
3195 fCalibraMode->SetDetChamb0(1);
0c349049 3196 ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
3197 AliInfo(Form("Total number of Xbins: %d",ntotal1));
170c35f1 3198
0c349049 3199 CreatePH2d(ntotal1);
170c35f1 3200 */
3201
3202 CreatePH2d(540);
3203
3204 return fPH2d;
3205}
3206//_____________________________________________________________________
3207TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3208{
3209 //
3210 // return pointer to TLinearFitter Calibration
3211 // if force is true create a new TLinearFitter if it doesn't exist allready
3212 //
170c35f1 3213
d0569428 3214 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3215 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3216 }
3217
3218 // if we are forced and TLinearFitter doesn't yet exist create it
170c35f1 3219
d0569428 3220 // new TLinearFitter
3221 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3222 fLinearFitterArray.AddAt(linearfitter,detector);
3223 return linearfitter;
170c35f1 3224}
170c35f1 3225
170c35f1 3226//_____________________________________________________________________
3227void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3228{
3229 //
3230 // FillCH2d: Marian style
3231 //
3232
3233 //skip simply the value out of range
3234 if((y>=300.0) || (y<0.0)) return;
3235
3236 //Calcul the y place
3237 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3238 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3239
3240 //Fill
3241 fCH2d->GetArray()[place]++;
3242
3243}
3a0f6479 3244
3245//____________________________________________________________________________
3246void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3247{
3248 //
3249 // Analyse array of linear fitter because can not be written
3250 // Store two arrays: one with the param the other one with the error param + number of entries
3251 //
3252
3253 for(Int_t k = 0; k < 540; k++){
3254 TLinearFitter *linearfitter = GetLinearFitter(k);
3255 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3256 TVectorD *par = new TVectorD(2);
3257 TVectorD pare = TVectorD(2);
3258 TVectorD *parE = new TVectorD(3);
3259 linearfitter->Eval();
3260 linearfitter->GetParameters(*par);
3261 linearfitter->GetErrors(pare);
3262 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3263 (*parE)[0] = pare[0]*ppointError;
3264 (*parE)[1] = pare[1]*ppointError;
3265 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3266 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3267 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3268 }
3269 }
3270}
d0569428 3271//________________________________________________________________________________
3272Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const
3273{
3274 // Partie entiere of the (x+0.5)
3275
3276 int i;
3277 if (x >= (-0.5)) {
3278 i = int(x + 0.5);
3279 //if (x + 0.5 == Float_t(i) && i & 1) i--;
3280 } else {
3281 i = int(x - 0.5);
3282 //if (x - 0.5 == Float_t(i) && i & 1) i++;
3283 if((x-0.5)==i) i++;
3284
3285 }
3286 return i;
3287}