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