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