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