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