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