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