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