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