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