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