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