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