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