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