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