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