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