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