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