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