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