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