]> 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
Adding macros to check the raw data (Marian)
[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
1246 for(Int_t j = 0; j < 21; j++){
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
1280 for(Int_t k = 0; k < 36; k++){
1281 for(Int_t j = 0; j < 21; j++){
1282 phvalue[j][k] = baseline;
3a0f6479 1283 }
1284 }
1285 }
e4db522f 1286
3a0f6479 1287 fDetectorPreviousTrack = idetector;
e4db522f 1288 fMCMPrevious = imcm;
1289 fROBPrevious = irob;
1290
3a0f6479 1291 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1292 if(nbtimebin == 0) return 0;
1293 if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
1294 fTimeMax = nbtimebin;
e4db522f 1295
1296 //baseline = rawStream->GetCommonAdditive(); // common additive baseline
1297
3a0f6479 1298 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
3a0f6479 1299 Int_t *signal = rawStream->GetSignals(); // current ADC signal
e4db522f 1300 Int_t col = (rawStream->GetCol())%18; // current COL MCM
1301
1302 if((col < 0) || (col >= 21)) return 0;
1303 if((imcm>=16) || (imcm < 0)) return 0;
1304
3a0f6479 1305 Int_t fin = TMath::Min(fTimeMax,(iTimeBin+3));
1306 Int_t n = 0;
1307 for(Int_t itime = iTimeBin; itime < fin; itime++){
e4db522f 1308 if(signal[n]> (baseline+3)) phvalue[col][itime] = signal[n];
3a0f6479 1309 n++;
1310 }
1311 }
e4db522f 1312
3a0f6479 1313 // fill the last one
1314 if(fDetectorPreviousTrack != -1){
e4db522f 1315
1316 // take the mean values and check the first time bin
1317 for(Int_t j = 0; j < 21; j++){
1318 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1319 else mean[j] = 0.0;
1320 if(phvalue[j][0] > 200.0) first[j] = 1;
1321 else first[j] = 0;
1322 }
1323
1324 // select
1325 for(Int_t j = 1; j < 20; j++){
1326 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])){
1327 select[j] = 1;
1328 }
1329 else select[j] = 0;
1330 }
1331
1332 // fill
1333 for(Int_t j = 1; j < 20; j++){
1334 if(select[j] == 1){
1335 withInput = 2;
1336 for(Int_t k = 0; k < fTimeMax; k++){
1337 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1338 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1339 }
1340 else{
1341 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1342 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1343 }
1344 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,(3*baseline),fTimeMax);
1345 }
1346 }
1347 }
1348 }
1349
1350 // reset
1351 for(Int_t k = 0; k < 36; k++){
1352 for(Int_t j = 0; j < 21; j++){
1353 phvalue[j][k] = baseline;
170c35f1 1354 }
1355 }
1356 }
3a0f6479 1357
170c35f1 1358 }
3a0f6479 1359 else{
1360
1361 while (rawStream->Next()) {
1362
1363 Int_t idetector = rawStream->GetDet(); // current detector
e4db522f 1364 Int_t imcm = rawStream->GetMCM(); // current MCM
1365 Int_t irob = rawStream->GetROB(); // current ROB
1366
1367 if(((fMCMPrevious != imcm) || (fDetectorPreviousTrack != idetector) || (fROBPrevious != irob)) && (fDetectorPreviousTrack != -1)){
1368
1369 // take the mean values and check the first time bin
1370 for(Int_t j = 0; j < 21; j++){
1371 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1372 else mean[j] = 0.0;
1373 if(phvalue[j][0] > 200.0) first[j] = 1;
1374 else first[j] = 0;
1375 }
1376
1377 // select
1378 for(Int_t j = 1; j < 20; j++){
1379 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])){
1380 select[j] = 1;
1381 }
1382 else select[j] = 0;
1383 }
1384
1385 // fill
1386 for(Int_t j = 1; j < 20; j++){
1387 if(select[j] == 1){
1388 withInput = 2;
1389 for(Int_t k = 0; k < fTimeMax; k++){
1390 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1391 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1392 }
1393 else{
1394 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1395 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1396 }
1397 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1398 }
1399 }
1400 }
1401 }
1402
1403 // reset
1404 for(Int_t k = 0; k < 36; k++){
1405 for(Int_t j = 0; j < 21; j++){
1406 phvalue[j][k] = baseline;
3a0f6479 1407 }
1408 }
1409 }
e4db522f 1410
3a0f6479 1411 fDetectorPreviousTrack = idetector;
e4db522f 1412 fMCMPrevious = imcm;
1413 fROBPrevious = irob;
1414
1415
1416
1417 //baseline = rawStream->GetCommonAdditive(); // common baseline
1418
3a0f6479 1419 nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
1420 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
3a0f6479 1421 Int_t *signal = rawStream->GetSignals(); // current ADC signal
e4db522f 1422 Int_t col = (rawStream->GetCol())%18; // current COL MCM
1423
3a0f6479 1424 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
1425 Int_t n = 0;
e4db522f 1426
1427 if((col < 0) || (col >= 21)) return 0;
1428 if((imcm>=16) || (imcm < 0)) return 0;
1429
3a0f6479 1430 for(Int_t itime = iTimeBin; itime < fin; itime++){
e4db522f 1431 if(signal[n]>13) phvalue[col][itime] = signal[n];
3a0f6479 1432 n++;
1433 }
1434 }
1435
1436 // fill the last one
1437 if(fDetectorPreviousTrack != -1){
e4db522f 1438
1439 // take the mean values and check the first time bin
1440 for(Int_t j = 0; j < 21; j++){
1441 if(TMath::RMS(fTimeMax,phvalue[j]) != 0.0) mean[j] = TMath::Mean(fTimeMax,phvalue[j]);
1442 else mean[j] = 0.0;
1443 if(phvalue[j][0] > 200.0) first[j] = 1;
1444 else first[j] = 0;
1445 }
1446
1447 // select
1448 for(Int_t j = 1; j < 20; j++){
1449 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])){
1450 select[j] = 1;
1451 }
1452 else select[j] = 0;
1453 }
1454
1455 // fill
1456 for(Int_t j = 1; j < 20; j++){
1457 if(select[j] == 1){
1458 withInput = 2;
1459 for(Int_t k = 0; k < fTimeMax; k++){
1460 if((phvalue[j][k] >= phvalue[j-1][k]) && (phvalue[j][k] >= phvalue[j+1][k])){
1461 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j-1][k]+phvalue[j][k]+phvalue[j+1][k]),fTimeMax);
1462 }
1463 else{
1464 if((j < 19) && (phvalue[j+1][k] >= phvalue[j][k]) && (phvalue[j+1][k] >= phvalue[j+2][k])){
1465 UpdateDAQ(fDetectorPreviousTrack,0,0,k,(phvalue[j][k]+phvalue[j+1][k]+phvalue[j+2][k]),fTimeMax);
1466 }
1467 else UpdateDAQ(fDetectorPreviousTrack,0,0,k,3*baseline,fTimeMax);
1468 }
1469 }
1470 }
1471 }
1472
1473 // reset
1474 for(Int_t k = 0; k < 36; k++){
1475 for(Int_t j = 0; j < 21; j++){
1476 phvalue[j][k] = baseline;
170c35f1 1477 }
1478 }
3a0f6479 1479 }
170c35f1 1480 }
1481
1482 return withInput;
e4db522f 1483
170c35f1 1484}
1485//_____________________________________________________________________
3a0f6479 1486Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
170c35f1 1487{
1488 //
1489 // Event processing loop - AliRawReader
1490 //
1491
1492
d4232bb6 1493 AliTRDrawStreamTB rawStream(rawReader);
170c35f1 1494
1495 rawReader->Select("TRD");
1496
3a0f6479 1497 return ProcessEventDAQ(&rawStream, nocheck);
170c35f1 1498}
1499//_________________________________________________________________________
3a0f6479 1500Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
170c35f1 1501#ifdef ALI_DATE
3a0f6479 1502 eventHeaderStruct *event,
1503 Bool_t nocheck
170c35f1 1504#else
3a0f6479 1505 eventHeaderStruct* /*event*/,
1506 Bool_t /*nocheck*/
170c35f1 1507
1508#endif
1509 )
1510{
1511 //
1512 // process date event
1513 //
1514#ifdef ALI_DATE
1515 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
3a0f6479 1516 Int_t result=ProcessEventDAQ(rawReader, nocheck);
170c35f1 1517 delete rawReader;
1518 return result;
1519#else
1520 Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1521 return 0;
1522#endif
1523
1524}
1525//____________Online trackling in AliTRDtrigger________________________________
1526Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1527{
1528 //
1529 // For the DAQ
1530 // Fill a simple average pulse height
1531 //
1532
1533 // Localisation of the Xbins involved
1534 //LocalisationDetectorXbins(det);
1535
1536 // Row and position in the pad groups
1537 //Int_t posr = 0;
1538 //if (fCalibraMode->GetNnZ(1) != 0) {
1539 // posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1540 //}
1541
1542 // Col of cluster and position in the pad groups
1543 //Int_t posc = 0;
1544 //if (fCalibraMode->GetNnRphi(1) != 0) {
1545 // posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1546 //}
1547
1548 //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1549
3a0f6479 1550 ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
170c35f1 1551
1552 return kTRUE;
1553
1554}
170c35f1 1555//____________Write_____________________________________________________
1556//_____________________________________________________________________
1557void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1558{
3a0f6479 1559 //
1560 // Write infos to file
1561 //
1562
170c35f1 1563 //For debugging
1564 if ( fDebugStreamer ) {
1565 delete fDebugStreamer;
1566 fDebugStreamer = 0x0;
55a288e5 1567 }
170c35f1 1568
55a288e5 1569 AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
170c35f1 1570 ,fNumberTrack
1571 ,fNumberUsedCh[0]
1572 ,fNumberUsedCh[1]
1573 ,fNumberUsedPh[0]
1574 ,fNumberUsedPh[1]));
1575
1576 TDirectory *backup = gDirectory;
1577 TString option;
1578
1579 if ( append )
1580 option = "update";
1581 else
1582 option = "recreate";
1583
1584 TFile f(filename,option.Data());
55a288e5 1585
1586 TStopwatch stopwatch;
1587 stopwatch.Start();
3a0f6479 1588 if(fVector2d) {
1589 f.WriteTObject(fCalibraVector);
1590 }
55a288e5 1591
170c35f1 1592 if (fCH2dOn ) {
55a288e5 1593 if (fHisto2d) {
170c35f1 1594 f.WriteTObject(fCH2d);
55a288e5 1595 }
55a288e5 1596 }
170c35f1 1597 if (fPH2dOn ) {
55a288e5 1598 if (fHisto2d) {
170c35f1 1599 f.WriteTObject(fPH2d);
55a288e5 1600 }
55a288e5 1601 }
170c35f1 1602 if (fPRF2dOn) {
55a288e5 1603 if (fHisto2d) {
170c35f1 1604 f.WriteTObject(fPRF2d);
55a288e5 1605 }
55a288e5 1606 }
3a0f6479 1607 if(fLinearFitterOn){
1608 AnalyseLinearFitter();
1609 f.WriteTObject(fLinearVdriftFit);
170c35f1 1610 }
3a0f6479 1611
170c35f1 1612 f.Close();
55a288e5 1613
170c35f1 1614 if ( backup ) backup->cd();
55a288e5 1615
170c35f1 1616 AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1617 ,stopwatch.RealTime(),stopwatch.CpuTime()));
55a288e5 1618}
170c35f1 1619//___________________________________________probe the histos__________________________________________________
55a288e5 1620Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1621{
1622 //
1623 // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1624 // debug mode with 2 for TH2I and 3 for TProfile2D
1625 // It gives a pointer to a Double_t[7] with the info following...
1626 // [0] : number of calibration groups with entries
1627 // [1] : minimal number of entries found
1628 // [2] : calibration group number of the min
1629 // [3] : maximal number of entries found
1630 // [4] : calibration group number of the max
1631 // [5] : mean number of entries found
3a0f6479 1632 // [6] : mean relative error
55a288e5 1633 //
1634
1635 Double_t *info = new Double_t[7];
1636
1637 // Number of Xbins (detectors or groups of pads)
170c35f1 1638 Int_t nbins = h->GetNbinsY(); //number of calibration groups
1639 Int_t nxbins = h->GetNbinsX(); //number of bins per histo
55a288e5 1640
1641 // Initialise
1642 Double_t nbwe = 0; //number of calibration groups with entries
1643 Double_t minentries = 0; //minimal number of entries found
1644 Double_t maxentries = 0; //maximal number of entries found
1645 Double_t placemin = 0; //calibration group number of the min
1646 Double_t placemax = -1; //calibration group number of the max
1647 Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1648 Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1649
1650 Double_t counter = 0;
1651
1652 //Debug
0c349049 1653 TH1F *nbEntries = 0x0;//distribution of the number of entries
1654 TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
1655 TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
55a288e5 1656
1657 // Beginning of the loop over the calibration groups
1658 for (Int_t idect = 0; idect < nbins; idect++) {
1659
170c35f1 1660 TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
55a288e5 1661 projch->SetDirectory(0);
1662
1663 // Number of entries for this calibration group
1664 Double_t nentries = 0.0;
1665 if((i%2) == 0){
170c35f1 1666 for (Int_t k = 0; k < nxbins; k++) {
1667 nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
55a288e5 1668 }
1669 }
1670 else{
170c35f1 1671 for (Int_t k = 0; k < nxbins; k++) {
1672 nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1673 if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1674 meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
55a288e5 1675 counter++;
1676 }
1677 }
1678 }
1679
1680 //Debug
1681 if(i > 1){
0c349049 1682 if((!((Bool_t)nbEntries)) && (nentries > 0)){
1683 nbEntries = new TH1F("Number of entries","Number of entries"
55a288e5 1684 ,100,(Int_t)nentries/2,nentries*2);
0c349049 1685 nbEntries->SetDirectory(0);
1686 nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
55a288e5 1687 ,nbins,0,nbins);
0c349049 1688 nbEntriesPerGroup->SetDirectory(0);
1689 nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
55a288e5 1690 ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
0c349049 1691 nbEntriesPerSp->SetDirectory(0);
55a288e5 1692 }
0c349049 1693 if(nbEntries){
1694 if(nentries > 0) nbEntries->Fill(nentries);
1695 nbEntriesPerGroup->Fill(idect+0.5,nentries);
1696 nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
55a288e5 1697 }
1698 }
1699
1700 //min amd max
1701 if(nentries > maxentries){
1702 maxentries = nentries;
1703 placemax = idect;
1704 }
1705 if(idect == 0) {
1706 minentries = nentries;
1707 }
1708 if(nentries < minentries){
1709 minentries = nentries;
1710 placemin = idect;
1711 }
1712 //nbwe
1713 if(nentries > 0) {
1714 nbwe++;
1715 meanstats += nentries;
1716 }
55a288e5 1717 }//calibration groups loop
1718
1719 if(nbwe > 0) meanstats /= nbwe;
1720 if(counter > 0) meanrelativerror /= counter;
1721
1722 AliInfo(Form("There are %f calibration groups with entries",nbwe));
1723 AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1724 AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1725 AliInfo(Form("The mean number of entries is %f",meanstats));
1726 if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1727
1728 info[0] = nbwe;
1729 info[1] = minentries;
1730 info[2] = placemin;
1731 info[3] = maxentries;
1732 info[4] = placemax;
1733 info[5] = meanstats;
1734 info[6] = meanrelativerror;
1735
1736 if(i > 1){
1737 gStyle->SetPalette(1);
1738 gStyle->SetOptStat(1111);
1739 gStyle->SetPadBorderMode(0);
1740 gStyle->SetCanvasColor(10);
1741 gStyle->SetPadLeftMargin(0.13);
1742 gStyle->SetPadRightMargin(0.01);
1743 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1744 stat->Divide(2,1);
1745 stat->cd(1);
0c349049 1746 nbEntries->Draw("");
55a288e5 1747 stat->cd(2);
0c349049 1748 nbEntriesPerSp->SetStats(0);
1749 nbEntriesPerSp->Draw("");
55a288e5 1750 TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1751 stat1->cd();
0c349049 1752 nbEntriesPerGroup->SetStats(0);
1753 nbEntriesPerGroup->Draw("");
55a288e5 1754 }
1755
1756 return info;
1757
1758}
170c35f1 1759//____________________________________________________________________________
1760Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1761{
1762 //
1763 // Return a Int_t[4] with:
1764 // 0 Mean number of entries
1765 // 1 median of number of entries
1766 // 2 rms of number of entries
1767 // 3 number of group with entries
1768 //
1769
1770 Double_t *stat = new Double_t[4];
1771 stat[3] = 0.0;
1772
1773 Int_t nbofgroups = CalculateTotalNumberOfBins(0);
1774 Double_t *weight = new Double_t[nbofgroups];
1775 Int_t *nonul = new Int_t[nbofgroups];
1776
1777 for(Int_t k = 0; k < nbofgroups; k++){
1778 if(fEntriesCH[k] > 0) {
1779 weight[k] = 1.0;
1780 nonul[(Int_t)stat[3]] = fEntriesCH[k];
1781 stat[3]++;
1782 }
1783 else weight[k] = 0.0;
1784 }
1785 stat[0] = TMath::Mean(nbofgroups,fEntriesCH,weight);
1786 stat[1] = TMath::Median(nbofgroups,fEntriesCH,weight);
1787 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1788
1789 return stat;
1790
1791}
1792//____________________________________________________________________________
1793Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1794{
1795 //
1796 // Return a Int_t[4] with:
1797 // 0 Mean number of entries
1798 // 1 median of number of entries
1799 // 2 rms of number of entries
1800 // 3 number of group with entries
1801 //
1802
1803 Double_t *stat = new Double_t[4];
1804 stat[3] = 0.0;
1805
1806 Int_t nbofgroups = 540;
1807 Double_t *weight = new Double_t[nbofgroups];
1808 Int_t *nonul = new Int_t[nbofgroups];
1809
1810 for(Int_t k = 0; k < nbofgroups; k++){
1811 if(fEntriesLinearFitter[k] > 0) {
1812 weight[k] = 1.0;
1813 nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1814 stat[3]++;
1815 }
1816 else weight[k] = 0.0;
1817 }
1818 stat[0] = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1819 stat[1] = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1820 stat[2] = TMath::RMS((Int_t)stat[3],nonul);
1821
1822 return stat;
1823
1824}
1825//_____________________________________________________________________________
1826void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1827{
1828 //
1829 // Should be between 0 and 6
1830 //
1831
1832 if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1833 AliInfo("The number of groups must be between 0 and 6!");
1834 }
1835 else {
1836 fNgroupprf = numberGroupsPRF;
1837 }
55a288e5 1838
170c35f1 1839}
55a288e5 1840//_____________________________________________________________________________
1841void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1842{
1843 //
1844 // Set the factor that will divide the deposited charge
1845 // to fit in the histo range [0,300]
1846 //
1847
1848 if (RelativeScale > 0.0) {
1849 fRelativeScale = RelativeScale;
1850 }
1851 else {
1852 AliInfo("RelativeScale must be strict positif!");
1853 }
1854
1855}
1856
1857//_____________________________________________________________________________
1858void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1859{
1860 //
1861 // Set the mode of calibration group in the z direction for the parameter i
1862 //
1863
1864 if ((Nz >= 0) &&
1865 (Nz < 5)) {
1866 fCalibraMode->SetNz(i, Nz);
1867 }
1868 else {
1869 AliInfo("You have to choose between 0 and 4");
1870 }
1871
1872}
1873
1874//_____________________________________________________________________________
1875void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1876{
1877 //
1878 // Set the mode of calibration group in the rphi direction for the parameter i
1879 //
1880
1881 if ((Nrphi >= 0) &&
1882 (Nrphi < 7)) {
1883 fCalibraMode->SetNrphi(i ,Nrphi);
1884 }
1885 else {
1886 AliInfo("You have to choose between 0 and 6");
1887 }
1888
1889}
55a288e5 1890//____________Protected Functions______________________________________________
1891//____________Create the 2D histo to be filled online__________________________
1892//
55a288e5 1893//_____________________________________________________________________________
1894void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1895{
1896 //
170c35f1 1897 // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1898 // If fNgroupprf is zero then no binning in tnp
55a288e5 1899 //
1900
1901 TString name("Nz");
1902 name += fCalibraMode->GetNz(2);
1903 name += "Nrphi";
1904 name += fCalibraMode->GetNrphi(2);
170c35f1 1905 name += "Ngp";
1906 name += fNgroupprf;
55a288e5 1907
170c35f1 1908 if(fNgroupprf != 0){
1909
1910 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1911 ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1912 fPRF2d->SetYTitle("Det/pad groups");
1913 fPRF2d->SetXTitle("Position x/W [pad width units]");
1914 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1915 fPRF2d->SetStats(0);
1916 }
1917 else{
1918 fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1919 ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1920 fPRF2d->SetYTitle("Det/pad groups");
1921 fPRF2d->SetXTitle("Position x/W [pad width units]");
1922 fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1923 fPRF2d->SetStats(0);
1924 }
55a288e5 1925
1926}
1927
1928//_____________________________________________________________________________
1929void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1930{
1931 //
1932 // Create the 2D histos
1933 //
1934
1935 TString name("Nz");
1936 name += fCalibraMode->GetNz(1);
1937 name += "Nrphi";
1938 name += fCalibraMode->GetNrphi(1);
170c35f1 1939
55a288e5 1940 fPH2d = new TProfile2D("PH2d",(const Char_t *) name
170c35f1 1941 ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1942 ,nn,0,nn);
1943 fPH2d->SetYTitle("Det/pad groups");
1944 fPH2d->SetXTitle("time [#mus]");
55a288e5 1945 fPH2d->SetZTitle("<PH> [a.u.]");
1946 fPH2d->SetStats(0);
1947
1948}
55a288e5 1949//_____________________________________________________________________________
1950void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1951{
1952 //
1953 // Create the 2D histos
1954 //
1955
1956 TString name("Nz");
1957 name += fCalibraMode->GetNz(0);
1958 name += "Nrphi";
1959 name += fCalibraMode->GetNrphi(0);
1960
1961 fCH2d = new TH2I("CH2d",(const Char_t *) name
170c35f1 1962 ,fNumberBinCharge,0,300,nn,0,nn);
1963 fCH2d->SetYTitle("Det/pad groups");
1964 fCH2d->SetXTitle("charge deposit [a.u]");
55a288e5 1965 fCH2d->SetZTitle("counts");
1966 fCH2d->SetStats(0);
1967 fCH2d->Sumw2();
1968
1969}
3a0f6479 1970
55a288e5 1971//____________Offine tracking in the AliTRDtracker_____________________________
1972void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1973{
1974 //
1975 // For the offline tracking or mcm tracklets
1976 // This function will be called in the functions UpdateHistogram...
1977 // to fill the info of a track for the relativ gain calibration
1978 //
1979
170c35f1 1980 Int_t nb = 0; // Nombre de zones traversees
1981 Int_t fd = -1; // Premiere zone non nulle
1982 Float_t totalcharge = 0.0; // Total charge for the supermodule histo
3a0f6479 1983
1984
55a288e5 1985
1986
1987 // See if the track goes through different zones
1988 for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1989 if (fAmpTotal[k] > 0.0) {
170c35f1 1990 totalcharge += fAmpTotal[k];
55a288e5 1991 nb++;
1992 if (nb == 1) {
1993 fd = k;
1994 }
1995 }
1996 }
55a288e5 1997
3a0f6479 1998
170c35f1 1999 switch (nb)
2000 {
2001 case 1:
55a288e5 2002 fNumberUsedCh[0]++;
170c35f1 2003 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2004 if (fHisto2d) {
170c35f1 2005 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2006 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2007 }
2008 if (fVector2d) {
3a0f6479 2009 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2010 }
170c35f1 2011 break;
2012 case 2:
55a288e5 2013 if ((fAmpTotal[fd] > 0.0) &&
170c35f1 2014 (fAmpTotal[fd+1] > 0.0)) {
55a288e5 2015 // One of the two very big
2016 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2017 if (fHisto2d) {
170c35f1 2018 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2019 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2020 }
2021 if (fVector2d) {
3a0f6479 2022 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2023 }
2024 fNumberUsedCh[1]++;
170c35f1 2025 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2026 }
2027 if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2028 if (fHisto2d) {
170c35f1 2029 FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
2030 //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
55a288e5 2031 }
2032 if (fVector2d) {
3a0f6479 2033 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
55a288e5 2034 }
2035 fNumberUsedCh[1]++;
170c35f1 2036 fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
55a288e5 2037 }
2038 }
55a288e5 2039 if (fCalibraMode->GetNfragZ(0) > 1) {
2040 if (fAmpTotal[fd] > 0.0) {
2041 if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2042 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2043 // One of the two very big
2044 if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2045 if (fHisto2d) {
170c35f1 2046 FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
2047 //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
55a288e5 2048 }
2049 if (fVector2d) {
3a0f6479 2050 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
55a288e5 2051 }
2052 fNumberUsedCh[1]++;
170c35f1 2053 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
55a288e5 2054 }
2055 if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2056 if (fHisto2d) {
170c35f1 2057 FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
2058 //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
55a288e5 2059 }
2060 fNumberUsedCh[1]++;
170c35f1 2061 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
55a288e5 2062 if (fVector2d) {
3a0f6479 2063 fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
55a288e5 2064 }
2065 }
2066 }
2067 }
2068 }
2069 }
170c35f1 2070 break;
2071 default: break;
55a288e5 2072 }
55a288e5 2073}
55a288e5 2074//____________Offine tracking in the AliTRDtracker_____________________________
2075void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2076{
2077 //
2078 // For the offline tracking or mcm tracklets
2079 // This function will be called in the functions UpdateHistogram...
2080 // to fill the info of a track for the drift velocity calibration
2081 //
2082
2083 Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
2084 Int_t fd1 = -1; // Premiere zone non nulle
2085 Int_t fd2 = -1; // Deuxieme zone non nulle
2086 Int_t k1 = -1; // Debut de la premiere zone
2087 Int_t k2 = -1; // Debut de la seconde zone
2088
2089 // See if the track goes through different zones
2090 for (Int_t k = 0; k < fTimeMax; k++) {
2091 if (fPHValue[k] > 0.0) {
2092 if (fd1 == -1) {
2093 fd1 = fPHPlace[k];
2094 k1 = k;
2095 }
2096 if (fPHPlace[k] != fd1) {
2097 if (fd2 == -1) {
2098 k2 = k;
2099 fd2 = fPHPlace[k];
2100 nb = 2;
2101 }
2102 if (fPHPlace[k] != fd2) {
2103 nb = 3;
2104 }
2105 }
2106 }
2107 }
170c35f1 2108
3a0f6479 2109
170c35f1 2110 switch(nb)
2111 {
2112 case 1:
2113 fNumberUsedPh[0]++;
2114 for (Int_t i = 0; i < fTimeMax; i++) {
2115 if (fHisto2d) {
2116 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2117 }
2118 if (fVector2d) {
3a0f6479 2119 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
55a288e5 2120 }
2121 }
170c35f1 2122 break;
2123 case 2:
2124 if ((fd1 == fd2+1) ||
2125 (fd2 == fd1+1)) {
2126 // One of the two fast all the think
2127 if (k2 > (k1+fDifference)) {
2128 //we choose to fill the fd1 with all the values
2129 fNumberUsedPh[1]++;
2130 for (Int_t i = 0; i < fTimeMax; i++) {
2131 if (fHisto2d) {
2132 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2133 }
2134 if (fVector2d) {
3a0f6479 2135 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2136 }
55a288e5 2137 }
170c35f1 2138 }
2139 if ((k2+fDifference) < fTimeMax) {
2140 //we choose to fill the fd2 with all the values
2141 fNumberUsedPh[1]++;
2142 for (Int_t i = 0; i < fTimeMax; i++) {
2143 if (fHisto2d) {
2144 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2145 }
55a288e5 2146 if (fVector2d) {
3a0f6479 2147 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
55a288e5 2148 }
170c35f1 2149 }
55a288e5 2150 }
2151 }
170c35f1 2152 // Two zones voisines sinon rien!
2153 if (fCalibraMode->GetNfragZ(1) > 1) {
2154 // Case 2
2155 if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2156 if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2157 // One of the two fast all the think
2158 if (k2 > (k1+fDifference)) {
2159 //we choose to fill the fd1 with all the values
2160 fNumberUsedPh[1]++;
2161 for (Int_t i = 0; i < fTimeMax; i++) {
2162 if (fHisto2d) {
2163 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2164 }
2165 if (fVector2d) {
3a0f6479 2166 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2167 }
55a288e5 2168 }
2169 }
170c35f1 2170 if ((k2+fDifference) < fTimeMax) {
2171 //we choose to fill the fd2 with all the values
2172 fNumberUsedPh[1]++;
2173 for (Int_t i = 0; i < fTimeMax; i++) {
2174 if (fHisto2d) {
2175 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2176 }
2177 if (fVector2d) {
3a0f6479 2178 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
170c35f1 2179 }
55a288e5 2180 }
2181 }
2182 }
2183 }
170c35f1 2184 // Two zones voisines sinon rien!
2185 // Case 3
2186 if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2187 if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2188 // One of the two fast all the think
2189 if (k2 > (k1 + fDifference)) {
2190 //we choose to fill the fd1 with all the values
2191 fNumberUsedPh[1]++;
2192 for (Int_t i = 0; i < fTimeMax; i++) {
2193 if (fHisto2d) {
2194 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2195 }
2196 if (fVector2d) {
3a0f6479 2197 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
170c35f1 2198 }
55a288e5 2199 }
2200 }
170c35f1 2201 if ((k2+fDifference) < fTimeMax) {
2202 //we choose to fill the fd2 with all the values
2203 fNumberUsedPh[1]++;
2204 for (Int_t i = 0; i < fTimeMax; i++) {
2205 if (fHisto2d) {
2206 fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2207 }
2208 if (fVector2d) {
3a0f6479 2209 fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
170c35f1 2210 }
55a288e5 2211 }
2212 }
2213 }
2214 }
2215 }
170c35f1 2216 break;
2217 default: break;
2218 }
55a288e5 2219}
170c35f1 2220//____________Offine tracking in the AliTRDtracker_____________________________
2221Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
55a288e5 2222{
2223 //
170c35f1 2224 // For the offline tracking
2225 // This function will be called in the functions UpdateHistogram...
2226 // to fill the find the parameter P1 of a track for the drift velocity calibration
55a288e5 2227 //
3a0f6479 2228
2229
170c35f1 2230 //Number of points: if less than 3 return kFALSE
0c349049 2231 Int_t npoints = fListClusters->GetEntriesFast();
2232 if(npoints <= 2) return kFALSE;
170c35f1 2233
2234 //Variables
2235 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2236 Double_t snp = 0.0; // sin angle in the plan yx track
2237 Double_t y = 0.0; // y clusters in the middle of the chamber
2238 Double_t z = 0.0; // z cluster in the middle of the chamber
2239 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2240 Double_t tnp = 0.0; // tan angle in the plan xy track
2241 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2242 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2243 Double_t pointError = 0.0; // error after straight line fit
2244 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
2245 Int_t snpright = 1; // if we took in the middle snp
2246 Int_t crossrow = 0; // if it crosses a pad row
2247 Double_t tiltingangle = 0; // tiltingangle of the pad
2248 Float_t dzdx = 0; // dz/dx now from dz/dl
2249 Int_t nbli = 0; // number linear fitter points
f162af62 2250 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
170c35f1 2251
2252 linearFitterTracklet.StoreData(kFALSE);
2253 linearFitterTracklet.ClearPoints();
2254
2255 //if more than one row
2256 Int_t rowp = -1; // if it crosses a pad row
2257
2258 //tiltingangle
2259 tiltingangle = padplane->GetTiltingAngle();
2260 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2261
2262 //Fill with points
0c349049 2263 for(Int_t k = 0; k < npoints; k++){
170c35f1 2264
2265 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2266 Double_t ycluster = cl->GetY();
2267 Int_t time = cl->GetLocalTimeBin();
2268 Double_t timeis = time/fSf;
2269 //See if cross two pad rows
2270 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2271 if(k==0) rowp = row;
2272 if(row != rowp) crossrow = 1;
2273 //Take in the middle of the chamber
3a0f6479 2274 //FollowBack
170c35f1 2275 if(time > (Int_t) 10) {
3a0f6479 2276 //Follow
2277 //if(time < (Int_t) 11) {
170c35f1 2278 z = cl->GetZ();
2279 y = cl->GetY();
2280 snp = fPar2[time];
2281 tgl = fPar3[time];
2282 }
2283 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2284 nbli++;
55a288e5 2285 }
3a0f6479 2286 //FollowBack
2287 if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
2288 //Follow
2289 //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
170c35f1 2290 if(nbli <= 2) return kFALSE;
2291
2292 // Do the straight line fit now
2293 TVectorD pars;
2294 linearFitterTracklet.Eval();
2295 linearFitterTracklet.GetParameters(pars);
2296 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2297 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2298 dydt = pars[1];
2299
2300 if( TMath::Abs(snp) < 1.){
2301 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2302 }
2303 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2304
2305 if(fDebugLevel > 0){
2306 if ( !fDebugStreamer ) {
2307 //debug stream
2308 TDirectory *backup = gDirectory;
3a0f6479 2309 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
170c35f1 2310 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2311 }
2312
2313 (* fDebugStreamer) << "VDRIFT0"<<
0c349049 2314 "npoints="<<npoints<<
170c35f1 2315 "\n";
2316
2317
2318 (* fDebugStreamer) << "VDRIFT"<<
2319 "snpright="<<snpright<<
0c349049 2320 "npoints="<<npoints<<
170c35f1 2321 "nbli="<<nbli<<
2322 "detector="<<detector<<
2323 "snp="<<snp<<
2324 "tnp="<<tnp<<
2325 "tgl="<<tgl<<
2326 "tnt="<<tnt<<
2327 "y="<<y<<
2328 "z="<<z<<
2329 "dydt="<<dydt<<
2330 "dzdx="<<dzdx<<
2331 "crossrow="<<crossrow<<
2332 "errorpar="<<errorpar<<
2333 "pointError="<<pointError<<
2334 "\n";
55a288e5 2335
55a288e5 2336 }
2337
0c349049 2338 if(npoints < fNumberClusters) return kFALSE;
170c35f1 2339 if(snpright == 0) return kFALSE;
2340 if(pointError >= 0.1) return kFALSE;
2341 if(crossrow == 1) return kFALSE;
2342
2343 if(fLinearFitterOn){
2344 //Add to the linear fitter of the detector
2345 if( TMath::Abs(snp) < 1.){
2346 Double_t x = tnp-dzdx*tnt;
2347 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2348 if(fLinearFitterDebugOn) {
3a0f6479 2349 fLinearVdriftFit->Update(detector,x,pars[1]);
170c35f1 2350 }
2351 fEntriesLinearFitter[detector]++;
2352 }
2353 }
2354 //AliInfo("End of FindP1TrackPH with success!")
55a288e5 2355 return kTRUE;
2356
2357}
170c35f1 2358//____________Offine tracking in the AliTRDtracker_____________________________
2359Bool_t AliTRDCalibraFillHisto::HandlePRF()
2360{
2361 //
2362 // For the offline tracking
2363 // Fit the tracklet with a line and take the position as reference for the PRF
2364 //
3a0f6479 2365
170c35f1 2366 //Number of points
0c349049 2367 Int_t npoints = fListClusters->GetEntriesFast(); // number of total points
2368 Int_t nb3pc = 0; // number of three pads clusters used for fit
170c35f1 2369 Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
2370
2371
2372 // To see the difference due to the fit
2373 Double_t *padPositions;
0c349049 2374 padPositions = new Double_t[npoints];
2375 for(Int_t k = 0; k < npoints; k++){
170c35f1 2376 padPositions[k] = 0.0;
2377 }
2378
2379
2380 //Find the position by a fit
2381 TLinearFitter fitter(2,"pol1");
2382 fitter.StoreData(kFALSE);
2383 fitter.ClearPoints();
0c349049 2384 for(Int_t k = 0; k < npoints; k++){
170c35f1 2385 //Take the cluster
2386 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2387 Short_t *signals = cl->GetSignals();
2388 Double_t time = cl->GetLocalTimeBin();
2389 //Calculate x if possible
2390 Float_t xcenter = 0.0;
2391 Bool_t echec = kTRUE;
2392 if((time<=7) || (time>=21)) continue;
2393 // Center 3 balanced: position with the center of the pad
2394 if ((((Float_t) signals[3]) > 0.0) &&
2395 (((Float_t) signals[2]) > 0.0) &&
2396 (((Float_t) signals[4]) > 0.0)) {
2397 echec = kFALSE;
2398 // Security if the denomiateur is 0
2399 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2400 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2401 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2402 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2403 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2404 }
2405 else {
2406 echec = kTRUE;
2407 }
2408 }
2409 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2410 if(echec) continue;
2411 //if no echec: calculate with the position of the pad
2412 // Position of the cluster
af26ce80 2413 Double_t padPosition = xcenter + cl->GetPadCol();
170c35f1 2414 padPositions[k] = padPosition;
0c349049 2415 nb3pc++;
170c35f1 2416 fitter.AddPoint(&time, padPosition,1);
2417 }//clusters loop
2418
0c349049 2419 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2420 if(nb3pc < 3) return kFALSE;
170c35f1 2421 fitter.Eval();
2422 TVectorD line(2);
2423 fitter.GetParameters(line);
2424 Float_t pointError = -1.0;
0c349049 2425 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
170c35f1 2426
2427
2428 // Now fill the PRF
0c349049 2429 for(Int_t k = 0; k < npoints; k++){
170c35f1 2430 //Take the cluster
2431 AliTRDcluster *cl = (AliTRDcluster *) fListClusters->At(k);
2432 Short_t *signals = cl->GetSignals(); // signal
2433 Double_t time = cl->GetLocalTimeBin(); // time bin
2434 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
af26ce80 2435 Float_t padPos = cl->GetPadCol(); // middle pad
170c35f1 2436 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2437 Float_t ycenter = 0.0; // relative center charge
2438 Float_t ymin = 0.0; // relative left charge
2439 Float_t ymax = 0.0; // relative right charge
2440 Double_t tgl = fPar3[(Int_t)time]; // dz/dl and not dz/dx
2441 Double_t pt = fPar4[(Int_t)time]; // pt
2442 Float_t dzdx = 0.0; // dzdx
2443
2444
2445 //Requiere simply two pads clusters at least
2446 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2447 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2448 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2449 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2450 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2451 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2452 }
2453
2454 //calibration group
2455 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2456 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
3a0f6479 2457 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
170c35f1 2458 Double_t snp = fPar2[(Int_t)time]; // sin angle in xy plan
2459 Float_t xcl = cl->GetY(); // y cluster
2460 Float_t qcl = cl->GetQ(); // charge cluster
2461 Int_t plane = GetPlane(detector); // plane
2462 Int_t chamber = GetChamber(detector); // chamber
2463 Double_t xdiff = dpad; // reconstructed position constant
2464 Double_t x = dpad; // reconstructed position moved
0c349049 2465 Float_t ep = pointError; // error of fit
170c35f1 2466 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2467 Float_t signal3 = (Float_t)signals[3]; // signal
2468 Float_t signal2 = (Float_t)signals[2]; // signal
2469 Float_t signal4 = (Float_t)signals[4]; // signal
2470 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2471 Float_t tnp = 0.0;
2472 if(TMath::Abs(snp) < 1.0){
2473 tnp = snp / (TMath::Sqrt(1-snp*snp));
2474 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2475 }
2476
55a288e5 2477
170c35f1 2478 if(fDebugLevel > 0){
2479 if ( !fDebugStreamer ) {
2480 //debug stream
2481 TDirectory *backup = gDirectory;
3a0f6479 2482 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
170c35f1 2483 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2484 }
2485
2486 (* fDebugStreamer) << "PRF0"<<
2487 "caligroup="<<caligroup<<
2488 "detector="<<detector<<
2489 "plane="<<plane<<
2490 "chamber="<<chamber<<
0c349049 2491 "npoints="<<npoints<<
2492 "Np="<<nb3pc<<
2493 "ep="<<ep<<
170c35f1 2494 "snp="<<snp<<
2495 "tnp="<<tnp<<
2496 "tgl="<<tgl<<
2497 "dzdx="<<dzdx<<
2498 "pt="<<pt<<
2499 "padPos="<<padPos<<
2500 "padPosition="<<padPositions[k]<<
2501 "padPosTracklet="<<padPosTracklet<<
2502 "x="<<x<<
2503 "ycenter="<<ycenter<<
2504 "ymin="<<ymin<<
2505 "ymax="<<ymax<<
2506 "xcl="<<xcl<<
2507 "qcl="<<qcl<<
2508 "signal1="<<signal1<<
2509 "signal2="<<signal2<<
2510 "signal3="<<signal3<<
2511 "signal4="<<signal4<<
2512 "signal5="<<signal5<<
2513 "time="<<time<<
2514 "\n";
2515 x = xdiff;
2516 Int_t type=0;
2517 Float_t y = ycenter;
2518 (* fDebugStreamer) << "PRFALL"<<
2519 "caligroup="<<caligroup<<
2520 "detector="<<detector<<
2521 "plane="<<plane<<
2522 "chamber="<<chamber<<
0c349049 2523 "npoints="<<npoints<<
2524 "Np="<<nb3pc<<
2525 "ep="<<ep<<
170c35f1 2526 "type="<<type<<
2527 "snp="<<snp<<
2528 "tnp="<<tnp<<
2529 "tgl="<<tgl<<
2530 "dzdx="<<dzdx<<
2531 "pt="<<pt<<
2532 "padPos="<<padPos<<
2533 "padPosition="<<padPositions[k]<<
2534 "padPosTracklet="<<padPosTracklet<<
2535 "x="<<x<<
2536 "y="<<y<<
2537 "xcl="<<xcl<<
2538 "qcl="<<qcl<<
2539 "signal1="<<signal1<<
2540 "signal2="<<signal2<<
2541 "signal3="<<signal3<<
2542 "signal4="<<signal4<<
2543 "signal5="<<signal5<<
2544 "time="<<time<<
2545 "\n";
2546 x=-(xdiff+1);
2547 y = ymin;
2548 type=-1;
2549 (* fDebugStreamer) << "PRFALL"<<
2550 "caligroup="<<caligroup<<
2551 "detector="<<detector<<
2552 "plane="<<plane<<
2553 "chamber="<<chamber<<
0c349049 2554 "npoints="<<npoints<<
2555 "Np="<<nb3pc<<
2556 "ep="<<ep<<
170c35f1 2557 "type="<<type<<
2558 "snp="<<snp<<
2559 "tnp="<<tnp<<
2560 "tgl="<<tgl<<
2561 "dzdx="<<dzdx<<
2562 "pt="<<pt<<
2563 "padPos="<<padPos<<
2564 "padPosition="<<padPositions[k]<<
2565 "padPosTracklet="<<padPosTracklet<<
2566 "x="<<x<<
2567 "y="<<y<<
2568 "xcl="<<xcl<<
2569 "qcl="<<qcl<<
2570 "signal1="<<signal1<<
2571 "signal2="<<signal2<<
2572 "signal3="<<signal3<<
2573 "signal4="<<signal4<<
2574 "signal5="<<signal5<<
2575 "time="<<time<<
2576 "\n";
2577 x=1-xdiff;
2578 y = ymax;
2579 type=1;
2580
2581 (* fDebugStreamer) << "PRFALL"<<
2582 "caligroup="<<caligroup<<
2583 "detector="<<detector<<
2584 "plane="<<plane<<
2585 "chamber="<<chamber<<
0c349049 2586 "npoints="<<npoints<<
2587 "Np="<<nb3pc<<
2588 "ep="<<ep<<
170c35f1 2589 "type="<<type<<
2590 "snp="<<snp<<
2591 "tnp="<<tnp<<
2592 "tgl="<<tgl<<
2593 "dzdx="<<dzdx<<
2594 "pt="<<pt<<
2595 "padPos="<<padPos<<
2596 "padPosition="<<padPositions[k]<<
2597 "padPosTracklet="<<padPosTracklet<<
2598 "x="<<x<<
2599 "y="<<y<<
2600 "xcl="<<xcl<<
2601 "qcl="<<qcl<<
2602 "signal1="<<signal1<<
2603 "signal2="<<signal2<<
2604 "signal3="<<signal3<<
2605 "signal4="<<signal4<<
2606 "signal5="<<signal5<<
2607 "time="<<time<<
2608 "\n";
2609
2610 }
2611
2612 // some cuts
0c349049 2613 if(npoints < fNumberClusters) continue;
2614 if(nb3pc <= 5) continue;
170c35f1 2615 if((time >= 21) || (time < 7)) continue;
2616 if(TMath::Abs(snp) >= 1.0) continue;
2617 if(qcl < 80) continue;
2618
2619 Bool_t echec = kFALSE;
2620 Double_t shift = 0.0;
2621 //Calculate the shift in x coresponding to this tnp
2622 if(fNgroupprf != 0.0){
2623 shift = -3.0*(fNgroupprf-1)-1.5;
2624 Double_t limithigh = -0.2*(fNgroupprf-1);
2625 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2626 else{
2627 while(tnp > limithigh){
2628 limithigh += 0.2;
2629 shift += 3.0;
2630 }
2631 }
2632 }
2633 if (fHisto2d && !echec) {
3a0f6479 2634 if(TMath::Abs(dpad) < 1.5) {
2635 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2636 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2637 }
2638 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
170c35f1 2639 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2640 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2641 }
3a0f6479 2642 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
170c35f1 2643 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2644 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2645 }
2646 }
2647 //Not equivalent anymore here!
2648 if (fVector2d && !echec) {
3a0f6479 2649 if(TMath::Abs(dpad) < 1.5) {
2650 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2651 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
170c35f1 2652 }
3a0f6479 2653 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2654 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2655 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2656 }
2657 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2658 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2659 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
170c35f1 2660 }
2661 }
2662 }
2663 return kTRUE;
2664
2665}
d0569428 2666//____________Offine tracking in the AliTRDtracker_____________________________
2667Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
55a288e5 2668{
2669 //
d0569428 2670 // For the offline tracking
2671 // This function will be called in the functions UpdateHistogram...
2672 // to fill the find the parameter P1 of a track for the drift velocity calibration
55a288e5 2673 //
2674
d0569428 2675
2676 //Number of points: if less than 3 return kFALSE
0c349049 2677 Int_t npoints = index1-index0;
2678 if(npoints <= 2) return kFALSE;
55a288e5 2679
d0569428 2680 //Variables
2681 TLinearFitter linearFitterTracklet = TLinearFitter(2,"pol1"); // TLinearFitter per tracklet
2682 Double_t snp = 0.0; // sin angle in the plan yx track
2683 Double_t y = 0.0; // y clusters in the middle of the chamber
2684 Double_t z = 0.0; // z cluster in the middle of the chamber
2685 Double_t dydt = 0.0; // dydt tracklet after straight line fit
2686 Double_t tnp = 0.0; // tan angle in the plan xy track
2687 Double_t tgl = 0.0; // dz/dl and not dz/dx!
2688 Double_t errorpar = 0.0; // error after straight line fit on dy/dt
2689 Double_t pointError = 0.0; // error after straight line fit
2690 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
2691 //Int_t snpright = 1; // if we took in the middle snp
2692 Int_t crossrow = 0; // if it crosses a pad row
2693 Double_t tiltingangle = 0; // tiltingangle of the pad
2694 Float_t dzdx = 0; // dz/dx now from dz/dl
2695 Int_t nbli = 0; // number linear fitter points
2696 AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
55a288e5 2697
d0569428 2698 linearFitterTracklet.StoreData(kFALSE);
2699 linearFitterTracklet.ClearPoints();
2700
2701 //if more than one row
2702 Int_t rowp = -1; // if it crosses a pad row
55a288e5 2703
d0569428 2704 //tiltingangle
2705 tiltingangle = padplane->GetTiltingAngle();
2706 Float_t tnt = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
55a288e5 2707
d0569428 2708 //Fill with points
0c349049 2709 for(Int_t k = 0; k < npoints; k++){
d0569428 2710
2711 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2712 Double_t ycluster = cl->GetY();
2713 Int_t time = cl->GetLocalTimeBin();
2714 Double_t timeis = time/fSf;
2715 //See if cross two pad rows
2716 Int_t row = padplane->GetPadRowNumber(cl->GetZ());
2717 if(k==0) rowp = row;
2718 if(row != rowp) crossrow = 1;
2719 //Take in the middle of the chamber
2720 //FollowBack
2721 //if(time > (Int_t) 10) {
2722 //Follow
2723 if(time < (Int_t) 11) {
2724 z = cl->GetZ();
2725 y = cl->GetY();
2726 }
2727 linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2728 nbli++;
2729 }
55a288e5 2730
d0569428 2731 // take now the snp, tnp and tgl from the track
2732 snp = t->GetSnpPlane(GetPlane(detector));
2733 tgl = t->GetTglPlane(GetPlane(detector));
170c35f1 2734
d0569428 2735 //FollowBack
2736 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0;
2737 //Follow
2738 //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0;
2739 if(nbli <= 2) return kFALSE;
170c35f1 2740
d0569428 2741 // Do the straight line fit now
2742 TVectorD pars;
2743 linearFitterTracklet.Eval();
2744 linearFitterTracklet.GetParameters(pars);
2745 pointError = TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2746 errorpar = linearFitterTracklet.GetParError(1)*pointError;
2747 dydt = pars[1];
2748
2749 if( TMath::Abs(snp) < 1.){
2750 tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2751 }
2752 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2753
2754 if(fDebugLevel > 0){
2755 if ( !fDebugStreamer ) {
2756 //debug stream
2757 TDirectory *backup = gDirectory;
2758 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2759 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2760 }
2761
2762
2763 (* fDebugStreamer) << "VDRIFT"<<
2764 //"snpright="<<snpright<<
0c349049 2765 "npoints="<<npoints<<
d0569428 2766 "nbli="<<nbli<<
2767 "detector="<<detector<<
2768 "snp="<<snp<<
2769 "tnp="<<tnp<<
2770 "tgl="<<tgl<<
2771 "tnt="<<tnt<<
2772 "y="<<y<<
2773 "z="<<z<<
2774 "dydt="<<dydt<<
2775 "dzdx="<<dzdx<<
2776 "crossrow="<<crossrow<<
2777 "errorpar="<<errorpar<<
2778 "pointError="<<pointError<<
2779 "\n";
2780
2781 }
2782
0c349049 2783 if(npoints < fNumberClusters) return kFALSE;
d0569428 2784 //if(snpright == 0) return kFALSE;
2785 if(pointError >= 0.1) return kFALSE;
2786 if(crossrow == 1) return kFALSE;
2787
2788 if(fLinearFitterOn){
2789 //Add to the linear fitter of the detector
2790 if( TMath::Abs(snp) < 1.){
2791 Double_t x = tnp-dzdx*tnt;
2792 (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2793 if(fLinearFitterDebugOn) {
2794 fLinearVdriftFit->Update(detector,x,pars[1]);
2795 }
2796 fEntriesLinearFitter[detector]++;
2797 }
2798 }
2799 //AliInfo("End of FindP1TrackPH with success!")
2800 return kTRUE;
2801
2802}
2803//____________Offine tracking in the AliTRDtracker_____________________________
2804Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2805{
2806 //
2807 // For the offline tracking
2808 // Fit the tracklet with a line and take the position as reference for the PRF
2809 //
2810
2811 //Number of points
0c349049 2812 Int_t npoints = index1-index0; // number of total points
2813 Int_t nb3pc = 0; // number of three pads clusters used for fit
d0569428 2814 Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
2815
2816
2817 // To see the difference due to the fit
2818 Double_t *padPositions;
0c349049 2819 padPositions = new Double_t[npoints];
2820 for(Int_t k = 0; k < npoints; k++){
d0569428 2821 padPositions[k] = 0.0;
2822 }
2823
2824
2825 //Find the position by a fit
2826 TLinearFitter fitter(2,"pol1");
2827 fitter.StoreData(kFALSE);
2828 fitter.ClearPoints();
0c349049 2829 for(Int_t k = 0; k < npoints; k++){
d0569428 2830 //Take the cluster
2831 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2832 Short_t *signals = cl->GetSignals();
2833 Double_t time = cl->GetLocalTimeBin();
2834 //Calculate x if possible
2835 Float_t xcenter = 0.0;
2836 Bool_t echec = kTRUE;
2837 if((time<=7) || (time>=21)) continue;
2838 // Center 3 balanced: position with the center of the pad
2839 if ((((Float_t) signals[3]) > 0.0) &&
2840 (((Float_t) signals[2]) > 0.0) &&
2841 (((Float_t) signals[4]) > 0.0)) {
2842 echec = kFALSE;
2843 // Security if the denomiateur is 0
2844 if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2845 ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2846 xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2847 / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2848 / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2849 }
2850 else {
2851 echec = kTRUE;
2852 }
2853 }
2854 if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2855 if(echec) continue;
2856 //if no echec: calculate with the position of the pad
2857 // Position of the cluster
af26ce80 2858 Double_t padPosition = xcenter + cl->GetPadCol();
d0569428 2859 padPositions[k] = padPosition;
0c349049 2860 nb3pc++;
d0569428 2861 fitter.AddPoint(&time, padPosition,1);
2862 }//clusters loop
2863
0c349049 2864 //printf("nb3pc %d, npoints %d\n",nb3pc,npoints);
2865 if(nb3pc < 3) return kFALSE;
d0569428 2866 fitter.Eval();
2867 TVectorD line(2);
2868 fitter.GetParameters(line);
2869 Float_t pointError = -1.0;
0c349049 2870 pointError = TMath::Sqrt(fitter.GetChisquare()/nb3pc);
d0569428 2871
2872 // Take the tgl and snp with the track t now
2873 Double_t tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
2874 Double_t snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
2875 Float_t dzdx = 0.0; // dzdx
2876 Float_t tnp = 0.0;
2877 if(TMath::Abs(snp) < 1.0){
2878 tnp = snp / (TMath::Sqrt(1-snp*snp));
2879 dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2880 }
2881
2882
2883 // Now fill the PRF
0c349049 2884 for(Int_t k = 0; k < npoints; k++){
d0569428 2885 //Take the cluster
2886 AliTRDcluster *cl = (AliTRDcluster *) t->GetCluster(k+index0);
2887 Short_t *signals = cl->GetSignals(); // signal
2888 Double_t time = cl->GetLocalTimeBin(); // time bin
2889 Float_t padPosTracklet = line[0]+line[1]*time; // reconstruct position from fit
af26ce80 2890 Float_t padPos = cl->GetPadCol(); // middle pad
d0569428 2891 Double_t dpad = padPosTracklet - padPos; // reconstruct position relative to middle pad from fit
2892 Float_t ycenter = 0.0; // relative center charge
2893 Float_t ymin = 0.0; // relative left charge
2894 Float_t ymax = 0.0; // relative right charge
2895
2896
2897
2898 //Requiere simply two pads clusters at least
2899 if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2900 ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2901 Float_t sum = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2902 if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2903 if(sum > 0.0) ymin = ((Float_t) signals[2])/ sum;
2904 if(sum > 0.0) ymax = ((Float_t) signals[4])/ sum;
2905 }
2906
2907 //calibration group
2908 Int_t *rowcol = CalculateRowCol(cl); // calcul col and row pad of the cluster
2909 Int_t grouplocal = CalculateCalibrationGroup(2,rowcol); // calcul the corresponding group
2910 Int_t caligroup = fCalibraMode->GetXbins(2)+ grouplocal; // calcul the corresponding group
2911 Float_t xcl = cl->GetY(); // y cluster
2912 Float_t qcl = cl->GetQ(); // charge cluster
2913 Int_t plane = GetPlane(detector); // plane
2914 Int_t chamber = GetChamber(detector); // chamber
2915 Double_t xdiff = dpad; // reconstructed position constant
2916 Double_t x = dpad; // reconstructed position moved
0c349049 2917 Float_t ep = pointError; // error of fit
d0569428 2918 Float_t signal1 = (Float_t)signals[1]; // signal at the border
2919 Float_t signal3 = (Float_t)signals[3]; // signal
2920 Float_t signal2 = (Float_t)signals[2]; // signal
2921 Float_t signal4 = (Float_t)signals[4]; // signal
2922 Float_t signal5 = (Float_t)signals[5]; // signal at the border
2923
2924
2925
2926 if(fDebugLevel > 0){
2927 if ( !fDebugStreamer ) {
2928 //debug stream
2929 TDirectory *backup = gDirectory;
2930 fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2931 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
2932 }
2933
2934 (* fDebugStreamer) << "PRF0"<<
2935 "caligroup="<<caligroup<<
2936 "detector="<<detector<<
2937 "plane="<<plane<<
2938 "chamber="<<chamber<<
0c349049 2939 "npoints="<<npoints<<
2940 "Np="<<nb3pc<<
2941 "ep="<<ep<<
d0569428 2942 "snp="<<snp<<
2943 "tnp="<<tnp<<
2944 "tgl="<<tgl<<
2945 "dzdx="<<dzdx<<
2946 "padPos="<<padPos<<
2947 "padPosition="<<padPositions[k]<<
2948 "padPosTracklet="<<padPosTracklet<<
2949 "x="<<x<<
2950 "ycenter="<<ycenter<<
2951 "ymin="<<ymin<<
2952 "ymax="<<ymax<<
2953 "xcl="<<xcl<<
2954 "qcl="<<qcl<<
2955 "signal1="<<signal1<<
2956 "signal2="<<signal2<<
2957 "signal3="<<signal3<<
2958 "signal4="<<signal4<<
2959 "signal5="<<signal5<<
2960 "time="<<time<<
2961 "\n";
2962 x = xdiff;
2963 Int_t type=0;
2964 Float_t y = ycenter;
2965 (* fDebugStreamer) << "PRFALL"<<
2966 "caligroup="<<caligroup<<
2967 "detector="<<detector<<
2968 "plane="<<plane<<
2969 "chamber="<<chamber<<
0c349049 2970 "npoints="<<npoints<<
2971 "Np="<<nb3pc<<
2972 "ep="<<ep<<
d0569428 2973 "type="<<type<<
2974 "snp="<<snp<<
2975 "tnp="<<tnp<<
2976 "tgl="<<tgl<<
2977 "dzdx="<<dzdx<<
2978 "padPos="<<padPos<<
2979 "padPosition="<<padPositions[k]<<
2980 "padPosTracklet="<<padPosTracklet<<
2981 "x="<<x<<
2982 "y="<<y<<
2983 "xcl="<<xcl<<
2984 "qcl="<<qcl<<
2985 "signal1="<<signal1<<
2986 "signal2="<<signal2<<
2987 "signal3="<<signal3<<
2988 "signal4="<<signal4<<
2989 "signal5="<<signal5<<
2990 "time="<<time<<
2991 "\n";
2992 x=-(xdiff+1);
2993 y = ymin;
2994 type=-1;
2995 (* fDebugStreamer) << "PRFALL"<<
2996 "caligroup="<<caligroup<<
2997 "detector="<<detector<<
2998 "plane="<<plane<<
2999 "chamber="<<chamber<<
0c349049 3000 "npoints="<<npoints<<
3001 "Np="<<nb3pc<<
3002 "ep="<<ep<<
d0569428 3003 "type="<<type<<
3004 "snp="<<snp<<
3005 "tnp="<<tnp<<
3006 "tgl="<<tgl<<
3007 "dzdx="<<dzdx<<
3008 "padPos="<<padPos<<
3009 "padPosition="<<padPositions[k]<<
3010 "padPosTracklet="<<padPosTracklet<<
3011 "x="<<x<<
3012 "y="<<y<<
3013 "xcl="<<xcl<<
3014 "qcl="<<qcl<<
3015 "signal1="<<signal1<<
3016 "signal2="<<signal2<<
3017 "signal3="<<signal3<<
3018 "signal4="<<signal4<<
3019 "signal5="<<signal5<<
3020 "time="<<time<<
3021 "\n";
3022 x=1-xdiff;
3023 y = ymax;
3024 type=1;
3025
3026 (* fDebugStreamer) << "PRFALL"<<
3027 "caligroup="<<caligroup<<
3028 "detector="<<detector<<
3029 "plane="<<plane<<
3030 "chamber="<<chamber<<
0c349049 3031 "npoints="<<npoints<<
3032 "Np="<<nb3pc<<
3033 "ep="<<ep<<
d0569428 3034 "type="<<type<<
3035 "snp="<<snp<<
3036 "tnp="<<tnp<<
3037 "tgl="<<tgl<<
3038 "dzdx="<<dzdx<<
3039 "padPos="<<padPos<<
3040 "padPosition="<<padPositions[k]<<
3041 "padPosTracklet="<<padPosTracklet<<
3042 "x="<<x<<
3043 "y="<<y<<
3044 "xcl="<<xcl<<
3045 "qcl="<<qcl<<
3046 "signal1="<<signal1<<
3047 "signal2="<<signal2<<
3048 "signal3="<<signal3<<
3049 "signal4="<<signal4<<
3050 "signal5="<<signal5<<
3051 "time="<<time<<
3052 "\n";
3053
3054 }
3055
3056 // some cuts
0c349049 3057 if(npoints < fNumberClusters) continue;
3058 if(nb3pc <= 5) continue;
d0569428 3059 if((time >= 21) || (time < 7)) continue;
3060 if(TMath::Abs(snp) >= 1.0) continue;
3061 if(qcl < 80) continue;
3062
3063 Bool_t echec = kFALSE;
3064 Double_t shift = 0.0;
3065 //Calculate the shift in x coresponding to this tnp
3066 if(fNgroupprf != 0.0){
3067 shift = -3.0*(fNgroupprf-1)-1.5;
3068 Double_t limithigh = -0.2*(fNgroupprf-1);
3069 if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
3070 else{
3071 while(tnp > limithigh){
3072 limithigh += 0.2;
3073 shift += 3.0;
3074 }
3075 }
3076 }
3077 if (fHisto2d && !echec) {
3078 if(TMath::Abs(dpad) < 1.5) {
3079 fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
3080 fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
3081 }
3082 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3083 fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
3084 fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
3085 }
3086 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3087 fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
3088 fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
3089 }
3090 }
3091 //Not equivalent anymore here!
3092 if (fVector2d && !echec) {
3093 if(TMath::Abs(dpad) < 1.5) {
3094 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
3095 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
3096 }
3097 if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
3098 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
3099 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
3100 }
3101 if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
3102 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
3103 fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
3104 }
3105 }
3106 }
3107 return kTRUE;
3108
3109}
3110//
3111//____________Some basic geometry function_____________________________________
3112//
3113//_____________________________________________________________________________
3114Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
3115{
3116 //
3117 // Reconstruct the plane number from the detector number
3118 //
3119
3120 return ((Int_t) (d % 6));
3121
3122}
3123
3124//_____________________________________________________________________________
3125Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
3126{
3127 //
3128 // Reconstruct the chamber number from the detector number
3129 //
3130 Int_t fgkNplan = 6;
3131
3132 return ((Int_t) (d % 30) / fgkNplan);
3133
3134}
3135
3136//_____________________________________________________________________________
3137Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3138{
3139 //
3140 // Reconstruct the sector number from the detector number
3141 //
3142 Int_t fg = 30;
3143
3144 return ((Int_t) (d / fg));
3145
3146}
3147//_____________________________________________________________________
3148TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3149{
3150 //
3151 // return pointer to fPH2d TProfile2D
3152 // create a new TProfile2D if it doesn't exist allready
3153 //
3154 if ( fPH2d )
3155 return fPH2d;
3156
3157 // Some parameters
3158 fTimeMax = nbtimebin;
3159 fSf = samplefrequency;
3160
3161 /*
3162 AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
3163 ,fCalibraMode->GetNz(1)
3164 ,fCalibraMode->GetNrphi(1)));
170c35f1 3165
3166 // Calcul the number of Xbins
0c349049 3167 Int_t ntotal1 = 0;
170c35f1 3168 fCalibraMode->ModePadCalibration(2,1);
3169 fCalibraMode->ModePadFragmentation(0,2,0,1);
3170 fCalibraMode->SetDetChamb2(1);
0c349049 3171 ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
170c35f1 3172 fCalibraMode->ModePadCalibration(0,1);
3173 fCalibraMode->ModePadFragmentation(0,0,0,1);
3174 fCalibraMode->SetDetChamb0(1);
0c349049 3175 ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
3176 AliInfo(Form("Total number of Xbins: %d",ntotal1));
170c35f1 3177
0c349049 3178 CreatePH2d(ntotal1);
170c35f1 3179 */
3180
3181 CreatePH2d(540);
3182
3183 return fPH2d;
3184}
3185//_____________________________________________________________________
3186TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3187{
3188 //
3189 // return pointer to TLinearFitter Calibration
3190 // if force is true create a new TLinearFitter if it doesn't exist allready
3191 //
170c35f1 3192
d0569428 3193 if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3194 return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3195 }
3196
3197 // if we are forced and TLinearFitter doesn't yet exist create it
170c35f1 3198
d0569428 3199 // new TLinearFitter
3200 TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3201 fLinearFitterArray.AddAt(linearfitter,detector);
3202 return linearfitter;
170c35f1 3203}
170c35f1 3204
170c35f1 3205//_____________________________________________________________________
3206void AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3207{
3208 //
3209 // FillCH2d: Marian style
3210 //
3211
3212 //skip simply the value out of range
3213 if((y>=300.0) || (y<0.0)) return;
3214
3215 //Calcul the y place
3216 Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3217 Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3218
3219 //Fill
3220 fCH2d->GetArray()[place]++;
3221
3222}
3a0f6479 3223
3224//____________________________________________________________________________
3225void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3226{
3227 //
3228 // Analyse array of linear fitter because can not be written
3229 // Store two arrays: one with the param the other one with the error param + number of entries
3230 //
3231
3232 for(Int_t k = 0; k < 540; k++){
3233 TLinearFitter *linearfitter = GetLinearFitter(k);
3234 if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3235 TVectorD *par = new TVectorD(2);
3236 TVectorD pare = TVectorD(2);
3237 TVectorD *parE = new TVectorD(3);
3238 linearfitter->Eval();
3239 linearfitter->GetParameters(*par);
3240 linearfitter->GetErrors(pare);
3241 Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3242 (*parE)[0] = pare[0]*ppointError;
3243 (*parE)[1] = pare[1]*ppointError;
3244 (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3245 ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3246 ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3247 }
3248 }
3249}
d0569428 3250//________________________________________________________________________________
3251Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const
3252{
3253 // Partie entiere of the (x+0.5)
3254
3255 int i;
3256 if (x >= (-0.5)) {
3257 i = int(x + 0.5);
3258 //if (x + 0.5 == Float_t(i) && i & 1) i--;
3259 } else {
3260 i = int(x - 0.5);
3261 //if (x - 0.5 == Float_t(i) && i & 1) i++;
3262 if((x-0.5)==i) i++;
3263
3264 }
3265 return i;
3266}