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