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