Replace AliTRDCalibra
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraFillHisto.cxx
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
26 // calibration (see AliTRDCalibraMode). 
27 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
28 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking 
29 // in the function "FollowBackProlongation" (AliTRDtracker)
30 // Per default the functions to fill are off.                                   
31 //                        
32 // Author:
33 //   R. Bailhache (R.Bailhache@gsi.de)
34 //                            
35 //////////////////////////////////////////////////////////////////////////////////////
36
37 #include <TTree.h>
38 #include <TProfile2D.h>
39 #include <TProfile.h>
40 #include <TFile.h>
41 #include <TChain.h>
42 #include <TStyle.h>
43 #include <TCanvas.h>
44 #include <TGraphErrors.h>
45 #include <TObjArray.h>
46 #include <TH1F.h>
47 #include <TH2I.h>
48 #include <TH2.h>
49 #include <TStopwatch.h>
50 #include <TMath.h>
51 #include <TDirectory.h>
52 #include <TROOT.h>
53
54 #include "AliLog.h"
55 #include "AliCDBManager.h"
56
57 #include "AliTRDCalibraFillHisto.h"
58 #include "AliTRDCalibraMode.h"
59 #include "AliTRDCalibraVector.h"
60 #include "AliTRDcalibDB.h"
61 #include "AliTRDCommonParam.h"
62 #include "AliTRDmcmTracklet.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcluster.h"
65 #include "AliTRDtrack.h"
66
67
68 ClassImp(AliTRDCalibraFillHisto)
69
70 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
71 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
72
73 //_____________singleton implementation_________________________________________________
74 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
75 {
76   //
77   // Singleton implementation
78   //
79
80   if (fgTerminated != kFALSE) {
81     return 0;
82   }
83
84   if (fgInstance == 0) {
85     fgInstance = new AliTRDCalibraFillHisto();
86   }
87
88   return fgInstance;
89
90 }
91
92 //______________________________________________________________________________________
93 void AliTRDCalibraFillHisto::Terminate()
94 {
95   //
96   // Singleton implementation
97   // Deletes the instance of this class
98   //
99
100   fgTerminated = kTRUE;
101
102   if (fgInstance != 0) {
103     delete fgInstance;
104     fgInstance = 0;
105   }
106
107 }
108
109 //______________________________________________________________________________________
110 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
111   :TObject()
112   ,fMITracking(kFALSE)
113   ,fMcmTracking(kFALSE)
114   ,fMcmCorrectAngle(kFALSE)
115   ,fCH2dOn(kFALSE)
116   ,fPH2dOn(kFALSE)
117   ,fPRF2dOn(kFALSE)
118   ,fHisto2d(kFALSE)
119   ,fVector2d(kFALSE)
120   ,fRelativeScale(0)
121   ,fCountRelativeScale(0)
122   ,fRelativeScaleAuto(kFALSE)
123   ,fThresholdClusterPRF1(0.0)
124   ,fThresholdClusterPRF2(0.0)
125   ,fCenterOfflineCluster(kFALSE)
126   ,fWriteName(0)
127   ,fDetectorAliTRDtrack(kFALSE)
128   ,fChamberAliTRDtrack(-1)
129   ,fDetectorPreviousTrack(-1)
130   ,fGoodTrack(kTRUE)
131   ,fAmpTotal(0x0)
132   ,fPHPlace(0x0)
133   ,fPHValue(0x0)
134   ,fNumberClusters(0)
135   ,fProcent(0.0)
136   ,fDifference(0)
137   ,fNumberTrack(0)
138   ,fTimeMax(0)
139   ,fSf(0.0)
140   ,fNumberBinCharge(0)
141   ,fNumberBinPRF(0)
142   ,fCalibraVector(0)
143   ,fPH2d(0x0)
144   ,fPRF2d(0x0)
145   ,fCH2d(0x0)
146 {
147   //
148   // Default constructor
149   //
150
151   fCalibraMode = new AliTRDCalibraMode();
152
153   // Write
154   for (Int_t i = 0; i < 3; i++) {
155     fWrite[i]     = kFALSE;
156   }
157
158   // Init
159   Init();
160   
161 }
162
163 //______________________________________________________________________________________
164 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
165   :TObject(c)
166   ,fMITracking(kFALSE)
167   ,fMcmTracking(kFALSE)
168   ,fMcmCorrectAngle(kFALSE)
169   ,fCH2dOn(kFALSE)
170   ,fPH2dOn(kFALSE)
171   ,fPRF2dOn(kFALSE)
172   ,fHisto2d(kFALSE)
173   ,fVector2d(kFALSE)
174   ,fRelativeScale(0)
175   ,fCountRelativeScale(0)
176   ,fRelativeScaleAuto(kFALSE)
177   ,fThresholdClusterPRF1(0.0)
178   ,fThresholdClusterPRF2(0.0)
179   ,fCenterOfflineCluster(kFALSE)
180   ,fWriteName(0)
181   ,fDetectorAliTRDtrack(kFALSE)
182   ,fChamberAliTRDtrack(-1)
183   ,fDetectorPreviousTrack(-1)
184   ,fGoodTrack(kTRUE)
185   ,fAmpTotal(0x0)
186   ,fPHPlace(0x0)
187   ,fPHValue(0x0)
188   ,fNumberClusters(0)
189   ,fProcent(0.0)
190   ,fDifference(0)
191   ,fNumberTrack(0)
192   ,fTimeMax(0)
193   ,fSf(0.0)
194   ,fNumberBinCharge(0)
195   ,fNumberBinPRF(0)
196   ,fCalibraVector(0)
197   ,fPH2d(0x0)
198   ,fPRF2d(0x0)
199   ,fCH2d(0x0) 
200 {
201   //
202   // Copy constructor
203   //
204
205 }
206
207 //____________________________________________________________________________________
208 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
209 {
210   //
211   // AliTRDCalibraFillHisto destructor
212   //
213
214   ClearHistos();
215   
216 }
217
218 //_____________________________________________________________________________
219 void AliTRDCalibraFillHisto::Destroy() 
220 {
221   //
222   // Delete instance 
223   //
224
225   if (fgInstance) {
226     delete fgInstance;
227     fgInstance = 0x0;
228   }
229
230 }
231
232 //_____________________________________________________________________________
233 void AliTRDCalibraFillHisto::ClearHistos() 
234 {
235   //
236   // Delete the histos
237   //
238
239   if (fPH2d) {
240     delete fPH2d;
241     fPH2d  = 0x0;
242   }
243   if (fCH2d) {
244     delete fCH2d;
245     fCH2d  = 0x0;
246   }
247   if (fPRF2d) {
248     delete fPRF2d;
249     fPRF2d = 0x0;
250   }
251
252 }
253
254 //_____________________________________________________________________________
255 void AliTRDCalibraFillHisto::Init() 
256 {
257   //
258   // Init some default values
259   //
260
261   // How to fill the 2D
262   fThresholdClusterPRF1 = 2.0;
263   fThresholdClusterPRF2 = 15.0;
264   
265   // Store the Info
266   fNumberBinCharge      = 100;
267   fNumberBinPRF         = 40;
268   
269   // Write
270   fWriteName            = "TRD.calibration.root";
271   
272   // Internal variables
273   
274   // Fill the 2D histos in the offline tracking
275   fDetectorPreviousTrack = -1;
276   fChamberAliTRDtrack    = -1;
277   fGoodTrack             = kTRUE;
278
279   fProcent               = 6.0;
280   fDifference            = 17;
281   fNumberClusters        = 18;
282   fNumberTrack           = 0;
283   fNumberUsedCh[0]       = 0;
284   fNumberUsedCh[1]       = 0;
285   fNumberUsedPh[0]       = 0;
286   fNumberUsedPh[1]       = 0;
287  
288 }
289
290 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
291 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
292 {
293   //
294   // For the offline tracking
295   // This function will be called in the function AliReconstruction::Run() 
296   // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE, 
297   //
298
299   // DB Setting
300   // Get cal
301   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
302   if (!cal) {
303     AliInfo("Could not get calibDB");
304     return kFALSE;
305   }
306   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
307   if (!parCom) {
308     AliInfo("Could not get CommonParam");
309     return kFALSE;
310   }
311
312   // Some parameters
313   fTimeMax = cal->GetNumberOfTimeBins();
314   fSf      = parCom->GetSamplingFrequency();
315   if (fRelativeScaleAuto) {
316     fRelativeScale = 0;
317   }
318   else {
319     fRelativeScale = 20;
320   }
321
322   //If vector method On initialised all the stuff
323   if(fVector2d){
324     fCalibraVector = new AliTRDCalibraVector();
325   }
326
327
328   // Create the 2D histos corresponding to the pad groupCalibration mode
329   if (fCH2dOn) {
330
331     AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
332                 ,fCalibraMode->GetNz(0)
333                 ,fCalibraMode->GetNrphi(0)));
334     
335     // Calcul the number of Xbins
336     Int_t Ntotal0 = 0;
337     fCalibraMode->ModePadCalibration(2,0);
338     fCalibraMode->ModePadFragmentation(0,2,0,0);
339     fCalibraMode->SetDetChamb2(0);
340     Ntotal0 += 6 * 18 * fCalibraMode->GetDetChamb2(0);
341     fCalibraMode->ModePadCalibration(0,0);
342     fCalibraMode->ModePadFragmentation(0,0,0,0);
343     fCalibraMode->SetDetChamb0(0);
344     Ntotal0 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(0);
345     AliInfo(Form("Total number of Xbins: %d",Ntotal0));
346
347     // Create the 2D histo
348     if (fHisto2d) {
349       CreateCH2d(Ntotal0);
350     }
351     if (fVector2d) {
352       fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
353     }
354
355     // Variable
356     fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
357     for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
358       fAmpTotal[k] = 0.0;
359     } 
360
361   }
362
363   if (fPH2dOn) {
364
365     AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
366                 ,fCalibraMode->GetNz(1)
367                 ,fCalibraMode->GetNrphi(1)));
368     
369     // Calcul the number of Xbins
370     Int_t Ntotal1 = 0;
371     fCalibraMode->ModePadCalibration(2,1);
372     fCalibraMode->ModePadFragmentation(0,2,0,1);
373     fCalibraMode->SetDetChamb2(1);
374     Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
375     fCalibraMode->ModePadCalibration(0,1);
376     fCalibraMode->ModePadFragmentation(0,0,0,1);
377     fCalibraMode->SetDetChamb0(1);
378     Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
379     AliInfo(Form("Total number of Xbins: %d",Ntotal1));
380
381     // Create the 2D histo
382     if (fHisto2d) {
383       CreatePH2d(Ntotal1);
384     }
385     if (fVector2d) {
386       fCalibraVector->SetTimeMax(fTimeMax);
387     }
388    
389     // Variable
390     fPHPlace = new Short_t[fTimeMax];
391     for (Int_t k = 0; k < fTimeMax; k++) {
392       fPHPlace[k] = -1;
393     } 
394     fPHValue = new Float_t[fTimeMax];
395     for (Int_t k = 0; k < fTimeMax; k++) {
396       fPHValue[k] = 0.0;
397     }
398
399   }
400
401   if (fPRF2dOn) {
402
403     AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
404                 ,fCalibraMode->GetNz(2)
405                 ,fCalibraMode->GetNrphi(2)));
406     
407     // Calcul the number of Xbins
408     Int_t Ntotal2 = 0;
409     fCalibraMode->ModePadCalibration(2,2);
410     fCalibraMode->ModePadFragmentation(0,2,0,2);
411     fCalibraMode->SetDetChamb2(2);
412     Ntotal2 += 6 * 18 * fCalibraMode->GetDetChamb2(2);
413     fCalibraMode->ModePadCalibration(0,2);
414     fCalibraMode->ModePadFragmentation(0,0,0,2);
415     fCalibraMode->SetDetChamb0(2);
416     Ntotal2 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(2);
417     AliInfo(Form("Total number of Xbins: %d",Ntotal2));
418
419     // Create the 2D histo
420     if (fHisto2d) {
421       CreatePRF2d(Ntotal2);
422     }
423     if (fVector2d) {
424       fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
425     }
426   
427   }
428
429   return kTRUE;
430
431 }
432
433 //____________Functions for filling the histos in the code_____________________
434
435 //____________Offine tracking in the AliTRDtracker_____________________________
436 Bool_t AliTRDCalibraFillHisto::ResetTrack()
437 {
438   //
439   // For the offline tracking
440   // This function will be called in the function
441   // AliTRDtracker::FollowBackPropagation() at the beginning 
442   // Reset the parameter to know we have a new TRD track
443   //
444   
445   fDetectorAliTRDtrack = kFALSE;
446   return kTRUE;
447
448 }
449
450 //____________Offline tracking in the AliTRDtracker____________________________
451 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
452 {
453   //
454   // For the offline tracking
455   // This function will be called in the function
456   // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
457   // of TRD tracks 
458   // Fill the 2D histos or the vectors with the info of the clusters at
459   // the end of a detectors if the track is "good"
460   //
461
462   // Get the parameter object
463   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
464   if (!parCom) {
465     AliInfo("Could not get CommonParam");
466     return kFALSE;
467   }
468
469   // Get the parameter object
470   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
471   if (!cal) {
472     AliInfo("Could not get calibDB");
473     return kFALSE;
474   }
475  
476   // Localisation of the detector
477   Int_t detector = cl->GetDetector();
478   Int_t chamber  = GetChamber(detector);
479   Int_t plane    = GetPlane(detector);
480
481   // Fill the infos for the previous clusters if not the same
482   // detector anymore or if not the same track
483   if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) && 
484       (fDetectorPreviousTrack != -1)) {
485     
486     fNumberTrack++;   
487     
488     // If the same track, then look if the previous detector is in
489     // the same plane, if yes: not a good track
490     if (fDetectorAliTRDtrack && 
491         (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
492       fGoodTrack = kFALSE;
493     }
494
495     // Fill only if the track doesn't touch a masked pad or doesn't
496     // appear in the middle (fGoodTrack)
497     if (fGoodTrack) {
498
499       // Gain calibration
500       if (fCH2dOn) {
501         FillTheInfoOfTheTrackCH();
502       }
503       
504       // PH calibration
505       if (fPH2dOn) {
506         FillTheInfoOfTheTrackPH();    
507       }
508       
509     } // if a good track
510     
511     ResetfVariables();
512     
513   } // Fill at the end the charge
514   
515   // Calcul the position of the detector
516   if (detector != fDetectorPreviousTrack) {
517     LocalisationDetectorXbins(detector);
518   }
519
520   // Reset the good track for the PRF
521   Bool_t good = kTRUE;
522   
523   // Localisation of the cluster
524   Double_t pos[3] = { 0.0, 0.0, 0.0 };
525   pos[0] = cl->GetX();
526   pos[1] = cl->GetY();
527   pos[2] = cl->GetZ();
528   Int_t    time   = cl->GetLocalTimeBin();
529   
530   // Reset the detector
531   fDetectorPreviousTrack = detector;
532   fDetectorAliTRDtrack   = kTRUE;
533   
534   // Position of the cluster
535   AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
536   Int_t    row        = padplane->GetPadRowNumber(pos[2]);
537   Double_t offsetz    = padplane->GetPadRowOffset(row,pos[2]);
538   Double_t offsettilt = padplane->GetTiltOffset(offsetz);
539   Int_t    col        = padplane->GetPadColNumber(pos[1]+offsettilt);
540   
541   // See if we are not near a masked pad
542   if (!IsPadOn(detector,col,row)) {
543     good       = kFALSE;
544     fGoodTrack = kFALSE;
545   }
546
547   if (col > 0) {
548     if (!IsPadOn(detector,col-1,row)) {
549       fGoodTrack = kFALSE;
550       good       = kFALSE;
551     }
552   }
553
554   if (col < 143) {
555     if (!IsPadOn(detector,col+1,row)) {
556       fGoodTrack = kFALSE;
557       good       = kFALSE;
558     }
559   }
560
561   // Row of the cluster and position in the pad groups
562   Int_t posr[3] = { 0, 0, 0 };
563   if ((fCH2dOn)  && (fCalibraMode->GetNnZ(0) != 0)) {
564     posr[0] = (Int_t) row / fCalibraMode->GetNnZ(0);
565   }
566   if ((fPH2dOn)  && (fCalibraMode->GetNnZ(1) != 0)) {
567     posr[1] = (Int_t) row / fCalibraMode->GetNnZ(1);
568   }
569   if ((fPRF2dOn) && (fCalibraMode->GetNnZ(2) != 0)) {
570     posr[2] = (Int_t) row / fCalibraMode->GetNnZ(2);
571   }  
572       
573   // Col of the cluster and position in the pad groups
574   Int_t posc[3] = { 0, 0, 0 };
575   if ((fCH2dOn)  && (fCalibraMode->GetNnRphi(0) != 0)) {
576     posc[0] = (Int_t) col / fCalibraMode->GetNnRphi(0);
577   }
578   if ((fPH2dOn)  && (fCalibraMode->GetNnRphi(1) != 0)) {
579     posc[1] = (Int_t) col / fCalibraMode->GetNnRphi(1);
580   }
581   if ((fPRF2dOn) && (fCalibraMode->GetNnRphi(2) != 0)) {
582     posc[2] = (Int_t) col / fCalibraMode->GetNnRphi(2);
583   }
584
585   // Charge in the cluster
586   // For the moment take the abs
587   Float_t  q        = TMath::Abs(cl->GetQ());
588   Short_t  *signals = cl->GetSignals();
589  
590   // Correction due to the track angle
591   Float_t correction    = 1.0;
592   Float_t normalisation = 6.67;
593   if ((q >0) && (t->GetNdedx() > 0)) {
594     correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
595   }
596
597   // Fill the fAmpTotal with the charge
598   if (fCH2dOn) {
599     fAmpTotal[(Int_t) (posc[0]*fCalibraMode->GetNfragZ(0)+posr[0])] += correction;
600   }
601
602   // Fill the fPHPlace and value
603   if (fPH2dOn) {
604     fPHPlace[time] = posc[1]*fCalibraMode->GetNfragZ(1)+posr[1];
605     fPHValue[time] = correction;
606   }
607
608   // Fill direct the PRF
609   if ((fPRF2dOn) && (good)) {
610
611     Float_t yminus  = 0.0;
612     Float_t xcenter = 0.0;
613     Float_t ycenter = 0.0;
614     Float_t ymax    = 0.0;
615     Bool_t  echec   = kFALSE;
616     
617     if ((cl->From3pad()) && (!cl->IsUsed())) { 
618          
619       // Center 3 balanced and cut on the cluster shape
620       if ((((Float_t) signals[3]) > fThresholdClusterPRF2) && 
621           (((Float_t) signals[2]) > fThresholdClusterPRF2) && 
622           (((Float_t) signals[4]) > fThresholdClusterPRF2) && 
623           (((Float_t) signals[1]) < fThresholdClusterPRF1) && 
624           (((Float_t) signals[5]) < fThresholdClusterPRF1) && 
625           ((((Float_t) signals[2])*((Float_t) signals[4])/(((Float_t) signals[3])*((Float_t) signals[3]))) < 0.06)) {
626
627
628         //First calculate the position of the cluster and the y
629         //echec enables to repair cases where it fails
630         //
631         // Col correspond to signals[3]
632         if (fCenterOfflineCluster) {
633           xcenter = cl->GetCenter();
634         }
635         else {
636           // Security if the denomiateur is 0 
637           if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / 
638                            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
639             xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
640                           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) 
641                                       / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
642           }
643           else {
644             echec = kTRUE;
645           }
646         }
647         //after having calculating the position calculate the y
648         if (TMath::Abs(xcenter) < 0.5) {
649           ycenter = (Float_t) (((Float_t) signals[3]) 
650                             / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
651           yminus  = (Float_t) (((Float_t) signals[2]) 
652                             / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
653           ymax    = (Float_t) (((Float_t) signals[4]) 
654                             / (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
655           //If the charge of the cluster is too far away from the corrected one cut
656           if ((TMath::Abs(((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4])) - q) > 10.0)) {
657             echec = kTRUE;
658           }
659         }
660         else {
661           echec = kTRUE;
662         }
663       
664         //Then Fill the histo if no echec
665         //
666         // Fill only if it is in the drift region!
667         if ((((Float_t) (((Float_t) time) / fSf)) > 0.3) && (!echec)) {
668           if (fHisto2d) {
669             fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),xcenter,ycenter);
670             fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),-(xcenter+1.0),yminus);
671             fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),1.0-xcenter,ymax);
672           }
673           if (fVector2d) {
674             fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],xcenter,ycenter);
675             fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],-(xcenter+1.0),yminus);
676             fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],1.0-xcenter,ymax);
677           }
678         } // If in the drift region
679       } // center 3 balanced and cut on the cluster shape
680     } // Cluster isole
681   } // PRF2dOn  
682   
683   return kTRUE;
684   
685 }
686
687 //____________Online trackling in AliTRDtrigger________________________________
688 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
689 {
690   //
691   // For the tracking
692   // This function will be called in the function AliTRDtrigger::TestTracklet
693   // before applying the pt cut on the tracklets 
694   // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
695   //
696   
697   // Localisation of the Xbins involved
698   Int_t idect = trk->GetDetector();
699   LocalisationDetectorXbins(idect);
700
701   // Get the parameter object
702   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
703   if (!cal) {
704     AliInfo("Could not get calibDB");
705     return kFALSE;
706   }
707    
708   // Reset
709   ResetfVariables();
710  
711   // Row of the tracklet and position in the pad groups
712   Int_t row     = trk->GetRow();
713   Int_t posr[3] = { 0, 0, 0 };
714   if ((fCH2dOn)  && (fCalibraMode->GetNnZ(0) != 0)) {
715     posr[0] = (Int_t) row / fCalibraMode->GetNnZ(0);
716   }
717   if ((fPH2dOn)  && (fCalibraMode->GetNnZ(1) != 0)) {
718     posr[1] = (Int_t) row / fCalibraMode->GetNnZ(1);
719   }
720   if ((fPRF2dOn) && (fCalibraMode->GetNnZ(2) != 0)) {
721     posr[2] = (Int_t) row / fCalibraMode->GetNnZ(2);
722   }
723  
724   // Eventuelle correction due to track angle in z direction
725   Float_t correction = 1.0;
726   if (fMcmCorrectAngle) {
727     Float_t z = trk->GetRowz();
728     Float_t r = trk->GetTime0();
729     correction = r / TMath::Sqrt((r*r+z*z));
730   }
731
732   // Boucle sur les clusters
733   // Condition on number of cluster: don't come from the middle of the detector
734   if (trk->GetNclusters() >= fNumberClusters) {
735
736     for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
737
738       Float_t amp[3] = { 0.0, 0.0, 0.0 };
739       Int_t   time   = trk->GetClusterTime(icl);
740       Int_t   col    = trk->GetClusterCol(icl);
741             
742       amp[0] = trk->GetClusterADC(icl)[0] * correction;
743       amp[1] = trk->GetClusterADC(icl)[1] * correction;
744       amp[2] = trk->GetClusterADC(icl)[2] * correction;
745
746       
747       if ((amp[0] < 0.0) || 
748           (amp[1] < 0.0) || 
749           (amp[2] < 0.0)) {
750         continue;
751       }
752
753       // Col of cluster and position in the pad groups
754       Int_t posc[3] = { 0, 0, 0 };
755       if ((fCH2dOn)  && (fCalibraMode->GetNnRphi(0) != 0)) {
756         posc[0] = (Int_t) col / fCalibraMode->GetNnRphi(0);
757       }
758       if ((fPH2dOn)  && (fCalibraMode->GetNnRphi(1) != 0)) {
759         posc[1] = (Int_t) col / fCalibraMode->GetNnRphi(1);
760       }
761       if ((fPRF2dOn) && (fCalibraMode->GetNnRphi(2) != 0)) {
762         posc[2] = (Int_t) col / fCalibraMode->GetNnRphi(2);
763       }
764
765       // See if we are not near a masked pad
766       Bool_t good = kTRUE;
767       if (!IsPadOn(idect,col,row)) {
768         fGoodTrack = kFALSE;
769         good       = kFALSE;
770       }
771
772       if (col >   0) {
773         if (!IsPadOn(idect,col-1,row)) {
774           fGoodTrack = kFALSE;
775           good       = kFALSE;
776         }
777       }
778       
779       if (col < 143) {
780         if (!IsPadOn(idect,col+1,row)) {
781           fGoodTrack = kFALSE;
782           good       = kFALSE;
783         }
784       }
785
786       // Total spectrum
787       if (fPH2dOn) {
788         fPHPlace[time] = posc[1] * fCalibraMode->GetNfragZ(1) + posr[1];
789       }
790
791       if (fCH2dOn) {
792         fAmpTotal[(Int_t) (posc[0]*fCalibraMode->GetNfragZ(0)+posr[0])] += (Float_t) (amp[0]+amp[1]+amp[2]);
793       }
794       if (fPH2dOn) {
795         fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
796       }
797       
798             
799       // Fill PRF direct
800       if (fPRF2dOn && good) {
801         
802         if ((amp[0] > fThresholdClusterPRF2) && 
803             (amp[1] > fThresholdClusterPRF2) && 
804             (amp[2] > fThresholdClusterPRF2) && 
805             ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
806         
807           // Security of the denomiateur is 0
808           if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1]))) 
809              / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
810             Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
811                                   / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
812             Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
813
814             if (TMath::Abs(xcenter) < 0.5) {
815               Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
816               Float_t ymax   = amp[2] / (amp[0]+amp[1]+amp[2]);
817               // Fill only if it is in the drift region!
818               if (((Float_t) time / fSf) > 0.3) {
819                 if (fHisto2d) {
820                   fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),xcenter,ycenter);
821                   fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),-(xcenter+1.0),yminus);
822                   fPRF2d->Fill((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]+0.5),(1.0-xcenter),ymax);
823                 }
824                 if (fVector2d) {
825                   fCalibraVector->UpdateVectorPRF((fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2]),xcenter,ycenter);
826                   fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],-(xcenter+1.0),yminus);
827                   fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+posc[2]*fCalibraMode->GetNfragZ(2)+posr[2],(1.0-xcenter),ymax);
828                 }
829               }//in the drift region 
830             }//in the middle
831           }//denominateur security
832         }//cluster shape and thresholds
833       }//good and PRF On
834       
835     } // Boucle clusters
836     
837     // Fill the charge
838     if (fCH2dOn && fGoodTrack) {
839       FillTheInfoOfTheTrackCH();
840     }
841
842     // PH calibration
843     if (fPH2dOn && fGoodTrack) {
844       FillTheInfoOfTheTrackPH();        
845     }
846
847     fNumberTrack++;
848         
849   } // Condition on number of clusters
850
851   return kTRUE;
852   
853 }
854
855 //_____________________________________________________________________________
856 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
857 {
858   //
859   // Look in the choosen database if the pad is On.
860   // If no the track will be "not good"
861   //
862
863   // Get the parameter object
864   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
865   if (!cal) {
866     AliInfo("Could not get calibDB");
867     return kFALSE;
868   }
869   
870   if (!cal->IsChamberInstalled(detector)     || 
871        cal->IsChamberMasked(detector)        ||
872        cal->IsPadMasked(detector,col,row)) {
873     return kFALSE;
874   }
875   else {
876     return kTRUE;
877   }
878   
879 }
880
881 //____________Functions for plotting the 2D____________________________________
882
883 //_____________________________________________________________________________
884 void AliTRDCalibraFillHisto::Plot2d()
885 {
886   //
887   // Plot the 2D histos 
888   //
889  
890   if (fPH2dOn) {
891     TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
892     cph2d->cd();
893     fPH2d->Draw("LEGO");
894   }
895   if (fCH2dOn) {
896     TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
897     cch2d->cd();
898     fCH2d->Draw("LEGO");
899   }
900   if (fPRF2dOn) {
901     TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
902     cPRF2d->cd();
903     fPRF2d->Draw("LEGO");
904   }
905
906 }
907
908 //____________Writing the 2D___________________________________________________
909
910 //_____________________________________________________________________________
911 Bool_t AliTRDCalibraFillHisto::Write2d()
912 {
913   //
914   // Write the 2D histograms or the vectors converted in trees in the file
915   // "TRD.calibration.root" 
916   //
917   
918   TFile *fout = TFile::Open(fWriteName,"RECREATE");
919   // Check if the file could be opened
920   if (!fout || !fout->IsOpen()) {
921     AliInfo("No File found!");
922     return kFALSE;
923   }
924   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
925               ,fNumberTrack
926               ,fNumberUsedCh[0]
927               ,fNumberUsedCh[1]
928               ,fNumberUsedPh[0]
929               ,fNumberUsedPh[1]));
930   
931   TStopwatch stopwatch;
932   stopwatch.Start();
933   AliInfo("Write2d");
934
935   if ((fCH2dOn ) && (fWrite[0])) {
936     if (fHisto2d) {
937       fout->WriteTObject(fCH2d);
938     }
939     if (fVector2d) {
940       TString name("Nz");
941       name += fCalibraMode->GetNz(0);
942       name += "Nrphi";
943       name += fCalibraMode->GetNrphi(0);
944       TTree *treeCH2d = fCalibraVector->ConvertVectorCTTreeHisto(fCalibraVector->GetVectorCH(),fCalibraVector->GetPlaCH(),"treeCH2d",(const char *) name);
945       fout->WriteTObject(treeCH2d);
946     }
947   }
948   if ((fPH2dOn ) && (fWrite[1])) {
949     if (fHisto2d) {
950       fout->WriteTObject(fPH2d);
951     }
952     if (fVector2d) {
953       TString name("Nz");
954       name += fCalibraMode->GetNz(1);
955       name += "Nrphi";
956       name += fCalibraMode->GetNrphi(1);
957       TTree *treePH2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPH(),fCalibraVector->GetPlaPH(),"treePH2d",(const char *) name);
958       fout->WriteTObject(treePH2d);
959     }
960   }
961   if ((fPRF2dOn ) && (fWrite[2])) {
962     if (fHisto2d) {
963       fout->WriteTObject(fPRF2d);
964     }
965     if (fVector2d) {
966       TString name("Nz");
967       name += fCalibraMode->GetNz(2);
968       name += "Nrphi";
969       name += fCalibraMode->GetNrphi(2);
970       TTree *treePRF2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPRF(),fCalibraVector->GetPlaPRF(),"treePRF2d",(const char *) name);
971       fout->WriteTObject(treePRF2d);
972     }
973   }
974   
975   fout->Close();
976   
977   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
978               ,stopwatch.RealTime(),stopwatch.CpuTime()));
979
980   return kTRUE;
981   
982 }
983
984 //____________Probe the histos_________________________________________________
985 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
986 {
987   //
988   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
989   // debug mode with 2 for TH2I and 3 for TProfile2D
990   // It gives a pointer to a Double_t[7] with the info following...
991   // [0] : number of calibration groups with entries
992   // [1] : minimal number of entries found
993   // [2] : calibration group number of the min
994   // [3] : maximal number of entries found
995   // [4] : calibration group number of the max
996   // [5] : mean number of entries found
997   // [6] : mean relativ error
998   //
999
1000   Double_t *info = new Double_t[7];
1001    
1002   // Number of Xbins (detectors or groups of pads)
1003   Int_t    nbins   = h->GetNbinsX(); //number of calibration groups
1004   Int_t    nybins  = h->GetNbinsY(); //number of bins per histo
1005
1006   // Initialise
1007   Double_t nbwe = 0; //number of calibration groups with entries
1008   Double_t minentries = 0; //minimal number of entries found
1009   Double_t maxentries = 0; //maximal number of entries found
1010   Double_t placemin = 0; //calibration group number of the min
1011   Double_t placemax = -1; //calibration group number of the max
1012   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1013   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1014
1015   Double_t counter = 0;
1016
1017   //Debug
1018   TH1F *NbEntries = 0x0;//distribution of the number of entries
1019   TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1020   TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1021     
1022   // Beginning of the loop over the calibration groups 
1023   for (Int_t idect = 0; idect < nbins; idect++) {
1024
1025     TH1I *projch = (TH1I *) h->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
1026     projch->SetDirectory(0);
1027     
1028     // Number of entries for this calibration group
1029     Double_t nentries = 0.0;
1030     if((i%2) == 0){
1031       for (Int_t k = 0; k < nybins; k++) {
1032         nentries += h->GetBinContent(h->GetBin(idect+1,k+1));
1033       }
1034     }
1035     else{
1036       for (Int_t k = 0; k < nybins; k++) {
1037         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(idect+1,k+1));
1038         if(h->GetBinContent(h->GetBin(idect+1,k+1)) != 0) {
1039           meanrelativerror += (h->GetBinError(h->GetBin(idect+1,k+1))
1040                             / (TMath::Abs(h->GetBinContent(h->GetBin(idect+1,k+1)))));
1041           counter++;
1042         } 
1043       }
1044     }
1045
1046     //Debug
1047     if(i > 1){
1048       if((!((Bool_t)NbEntries)) && (nentries > 0)){
1049         NbEntries = new TH1F("Number of entries","Number of entries"
1050                                ,100,(Int_t)nentries/2,nentries*2);
1051         NbEntries->SetDirectory(0);
1052         NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1053                                ,nbins,0,nbins);
1054         NbEntriesPerGroup->SetDirectory(0);
1055         NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1056                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1057         NbEntriesPerSp->SetDirectory(0);
1058       }
1059       if(NbEntries){
1060         if(nentries > 0) NbEntries->Fill(nentries);
1061         NbEntriesPerGroup->Fill(idect+0.5,nentries);
1062         NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1063       }
1064     }
1065
1066     //min amd max
1067     if(nentries > maxentries){
1068       maxentries = nentries;
1069       placemax = idect;
1070     }
1071     if(idect == 0) {
1072       minentries = nentries;
1073     }
1074     if(nentries < minentries){
1075       minentries = nentries;
1076       placemin = idect;
1077     }
1078     //nbwe
1079     if(nentries > 0) {
1080       nbwe++;
1081       meanstats += nentries;
1082     }
1083
1084   }//calibration groups loop
1085   
1086   if(nbwe > 0) meanstats /= nbwe;
1087   if(counter > 0) meanrelativerror /= counter;
1088
1089   AliInfo(Form("There are %f calibration groups with entries",nbwe));
1090   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1091   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1092   AliInfo(Form("The mean number of entries is %f",meanstats));
1093   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1094
1095   info[0] = nbwe;
1096   info[1] = minentries;
1097   info[2] = placemin;
1098   info[3] = maxentries;
1099   info[4] = placemax;
1100   info[5] = meanstats;
1101   info[6] = meanrelativerror;
1102
1103   if(i > 1){
1104     gStyle->SetPalette(1);
1105     gStyle->SetOptStat(1111);
1106     gStyle->SetPadBorderMode(0);
1107     gStyle->SetCanvasColor(10);
1108     gStyle->SetPadLeftMargin(0.13);
1109     gStyle->SetPadRightMargin(0.01);
1110     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1111     stat->Divide(2,1);
1112     stat->cd(1);
1113     NbEntries->Draw("");
1114     stat->cd(2);
1115     NbEntriesPerSp->SetStats(0);
1116     NbEntriesPerSp->Draw("");
1117     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1118     stat1->cd();
1119     NbEntriesPerGroup->SetStats(0);
1120     NbEntriesPerGroup->Draw("");
1121   }
1122
1123   return info;
1124
1125 }
1126
1127 //_____________________________________________________________________________
1128 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1129 {
1130   //
1131   // Set the factor that will divide the deposited charge
1132   // to fit in the histo range [0,300]
1133   //
1134  
1135   if (RelativeScale > 0.0) {
1136     fRelativeScale = RelativeScale;
1137   } 
1138   else {
1139     AliInfo("RelativeScale must be strict positif!");
1140   }
1141
1142
1143
1144 //_____________________________________________________________________________
1145 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1146 {
1147   //
1148   // Set the mode of calibration group in the z direction for the parameter i
1149   // 
1150
1151   if ((Nz >= 0) && 
1152       (Nz <  5)) {
1153     fCalibraMode->SetNz(i, Nz); 
1154   }
1155   else { 
1156     AliInfo("You have to choose between 0 and 4");
1157   }
1158
1159 }
1160
1161 //_____________________________________________________________________________
1162 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1163 {
1164   //
1165   // Set the mode of calibration group in the rphi direction for the parameter i
1166   //
1167  
1168   if ((Nrphi >= 0) && 
1169       (Nrphi <  7)) {
1170     fCalibraMode->SetNrphi(i ,Nrphi); 
1171   }
1172   else {
1173     AliInfo("You have to choose between 0 and 6");
1174   }
1175
1176 }
1177
1178 //____________Protected Functions______________________________________________
1179 //____________Create the 2D histo to be filled online__________________________
1180 //
1181
1182 //_____________________________________________________________________________
1183 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1184 {
1185   //
1186   // Create the 2D histos
1187   //
1188
1189   TString name("Nz");
1190   name += fCalibraMode->GetNz(2);
1191   name += "Nrphi";
1192   name += fCalibraMode->GetNrphi(2);
1193
1194   fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1195                                  ,nn,0,nn,fNumberBinPRF,-1.5,1.5);
1196   fPRF2d->SetXTitle("Det/pad groups");
1197   fPRF2d->SetYTitle("Position x/W [pad width units]");
1198   fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1199   fPRF2d->SetStats(0);
1200
1201 }
1202
1203 //_____________________________________________________________________________
1204 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1205 {
1206   //
1207   // Create the 2D histos
1208   //
1209
1210   TString name("Nz");
1211   name += fCalibraMode->GetNz(1);
1212   name += "Nrphi";
1213   name += fCalibraMode->GetNrphi(1);
1214
1215   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1216                                ,nn,0,nn,fTimeMax
1217                                ,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf);
1218   fPH2d->SetXTitle("Det/pad groups");
1219   fPH2d->SetYTitle("time [#mus]");
1220   fPH2d->SetZTitle("<PH> [a.u.]");
1221   fPH2d->SetStats(0);
1222
1223 }
1224
1225 //_____________________________________________________________________________
1226 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1227 {
1228   //
1229   // Create the 2D histos
1230   //
1231
1232   TString name("Nz");
1233   name += fCalibraMode->GetNz(0);
1234   name += "Nrphi";
1235   name += fCalibraMode->GetNrphi(0);
1236
1237   fCH2d = new TH2I("CH2d",(const Char_t *) name
1238                          ,nn,0,nn,fNumberBinCharge,0,300);
1239   fCH2d->SetXTitle("Det/pad groups");
1240   fCH2d->SetYTitle("charge deposit [a.u]");
1241   fCH2d->SetZTitle("counts");
1242   fCH2d->SetStats(0);
1243   fCH2d->Sumw2();
1244
1245 }
1246
1247 //____________Offine tracking in the AliTRDtracker_____________________________
1248 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1249 {
1250   //
1251   // For the offline tracking or mcm tracklets
1252   // This function will be called in the functions UpdateHistogram... 
1253   // to fill the info of a track for the relativ gain calibration
1254   //
1255         
1256   Int_t nb =  0; // Nombre de zones traversees
1257   Int_t fd = -1; // Premiere zone non nulle
1258   
1259   
1260   // See if the track goes through different zones
1261   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1262     if (fAmpTotal[k] > 0.0) {
1263       nb++;
1264       if (nb == 1) {
1265         fd = k;
1266       }
1267     }
1268   }
1269  
1270   // If automatic scale
1271   if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
1272     // Take only the one zone track
1273     if (nb == 1) {
1274       fRelativeScale += fAmpTotal[fd] * 0.014 * 0.01;
1275       fCountRelativeScale++;
1276     }
1277   }
1278
1279   // We fill the CH2d after having scale with the first 100
1280   if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
1281     // Case of track with only one zone
1282     if (nb == 1) {
1283       if (fHisto2d) {
1284         fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1285       }
1286       if (fVector2d) {
1287         fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1288       }
1289     } // Case 1 zone
1290     // Case of track with two zones
1291     if (nb == 2) {
1292       // Two zones voisines sinon rien!
1293       if ((fAmpTotal[fd]   > 0.0) && 
1294           (fAmpTotal[fd+1] > 0.0)) {
1295         // One of the two very big
1296         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1297           if (fHisto2d) {
1298             fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1299           }
1300           if (fVector2d) {
1301             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1302           }
1303         }
1304         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd])  {
1305           if (fHisto2d) {
1306             fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
1307           }
1308           if (fVector2d) {
1309             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd+1]/fRelativeScale);
1310           }
1311         }
1312       }
1313     } // Case 2 zones
1314   }
1315
1316   // Fill with no automatic scale
1317   if (!fRelativeScaleAuto) {
1318     // Case of track with only one zone
1319     if (nb == 1) {
1320       fNumberUsedCh[0]++;
1321       if (fHisto2d) {
1322         fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1323       }
1324       if (fVector2d) {
1325         fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1326       }
1327     } // Case 1 zone
1328     // Case of track with two zones
1329     if (nb == 2) {
1330       // Two zones voisines sinon rien!
1331       // Case 1
1332       if ((fAmpTotal[fd]   > 0.0) && 
1333           (fAmpTotal[fd+1] > 0.0)) {
1334         // One of the two very big
1335         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1336           if (fHisto2d) {
1337             fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1338           }
1339           if (fVector2d) {
1340             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1341           }
1342           fNumberUsedCh[1]++;
1343         }
1344         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1345           if (fHisto2d) {
1346             fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
1347           }
1348           if (fVector2d) {
1349             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1350           }
1351           fNumberUsedCh[1]++;
1352         }
1353       }
1354       // Case 2
1355       if (fCalibraMode->GetNfragZ(0) > 1) {
1356         if (fAmpTotal[fd] > 0.0) {
1357           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1358             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1359               // One of the two very big
1360               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1361                 if (fHisto2d) {
1362                   fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+0.5,fAmpTotal[fd]/fRelativeScale);
1363                 }
1364                 if (fVector2d) {
1365                   fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1366                 }
1367                 fNumberUsedCh[1]++;
1368               }
1369               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1370                 if (fHisto2d) {
1371                   fCH2d->Fill(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)
1372                             + 0.5,fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1373                 }
1374                 fNumberUsedCh[1]++;
1375                 if (fVector2d) {
1376                   fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)
1377                                                 ,fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1378                 }
1379               }
1380             }
1381           }
1382         }
1383       }
1384     } // Case 2 zones
1385
1386   }
1387
1388 }
1389
1390 //____________Offine tracking in the AliTRDtracker_____________________________
1391 void AliTRDCalibraFillHisto::ResetfVariables()
1392 {
1393   //
1394   // Reset values of fAmpTotal, fPHValue and fPHPlace for
1395   // the updateHistogram... functions
1396   //
1397
1398   // Reset the good track
1399   fGoodTrack = kTRUE;
1400   
1401   // Reset the fAmpTotal where we put value
1402   if (fCH2dOn) {
1403     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1404       fAmpTotal[k] = 0.0;
1405     }
1406   }
1407   
1408   // Reset the fPHValue
1409   if (fPH2dOn) {
1410     for (Int_t k = 0; k < fTimeMax; k++) {
1411       fPHValue[k] = 0.0;
1412       fPHPlace[k] = -1;
1413     }
1414   }
1415
1416 }
1417
1418 //____________Offine tracking in the AliTRDtracker_____________________________
1419 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1420 {
1421   //
1422   // For the offline tracking or mcm tracklets
1423   // This function will be called in the functions UpdateHistogram... 
1424   // to fill the info of a track for the drift velocity  calibration
1425   //
1426     
1427   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1428   Int_t fd1 = -1; // Premiere zone non nulle
1429   Int_t fd2 = -1; // Deuxieme zone non nulle
1430   Int_t k1  = -1; // Debut de la premiere zone
1431   Int_t k2  = -1; // Debut de la seconde zone
1432
1433   // See if the track goes through different zones
1434   for (Int_t k = 0; k < fTimeMax; k++) {
1435     if (fPHValue[k] > 0.0) {
1436       if (fd1 == -1) {
1437         fd1 = fPHPlace[k];
1438         k1  = k;              
1439       }
1440       if (fPHPlace[k] != fd1) {
1441         if (fd2 == -1) {
1442           k2  = k;
1443           fd2 = fPHPlace[k];
1444           nb  = 2;
1445         }
1446         if (fPHPlace[k] != fd2) {
1447           nb = 3;
1448         }
1449       }
1450     }
1451   }
1452   
1453   // Fill 
1454   // Case of track with only one zone
1455   if (nb == 1) {
1456     fNumberUsedPh[0]++;
1457     //fd1 is the only zone
1458     for (Int_t i = 0; i < fTimeMax; i++) {
1459       if (fHisto2d) {
1460         fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1461       }
1462       if (fVector2d) {
1463         fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1464       }
1465     }
1466   } // Case 1 zone
1467   // Case of track with two zones
1468   if (nb == 2) {
1469     // Two zones voisines sinon rien!
1470     // Case 1
1471     if ((fd1 == fd2+1) || 
1472         (fd2 == fd1+1)) {
1473       // One of the two fast all the think
1474       if (k2 > (k1+fDifference)) {
1475         //we choose to fill the fd1 with all the values
1476         fNumberUsedPh[1]++;
1477         for (Int_t i = 0; i < fTimeMax; i++) {
1478           if (fHisto2d) {
1479             fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1480           }
1481           if (fVector2d) {
1482             fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1483           }
1484         }
1485       }
1486       if ((k2+fDifference) < fTimeMax) {
1487         //we choose to fill the fd2 with all the values
1488         fNumberUsedPh[1]++;
1489         for (Int_t i = 0; i < fTimeMax; i++) {
1490           if (fHisto2d) {
1491             fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1492           }
1493           if (fVector2d) {
1494             fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1495           }
1496         }
1497       }
1498     }
1499     // Two zones voisines sinon rien!
1500     if (fCalibraMode->GetNfragZ(1) > 1) {
1501       // Case 2
1502       if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1503         if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1504           // One of the two fast all the think
1505           if (k2 > (k1+fDifference)) {
1506             //we choose to fill the fd1 with all the values
1507             fNumberUsedPh[1]++;
1508             for (Int_t i = 0; i < fTimeMax; i++) {
1509               if (fHisto2d) {
1510                 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1511               }
1512               if (fVector2d) {
1513                 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1514               }
1515             }
1516           }
1517           if ((k2+fDifference) < fTimeMax) {
1518             //we choose to fill the fd2 with all the values
1519             fNumberUsedPh[1]++;
1520             for (Int_t i = 0; i < fTimeMax; i++) {
1521               if (fHisto2d) {
1522                 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1523               }
1524               if (fVector2d) {
1525                 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1526               }
1527             }
1528           }
1529         }
1530       }
1531       // Two zones voisines sinon rien!
1532       // Case 3
1533       if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1534         if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1535           // One of the two fast all the think
1536           if (k2 > (k1 + fDifference)) {
1537             //we choose to fill the fd1 with all the values
1538             fNumberUsedPh[1]++;
1539             for (Int_t i = 0; i < fTimeMax; i++) {
1540               if (fHisto2d) {
1541                 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1542               }
1543               if (fVector2d) {
1544                 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1545               }
1546             }
1547           }
1548           if ((k2+fDifference) < fTimeMax) {
1549             //we choose to fill the fd2 with all the values
1550             fNumberUsedPh[1]++;
1551             for (Int_t i = 0; i < fTimeMax; i++) {
1552               if (fHisto2d) {
1553                 fPH2d->Fill((fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
1554               }
1555               if (fVector2d) {
1556                 fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1557               }
1558             }
1559           }
1560         }
1561       }
1562     }
1563
1564   } // case 2 zones
1565
1566 }
1567
1568 //____________Set the pad calibration variables for the detector_______________
1569 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1570 {
1571   //
1572   // For the detector calcul the first Xbins and set the number of row
1573   // and col pads per calibration groups, the number of calibration
1574   // groups in the detector.
1575   //
1576   
1577   // first Xbins of the detector
1578   if (fCH2dOn) {
1579     fCalibraMode->CalculXBins(detector,0);
1580   }
1581   if (fPH2dOn) {
1582     fCalibraMode->CalculXBins(detector,1);
1583   }
1584   if (fPRF2dOn) {
1585     fCalibraMode->CalculXBins(detector,2);
1586   }
1587
1588   // fragmentation of idect
1589   for (Int_t i = 0; i < 3; i++) {
1590     fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
1591     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
1592                        , (Int_t) GetChamber(detector)
1593                        , (Int_t) GetSector(detector),i);
1594   }
1595   
1596   return kTRUE;
1597
1598 }
1599
1600 //
1601 //____________Some basic geometry function_____________________________________
1602 //
1603
1604 //_____________________________________________________________________________
1605 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
1606 {
1607   //
1608   // Reconstruct the plane number from the detector number
1609   //
1610
1611   return ((Int_t) (d % 6));
1612
1613 }
1614
1615 //_____________________________________________________________________________
1616 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
1617 {
1618   //
1619   // Reconstruct the chamber number from the detector number
1620   //
1621   Int_t fgkNplan = 6;
1622
1623   return ((Int_t) (d % 30) / fgkNplan);
1624
1625 }
1626
1627 //_____________________________________________________________________________
1628 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
1629 {
1630   //
1631   // Reconstruct the sector number from the detector number
1632   //
1633   Int_t fg = 30;
1634
1635   return ((Int_t) (d / fg));
1636
1637 }