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