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