Move pad planes from AliTRDCommomParam to AliTRDgeometry
[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 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 <TTree.h>
37 #include <TProfile2D.h>
38 #include <TProfile.h>
39 #include <TFile.h>
40 #include <TChain.h>
41 #include <TStyle.h>
42 #include <TCanvas.h>
43 #include <TGraphErrors.h>
44 #include <TObjArray.h>
45 #include <TH1F.h>
46 #include <TH2I.h>
47 #include <TH2.h>
48 #include <TStopwatch.h>
49 #include <TMath.h>
50 #include <TDirectory.h>
51 #include <TROOT.h>
52 #include <TTreeStream.h>
53 #include <TVectorD.h>
54
55 #include "AliLog.h"
56 #include "AliCDBManager.h"
57 #include "AliRawReader.h"
58 #include "AliRawReaderDate.h"
59
60 #include "AliTRDCalibraFillHisto.h"
61 #include "AliTRDCalibraMode.h"
62 #include "AliTRDCalibraVector.h"
63 #include "AliTRDcalibDB.h"
64 #include "AliTRDCommonParam.h"
65 #include "AliTRDmcmTracklet.h"
66 #include "AliTRDpadPlane.h"
67 #include "AliTRDcluster.h"
68 #include "AliTRDtrack.h"
69 #include "AliTRDRawStream.h"
70 #include "AliTRDgeometry.h"
71
72 #ifdef ALI_DATE
73 #include "event.h"
74 #endif
75
76 ClassImp(AliTRDCalibraFillHisto)
77
78 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
79 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
80
81 //_____________singleton implementation_________________________________________________
82 AliTRDCalibraFillHisto *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 //______________________________________________________________________________________
101 void 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 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
118   :TObject()
119   ,fGeo(0)
120   ,fMITracking(kFALSE)
121   ,fMcmTracking(kFALSE)
122   ,fMcmCorrectAngle(kFALSE)
123   ,fCH2dOn(kFALSE)
124   ,fPH2dOn(kFALSE)
125   ,fPRF2dOn(kFALSE)
126   ,fHisto2d(kFALSE)
127   ,fVector2d(kFALSE)
128   ,fLinearFitterOn(kFALSE)
129   ,fLinearFitterDebugOn(kFALSE)
130   ,fRelativeScale(0)
131   ,fThresholdClusterPRF2(15.0)
132   ,fCalibraMode(new AliTRDCalibraMode())
133   ,fDebugStreamer(0)
134   ,fDebugLevel(0)
135   ,fDetectorAliTRDtrack(kFALSE)
136   ,fDetectorPreviousTrack(-1)
137   ,fNumberClusters(18)
138   ,fProcent(6.0)
139   ,fDifference(17)
140   ,fNumberTrack(0)
141   ,fTimeMax(0)
142   ,fSf(10.0)
143   ,fNumberBinCharge(100)
144   ,fNumberBinPRF(40)
145   ,fNgroupprf(0)
146   ,fListClusters(new TObjArray()) 
147   ,fPar0(0x0)
148   ,fPar1(0x0)
149   ,fPar2(0x0)
150   ,fPar3(0x0)
151   ,fPar4(0x0)
152   ,fAmpTotal(0x0)
153   ,fPHPlace(0x0)
154   ,fPHValue(0x0)
155   ,fGoodTracklet(kTRUE)
156   ,fGoodTrack(kTRUE)
157   ,fEntriesCH(0x0)
158   ,fEntriesLinearFitter(0x0)
159   ,fCalibraVector(0x0)
160   ,fPH2d(0x0)
161   ,fPRF2d(0x0)
162   ,fCH2d(0x0)
163   ,fLinearFitterArray(0)
164   ,fLinearFitterHistoArray(0)
165 {
166   //
167   // Default constructor
168   //
169
170   //
171   // Init some default values
172   //
173
174   fNumberUsedCh[0]       = 0;
175   fNumberUsedCh[1]       = 0;
176   fNumberUsedPh[0]       = 0;
177   fNumberUsedPh[1]       = 0;
178
179   fGeo = new AliTRDgeometry();
180  
181 }
182
183 //______________________________________________________________________________________
184 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
185   :TObject(c)
186   ,fGeo(0)
187   ,fMITracking(c.fMITracking)
188   ,fMcmTracking(c.fMcmTracking)
189   ,fMcmCorrectAngle(c.fMcmCorrectAngle)
190   ,fCH2dOn(c.fCH2dOn)
191   ,fPH2dOn(c.fPH2dOn)
192   ,fPRF2dOn(c.fPRF2dOn)
193   ,fHisto2d(c.fHisto2d)
194   ,fVector2d(c.fVector2d)
195   ,fLinearFitterOn(c.fLinearFitterOn)
196   ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
197   ,fRelativeScale(c.fRelativeScale)
198   ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
199   ,fCalibraMode(0x0)
200   ,fDebugStreamer(0)
201   ,fDebugLevel(c.fDebugLevel)
202   ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
203   ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
204   ,fNumberClusters(c.fNumberClusters)
205   ,fProcent(c.fProcent)
206   ,fDifference(c.fDifference)
207   ,fNumberTrack(c.fNumberTrack)
208   ,fTimeMax(c.fTimeMax)
209   ,fSf(c.fSf)
210   ,fNumberBinCharge(c.fNumberBinCharge)
211   ,fNumberBinPRF(c.fNumberBinPRF)
212   ,fNgroupprf(c.fNgroupprf)
213   ,fListClusters(new TObjArray())
214   ,fPar0(c.fPar0)
215   ,fPar1(c.fPar1)
216   ,fPar2(c.fPar2)
217   ,fPar3(c.fPar3)
218   ,fPar4(c.fPar4)
219   ,fAmpTotal(c.fAmpTotal)
220   ,fPHPlace(c.fPHPlace)
221   ,fPHValue(c.fPHValue)
222   ,fGoodTracklet(c.fGoodTracklet)
223   ,fGoodTrack(c.fGoodTrack)
224   ,fEntriesCH(c.fEntriesCH)
225   ,fEntriesLinearFitter(fEntriesLinearFitter)
226   ,fCalibraVector(0x0)
227   ,fPH2d(0x0)
228   ,fPRF2d(0x0)
229   ,fCH2d(0x0)
230   ,fLinearFitterArray(0)
231   ,fLinearFitterHistoArray(0)
232 {
233   //
234   // Copy constructor
235   //
236   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
237   if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
238   if(c.fPH2d) {
239     fPH2d = new TProfile2D(*c.fPH2d);
240     fPH2d->SetDirectory(0);
241   }
242   if(c.fPRF2d) {
243     fPRF2d = new TProfile2D(*c.fPRF2d);
244     fPRF2d->SetDirectory(0);
245   }
246   if(c.fCH2d) {
247     fCH2d = new TH2I(*c.fCH2d);
248     fCH2d->SetDirectory(0);
249   }
250   if(c.fLinearFitterOn){
251     fLinearFitterArray.Expand(540);
252     for (Int_t cb = 0; cb < 540; cb++){
253       const TLinearFitter *linearFitter = (TLinearFitter*)c.fLinearFitterArray.UncheckedAt(cb);
254       if ( linearFitter != 0x0 ) fLinearFitterArray.AddAt(new TLinearFitter(*linearFitter), cb);
255     }
256   }
257   if(c.fLinearFitterDebugOn){
258     fLinearFitterHistoArray.Expand(540);
259     for (Int_t cb = 0; cb < 540; cb++){
260       const TH2F *linearfitterhisto = (TH2F*)c.fLinearFitterHistoArray.UncheckedAt(cb);
261       if ( linearfitterhisto != 0x0 ){
262         TH2F *hNew = new TH2F(*linearfitterhisto);
263         hNew->SetDirectory(0);
264         fLinearFitterHistoArray.AddAt(hNew,cb);
265       }
266     }
267   }
268   if (fGeo) {
269     delete fGeo;
270   }
271   fGeo = new AliTRDgeometry();
272
273 }
274
275 //____________________________________________________________________________________
276 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
277 {
278   //
279   // AliTRDCalibraFillHisto destructor
280   //
281
282   ClearHistos();
283   if ( fDebugStreamer ) delete fDebugStreamer;
284
285   if (fGeo) {
286     delete fGeo;
287   }
288   
289 }
290
291 //_____________________________________________________________________________
292 void AliTRDCalibraFillHisto::Destroy() 
293 {
294   //
295   // Delete instance 
296   //
297
298   if (fgInstance) {
299     delete fgInstance;
300     fgInstance = 0x0;
301   }
302
303 }
304
305 //_____________________________________________________________________________
306 void AliTRDCalibraFillHisto::ClearHistos() 
307 {
308   //
309   // Delete the histos
310   //
311
312   if (fPH2d) {
313     delete fPH2d;
314     fPH2d  = 0x0;
315   }
316   if (fCH2d) {
317     delete fCH2d;
318     fCH2d  = 0x0;
319   }
320   if (fPRF2d) {
321     delete fPRF2d;
322     fPRF2d = 0x0;
323   }
324
325 }
326 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
327 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
328 {
329   //
330   // For the offline tracking
331   // This function will be called in the function AliReconstruction::Run() 
332   // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE, 
333   //
334
335   // DB Setting
336   // Get cal
337   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
338   if (!cal) {
339     AliInfo("Could not get calibDB");
340     return kFALSE;
341   }
342   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
343   if (!parCom) {
344     AliInfo("Could not get CommonParam");
345     return kFALSE;
346   }
347
348   // Some parameters
349   fTimeMax = cal->GetNumberOfTimeBins();
350   fSf      = parCom->GetSamplingFrequency();
351   fRelativeScale = 20;
352
353   //Init the tracklet parameters
354   fPar0 = new Double_t[fTimeMax];
355   fPar1 = new Double_t[fTimeMax];
356   fPar2 = new Double_t[fTimeMax];
357   fPar3 = new Double_t[fTimeMax];
358   fPar4 = new Double_t[fTimeMax];
359   
360   for(Int_t k = 0; k < fTimeMax; k++){
361     fPar0[k] = 0.0;
362     fPar1[k] = 0.0;
363     fPar2[k] = 0.0;
364     fPar3[k] = 0.0;
365     fPar4[k] = 0.0;
366   }
367   
368   //If vector method On initialised all the stuff
369   if(fVector2d){
370     fCalibraVector = new AliTRDCalibraVector();
371     fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
372     fCalibraVector->SetTimeMax(fTimeMax);
373     if(fNgroupprf != 0) {
374       fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
375       fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
376     }
377     else {
378       fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
379       fCalibraVector->SetPRFRange(1.5);
380     }
381   }
382  
383   // Create the 2D histos corresponding to the pad groupCalibration mode
384   if (fCH2dOn) {
385
386     AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
387                 ,fCalibraMode->GetNz(0)
388                 ,fCalibraMode->GetNrphi(0)));
389     
390     // Calcul the number of Xbins
391     Int_t Ntotal0 = CalculateTotalNumberOfBins(0);
392     // Create the 2D histo
393     if (fHisto2d) {
394       CreateCH2d(Ntotal0);
395     }
396     // Variable
397     fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
398     for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
399       fAmpTotal[k] = 0.0;
400     } 
401     //Statistics
402     fEntriesCH = new Int_t[Ntotal0];
403     for(Int_t k = 0; k < Ntotal0; k++){
404       fEntriesCH[k] = 0;
405     }
406   }
407   if (fPH2dOn) {
408
409     AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
410                 ,fCalibraMode->GetNz(1)
411                 ,fCalibraMode->GetNrphi(1)));
412     
413     // Calcul the number of Xbins
414     Int_t Ntotal1 = CalculateTotalNumberOfBins(1);
415     // Create the 2D histo
416     if (fHisto2d) {
417       CreatePH2d(Ntotal1);
418     }
419     // Variable
420     fPHPlace = new Short_t[fTimeMax];
421     for (Int_t k = 0; k < fTimeMax; k++) {
422       fPHPlace[k] = -1;
423     } 
424     fPHValue = new Float_t[fTimeMax];
425     for (Int_t k = 0; k < fTimeMax; k++) {
426       fPHValue[k] = 0.0;
427     }
428   }
429   if (fLinearFitterOn) {
430     fLinearFitterArray.Expand(540);
431     fLinearFitterArray.SetName("ArrayLinearFitters");
432     fEntriesLinearFitter = new Int_t[540];
433     for(Int_t k = 0; k < 540; k++){
434       fEntriesLinearFitter[k] = 0;
435     }
436     if(fLinearFitterDebugOn) {
437       fLinearFitterHistoArray.Expand(540);
438       fLinearFitterHistoArray.SetName("ArrayHistos");
439     }
440   }
441
442   if (fPRF2dOn) {
443
444     AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
445                 ,fCalibraMode->GetNz(2)
446                 ,fCalibraMode->GetNrphi(2)));
447     
448     // Calcul the number of Xbins
449     Int_t Ntotal2 = CalculateTotalNumberOfBins(2);
450     // Create the 2D histo
451     if (fHisto2d) {
452       CreatePRF2d(Ntotal2);
453     }
454   }
455
456   return kTRUE;
457
458 }
459
460 //____________Functions for filling the histos in the code_____________________
461
462 //____________Offine tracking in the AliTRDtracker_____________________________
463 Bool_t AliTRDCalibraFillHisto::ResetTrack()
464 {
465   //
466   // For the offline tracking
467   // This function will be called in the function
468   // AliTRDtracker::FollowBackPropagation() at the beginning 
469   // Reset the parameter to know we have a new TRD track
470   //
471   
472   fDetectorAliTRDtrack = kFALSE;
473   fGoodTrack           = kTRUE;
474   return kTRUE;
475
476 }
477 //____________Offine tracking in the AliTRDtracker_____________________________
478 void AliTRDCalibraFillHisto::ResetfVariables()
479 {
480   //
481   // Reset values per tracklet
482   //
483
484   // Reset the list of clusters
485   fListClusters->Clear();
486
487   //Reset the tracklet parameters
488   for(Int_t k = 0; k < fTimeMax; k++){
489     fPar0[k] = 0.0;
490     fPar1[k] = 0.0;
491     fPar2[k] = 0.0;
492     fPar3[k] = 0.0; 
493     fPar4[k] = 0.0;
494   }
495
496     
497   //Reset good tracklet
498   fGoodTracklet = kTRUE;
499
500   // Reset the fPHValue
501   if (fPH2dOn) {
502     //Reset the fPHValue and fPHPlace
503     for (Int_t k = 0; k < fTimeMax; k++) {
504       fPHValue[k] = 0.0;
505       fPHPlace[k] = -1;
506     }
507   }
508
509   // Reset the fAmpTotal where we put value
510   if (fCH2dOn) {
511     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
512       fAmpTotal[k] = 0.0;
513     }
514   }
515
516 }
517 //____________Offline tracking in the AliTRDtracker____________________________
518 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
519 {
520   //
521   // For the offline tracking
522   // This function will be called in the function
523   // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
524   // of TRD tracks 
525   // Fill the 2D histos or the vectors with the info of the clusters at
526   // the end of a detectors if the track is "good"
527   //
528
529   // Localisation of the detector
530   Int_t detector = cl->GetDetector();
531
532   // Fill the infos for the previous clusters if not the same
533   // detector anymore or if not the same track
534   if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) && 
535       (fDetectorPreviousTrack != -1)) {
536     
537     fNumberTrack++;   
538     
539     // If the same track, then look if the previous detector is in
540     // the same plane, if yes: not a good track
541     if (fDetectorAliTRDtrack && 
542         (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
543       fGoodTrack = kFALSE;
544     }
545
546     // Fill only if the track doesn't touch a masked pad or doesn't
547     // appear in the middle (fGoodTrack)
548     if (fGoodTrack && fGoodTracklet) {
549
550       // drift velocity unables to cut bad tracklets 
551       Bool_t  pass = FindP1TrackPH();
552       
553       // Gain calibration
554       if (fCH2dOn) {
555         FillTheInfoOfTheTrackCH();
556       }
557       
558       // PH calibration
559       if (fPH2dOn) {
560         FillTheInfoOfTheTrackPH();    
561       }
562
563       if(pass && fPRF2dOn) HandlePRF();
564
565       
566     } // if a good track
567     
568     ResetfVariables();
569     
570   } // Fill at the end the charge
571   
572   // Calcul the position of the detector
573   if (detector != fDetectorPreviousTrack) {
574     LocalisationDetectorXbins(detector);
575   }
576  
577   // Reset the detector
578   fDetectorPreviousTrack = detector;
579   fDetectorAliTRDtrack   = kTRUE;
580
581   // Store the infos of the tracklets
582   AliTRDcluster *kcl = new AliTRDcluster(*cl);
583   fListClusters->Add((TObject *)kcl);
584   Int_t time = cl->GetLocalTimeBin();
585   fPar0[time] = t->GetY();
586   fPar1[time] = t->GetZ();
587   fPar2[time] = t->GetSnp();
588   fPar3[time] = t->GetTgl();
589   fPar4[time] = t->Get1Pt();
590
591   // Store the info bis of the tracklet
592   Int_t *rowcol   = CalculateRowCol(cl);
593   CheckGoodTracklet(detector,rowcol);
594   Int_t     group[2] = {0,0};
595   if(fCH2dOn)  group[0]  = CalculateCalibrationGroup(0,rowcol);
596   if(fPH2dOn)  group[1]  = CalculateCalibrationGroup(1,rowcol);
597   StoreInfoCHPH(cl,t,group);
598    
599   return kTRUE;
600   
601 }
602 //____________Online trackling in AliTRDtrigger________________________________
603 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
604 {
605   //
606   // For the tracking
607   // This function will be called in the function AliTRDtrigger::TestTracklet
608   // before applying the pt cut on the tracklets 
609   // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
610   //
611   
612   // Localisation of the Xbins involved
613   Int_t idect = trk->GetDetector();
614   LocalisationDetectorXbins(idect);
615
616   // Get the parameter object
617   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
618   if (!cal) {
619     AliInfo("Could not get calibDB");
620     return kFALSE;
621   }
622    
623   // Reset
624   ResetfVariables();
625  
626   // Row of the tracklet and position in the pad groups
627   Int_t *rowcol  = new Int_t[2];
628   rowcol[0]     = trk->GetRow();
629   Int_t group[3] = {-1,-1,-1};
630   
631   // Eventuelle correction due to track angle in z direction
632   Float_t correction = 1.0;
633   if (fMcmCorrectAngle) {
634     Float_t z = trk->GetRowz();
635     Float_t r = trk->GetTime0();
636     correction = r / TMath::Sqrt((r*r+z*z));
637   }
638
639   // Boucle sur les clusters
640   // Condition on number of cluster: don't come from the middle of the detector
641   if (trk->GetNclusters() >= fNumberClusters) {
642
643     for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
644
645       Float_t amp[3] = { 0.0, 0.0, 0.0 };
646       Int_t   time   = trk->GetClusterTime(icl);
647       rowcol[1]      = trk->GetClusterCol(icl);
648             
649       amp[0] = trk->GetClusterADC(icl)[0] * correction;
650       amp[1] = trk->GetClusterADC(icl)[1] * correction;
651       amp[2] = trk->GetClusterADC(icl)[2] * correction;
652
653       
654       if ((amp[0] < 0.0) || 
655           (amp[1] < 0.0) || 
656           (amp[2] < 0.0)) {
657         continue;
658       }
659
660       // Col of cluster and position in the pad groups
661       if(fCH2dOn)  {
662         group[0] = CalculateCalibrationGroup(0,rowcol);
663         fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
664       }
665       if(fPH2dOn)  {
666         group[1] = CalculateCalibrationGroup(1,rowcol);
667         fPHPlace[time] = group[1];
668         fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
669       }
670       if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
671
672       // See if we are not near a masked pad fGoodTracklet
673       CheckGoodTracklet(idect,rowcol);
674                
675       // Fill PRF direct without tnp bins...only for monitoring...
676       if (fPRF2dOn && fGoodTracklet) {
677         
678         if ((amp[0] > fThresholdClusterPRF2) && 
679             (amp[1] > fThresholdClusterPRF2) && 
680             (amp[2] > fThresholdClusterPRF2) && 
681             ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
682         
683           // Security of the denomiateur is 0
684           if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1]))) 
685              / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
686             Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
687                                   / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
688             Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
689
690             if (TMath::Abs(xcenter) < 0.5) {
691               Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
692               Float_t ymax   = amp[2] / (amp[0]+amp[1]+amp[2]);
693               // Fill only if it is in the drift region!
694               if (((Float_t) time / fSf) > 0.3) {
695                 if (fHisto2d) {
696                   fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
697                   fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
698                   fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
699                 }
700                 if (fVector2d) {
701                   fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],xcenter,ycenter);
702                   fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],-(xcenter+1.0),yminus);
703                   fCalibraVector->UpdateVectorPRF(fCalibraMode->GetXbins(2)+group[2],(1.0-xcenter),ymax);
704                 }
705               }//in the drift region 
706             }//in the middle
707           }//denominateur security
708         }//cluster shape and thresholds
709       }//good and PRF On
710       
711     } // Boucle clusters
712     
713     // Fill the charge
714     if(fGoodTracklet){
715       if (fCH2dOn) FillTheInfoOfTheTrackCH();
716       if (fPH2dOn) FillTheInfoOfTheTrackPH();   
717     }
718
719     fNumberTrack++;
720         
721   } // Condition on number of clusters
722
723   return kTRUE;
724   
725 }
726 //_____________________________________________________________________________
727 Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
728 {
729   //
730   // Calculate the row and col number of the cluster
731   //
732
733
734   Int_t *rowcol = new Int_t[2];
735   rowcol[0] =  0;
736   rowcol[1] =  0;
737
738   // Localisation of the detector
739   Int_t detector = cl->GetDetector();
740   Int_t chamber  = GetChamber(detector);
741   Int_t plane    = GetPlane(detector);
742
743   // Localisation of the cluster
744   Double_t pos[3] = { 0.0, 0.0, 0.0 };
745   pos[0] = ((AliCluster *)cl)->GetX();
746   pos[1] = cl->GetY();
747   pos[2] = cl->GetZ();
748
749   // Position of the cluster
750   AliTRDpadPlane *padplane  = fGeo->GetPadPlane(plane,chamber);
751   Int_t    row              = padplane->GetPadRowNumber(pos[2]);
752   //Do not take from here because it was corrected from ExB already....
753   //Double_t offsetz         = padplane->GetPadRowOffset(row,pos[2]);
754   //Double_t offsettilt      = padplane->GetTiltOffset(offsetz);
755   //Int_t    col             = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
756   //Int_t    col             = padplane->GetPadColNumber(pos[1]+offsettilt);
757   Int_t    col               = cl->GetPad(); 
758
759   //return
760   rowcol[0]     = row;
761   rowcol[1]     = col; 
762   return rowcol;
763   
764 }
765 //_____________________________________________________________________________
766 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
767 {
768   //
769   // See if we are not near a masked pad
770   //
771
772   Int_t row = rowcol[0];
773   Int_t col = rowcol[1];
774
775   if (!IsPadOn(detector, col, row)) {
776     fGoodTracklet = kFALSE;
777   }
778
779   if (col > 0) {
780     if (!IsPadOn(detector, col-1, row)) {
781       fGoodTracklet = kFALSE;
782     }
783   }
784
785   if (col < 143) {
786     if (!IsPadOn(detector, col+1, row)) {
787       fGoodTracklet = kFALSE;
788     }
789   }
790   
791 }
792 //_____________________________________________________________________________
793 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
794 {
795   //
796   // Look in the choosen database if the pad is On.
797   // If no the track will be "not good"
798   //
799
800   // Get the parameter object
801   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
802   if (!cal) {
803     AliInfo("Could not get calibDB");
804     return kFALSE;
805   }
806   
807   if (!cal->IsChamberInstalled(detector)     || 
808        cal->IsChamberMasked(detector)        ||
809        cal->IsPadMasked(detector,col,row)) {
810     return kFALSE;
811   }
812   else {
813     return kTRUE;
814   }
815   
816 }
817 //_____________________________________________________________________________
818 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
819 {
820   //
821   // Calculate the calibration group number for i
822   //
823  
824   // Row of the cluster and position in the pad groups
825   Int_t posr = 0;
826   if (fCalibraMode->GetNnZ(i) != 0) {
827     posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
828   }
829  
830       
831   // Col of the cluster and position in the pad groups
832   Int_t posc = 0;
833   if (fCalibraMode->GetNnRphi(i) != 0) {
834     posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
835   }
836   
837   return posc*fCalibraMode->GetNfragZ(i)+posr;
838   
839 }
840 //____________________________________________________________________________________
841 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
842 {
843   //
844   // Calculate the total number of calibration groups
845   //
846   
847   Int_t Ntotal = 0;
848   fCalibraMode->ModePadCalibration(2,i);
849   fCalibraMode->ModePadFragmentation(0,2,0,i);
850   fCalibraMode->SetDetChamb2(i);
851   Ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
852   fCalibraMode->ModePadCalibration(0,i);
853   fCalibraMode->ModePadFragmentation(0,0,0,i);
854   fCalibraMode->SetDetChamb0(i);
855   Ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
856   AliInfo(Form("Total number of Xbins: %d",Ntotal));
857   return Ntotal;
858
859 }
860 //____________Set the pad calibration variables for the detector_______________
861 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
862 {
863   //
864   // For the detector calcul the first Xbins and set the number of row
865   // and col pads per calibration groups, the number of calibration
866   // groups in the detector.
867   //
868   
869   // first Xbins of the detector
870   if (fCH2dOn) {
871     fCalibraMode->CalculXBins(detector,0);
872   }
873   if (fPH2dOn) {
874     fCalibraMode->CalculXBins(detector,1);
875   }
876   if (fPRF2dOn) {
877     fCalibraMode->CalculXBins(detector,2);
878   }
879
880   // fragmentation of idect
881   for (Int_t i = 0; i < 3; i++) {
882     fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
883     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
884                        , (Int_t) GetChamber(detector)
885                        , (Int_t) GetSector(detector),i);
886   }
887   
888   return kTRUE;
889
890 }
891 //_____________________________________________________________________________
892 void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group)
893 {
894   //
895   // Store the infos in fAmpTotal, fPHPlace and fPHValue
896   //
897   
898   // Charge in the cluster
899   Float_t  q        = TMath::Abs(cl->GetQ());
900   Int_t    time     = cl->GetLocalTimeBin();
901    
902   // Correction due to the track angle
903   Float_t correction    = 1.0;
904   Float_t normalisation = 6.67;
905   if ((q >0) && (t->GetNdedx() > 0)) {
906     correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
907   }
908
909   // Fill the fAmpTotal with the charge
910   if (fCH2dOn) {
911     fAmpTotal[(Int_t) group[0]] += correction;
912   }
913
914   // Fill the fPHPlace and value
915   if (fPH2dOn) {
916     fPHPlace[time] = group[1];
917     fPHValue[time] = correction;
918   }
919   
920 }
921 //_____________________________________________________________________
922 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDRawStream *rawStream)
923 {
924   //
925   // Event Processing loop - AliTRDRawStream
926   //
927
928   Bool_t withInput = kFALSE;
929
930   Int_t phvalue[36];
931   //Int_t row[36];
932   //Int_t col[36];
933   for(Int_t k = 0; k < 36; k++){
934     phvalue[k] = 0;
935     //row[k]     = -1;
936     //col[36]    = -1;
937   }
938   fDetectorPreviousTrack = -1;
939   Int_t nbtimebin = 0;                                           
940
941   while (rawStream->Next()) {
942
943     Int_t idetector = rawStream->GetDet();                            //  current detector
944     if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
945       if(TMath::Mean(nbtimebin,phvalue)>0.0){
946         withInput = kTRUE;
947         for(Int_t k = 0; k < nbtimebin; k++){
948           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
949           phvalue[k] = 0;
950           //row[k]     = -1;
951           //col[k]     = -1;
952         }
953       }
954     }
955     fDetectorPreviousTrack = idetector;
956     nbtimebin         = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
957     Int_t iTimeBin    = rawStream->GetTimeBin();                       //  current time bin
958     //row[iTimeBin]   = rawStream->GetRow();                           //  current row
959     //col[iTimeBin]   = rawStream->GetCol();                           //  current col     
960     Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
961         
962     Int_t fin     = TMath::Min(nbtimebin,(iTimeBin+3));
963     Int_t n       = 0;
964     for(Int_t itime = iTimeBin; itime < fin; itime++){
965       // should extract baseline here!
966       if(signal[n]>5.0) phvalue[itime] = signal[n];
967       n++;
968     }
969   }
970   
971   // fill the last one
972   if(fDetectorPreviousTrack != -1){
973       if(TMath::Mean(nbtimebin,phvalue)>0.0){
974         withInput = kTRUE;
975         for(Int_t k = 0; k < nbtimebin; k++){
976           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
977           phvalue[k] = 0;
978           //row[k]     = -1;
979           //col[k]     = -1;
980         }
981       }
982   }
983   
984   return withInput;
985   
986 }
987 //_____________________________________________________________________
988 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
989 {
990   //
991   //  Event processing loop - AliRawReader
992   //
993
994
995   AliTRDRawStream rawStream(rawReader);
996
997   rawReader->Select("TRD");
998
999   return ProcessEventDAQ(&rawStream);
1000 }
1001 //_________________________________________________________________________
1002 Bool_t AliTRDCalibraFillHisto::ProcessEventDAQ(
1003 #ifdef ALI_DATE
1004                                    eventHeaderStruct *event
1005 #else
1006                                    eventHeaderStruct* /*event*/
1007             
1008 #endif 
1009                                    )
1010 {
1011   //
1012   //  process date event
1013   //
1014 #ifdef ALI_DATE
1015     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1016     Bool_t result=ProcessEventDAQ(rawReader);
1017     delete rawReader;
1018     return result;
1019 #else
1020     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1021     return 0;
1022 #endif
1023
1024 }
1025 //____________Online trackling in AliTRDtrigger________________________________
1026 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1027 {
1028   //
1029   // For the DAQ
1030   // Fill a simple average pulse height
1031   //
1032   
1033   // Localisation of the Xbins involved
1034   //LocalisationDetectorXbins(det);
1035
1036   // Row  and position in the pad groups
1037   //Int_t posr = 0;
1038   //if (fCalibraMode->GetNnZ(1) != 0) {
1039   //  posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1040   //}
1041  
1042   // Col of cluster and position in the pad groups
1043   //Int_t posc = 0;
1044   //if (fCalibraMode->GetNnRphi(1) != 0) {
1045   //  posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1046   //}
1047   
1048   //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);   
1049   
1050   ((TProfile2D *)GetPH2d(nbtimebins,fSf,kTRUE))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1051   
1052   return kTRUE;
1053   
1054 }
1055 //____________Functions for plotting the 2D____________________________________
1056
1057 //_____________________________________________________________________________
1058 void AliTRDCalibraFillHisto::Plot2d()
1059 {
1060   //
1061   // Plot the 2D histos 
1062   //
1063  
1064   if (fPH2dOn) {
1065     TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
1066     cph2d->cd();
1067     fPH2d->Draw("LEGO");
1068   }
1069   if (fCH2dOn) {
1070     TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
1071     cch2d->cd();
1072     fCH2d->Draw("LEGO");
1073   }
1074   if (fPRF2dOn) {
1075     TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
1076     cPRF2d->cd();
1077     fPRF2d->Draw("LEGO");
1078   }
1079
1080 }
1081 //_____________________Reset____________________________________________________
1082 //_______________________________________________________________________________
1083 void AliTRDCalibraFillHisto::ResetLinearFitter()
1084 {
1085   fLinearFitterArray.SetOwner();
1086   fLinearFitterArray.Clear();
1087   fLinearFitterHistoArray.SetOwner();
1088   fLinearFitterHistoArray.Clear();
1089 }
1090 //____________Write_____________________________________________________
1091 //_____________________________________________________________________
1092 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1093 {
1094     //
1095     //  Write infos to file
1096     //
1097
1098   //For debugging
1099   if ( fDebugStreamer ) {
1100     delete fDebugStreamer;
1101     fDebugStreamer = 0x0;
1102   }
1103
1104   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1105                ,fNumberTrack
1106                ,fNumberUsedCh[0]
1107                ,fNumberUsedCh[1]
1108                ,fNumberUsedPh[0]
1109                ,fNumberUsedPh[1]));
1110   
1111   TDirectory *backup = gDirectory;
1112   TString option;
1113   
1114   if ( append )
1115     option = "update";
1116   else
1117     option = "recreate";
1118   
1119   TFile f(filename,option.Data());
1120   
1121   TStopwatch stopwatch;
1122   stopwatch.Start();
1123
1124   if (fCH2dOn ) {
1125     if (fHisto2d) {
1126       f.WriteTObject(fCH2d);
1127     }
1128     if (fVector2d) {
1129       TString name("Nz");
1130       name += fCalibraMode->GetNz(0);
1131       name += "Nrphi";
1132       name += fCalibraMode->GetNrphi(0);
1133       TTree *treeCH2d = fCalibraVector->ConvertVectorCTTreeHisto(fCalibraVector->GetVectorCH(),fCalibraVector->GetPlaCH(),"treeCH2d",(const char *) name);
1134         f.WriteTObject(treeCH2d);
1135     }
1136   }
1137   if (fPH2dOn ) {
1138     if (fHisto2d) {
1139       f.WriteTObject(fPH2d);
1140     }
1141     if (fVector2d) {
1142       TString name("Nz");
1143       name += fCalibraMode->GetNz(1);
1144       name += "Nrphi";
1145       name += fCalibraMode->GetNrphi(1);
1146       TTree *treePH2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPH(),fCalibraVector->GetPlaPH(),"treePH2d",(const char *) name);
1147       f.WriteTObject(treePH2d);
1148     }
1149   }
1150   if (fPRF2dOn) {
1151     if (fHisto2d) {
1152         f.WriteTObject(fPRF2d);
1153     }
1154     if (fVector2d) {
1155       TString name("Nz");
1156       name += fCalibraMode->GetNz(2);
1157       name += "Nrphi";
1158       name += fCalibraMode->GetNrphi(2);
1159       name += "Ngp";
1160       name += fNgroupprf;
1161       TTree *treePRF2d = fCalibraVector->ConvertVectorPTreeHisto(fCalibraVector->GetVectorPRF(),fCalibraVector->GetPlaPRF(),"treePRF2d",(const char *) name);
1162       f.WriteTObject(treePRF2d);
1163     }
1164   }
1165   if(fLinearFitterOn && fLinearFitterDebugOn){
1166     f.WriteTObject(&fLinearFitterHistoArray);
1167   }
1168   
1169   f.Close();
1170   
1171   if ( backup ) backup->cd();
1172   
1173   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1174                ,stopwatch.RealTime(),stopwatch.CpuTime()));
1175 }
1176 //___________________________________________probe the histos__________________________________________________
1177 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1178 {
1179   //
1180   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1181   // debug mode with 2 for TH2I and 3 for TProfile2D
1182   // It gives a pointer to a Double_t[7] with the info following...
1183   // [0] : number of calibration groups with entries
1184   // [1] : minimal number of entries found
1185   // [2] : calibration group number of the min
1186   // [3] : maximal number of entries found
1187   // [4] : calibration group number of the max
1188   // [5] : mean number of entries found
1189   // [6] : mean relativ error
1190   //
1191
1192   Double_t *info = new Double_t[7];
1193    
1194   // Number of Xbins (detectors or groups of pads)
1195   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
1196   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
1197
1198   // Initialise
1199   Double_t nbwe = 0; //number of calibration groups with entries
1200   Double_t minentries = 0; //minimal number of entries found
1201   Double_t maxentries = 0; //maximal number of entries found
1202   Double_t placemin = 0; //calibration group number of the min
1203   Double_t placemax = -1; //calibration group number of the max
1204   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1205   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1206
1207   Double_t counter = 0;
1208
1209   //Debug
1210   TH1F *NbEntries = 0x0;//distribution of the number of entries
1211   TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1212   TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1213     
1214   // Beginning of the loop over the calibration groups 
1215   for (Int_t idect = 0; idect < nbins; idect++) {
1216
1217     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1218     projch->SetDirectory(0);
1219     
1220     // Number of entries for this calibration group
1221     Double_t nentries = 0.0;
1222     if((i%2) == 0){
1223       for (Int_t k = 0; k < nxbins; k++) {
1224         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1225       }
1226     }
1227     else{
1228       for (Int_t k = 0; k < nxbins; k++) {
1229         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1230         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1231           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1232           counter++;
1233         } 
1234       }
1235     }
1236
1237     //Debug
1238     if(i > 1){
1239       if((!((Bool_t)NbEntries)) && (nentries > 0)){
1240         NbEntries = new TH1F("Number of entries","Number of entries"
1241                                ,100,(Int_t)nentries/2,nentries*2);
1242         NbEntries->SetDirectory(0);
1243         NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1244                                ,nbins,0,nbins);
1245         NbEntriesPerGroup->SetDirectory(0);
1246         NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1247                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1248         NbEntriesPerSp->SetDirectory(0);
1249       }
1250       if(NbEntries){
1251         if(nentries > 0) NbEntries->Fill(nentries);
1252         NbEntriesPerGroup->Fill(idect+0.5,nentries);
1253         NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1254       }
1255     }
1256
1257     //min amd max
1258     if(nentries > maxentries){
1259       maxentries = nentries;
1260       placemax = idect;
1261     }
1262     if(idect == 0) {
1263       minentries = nentries;
1264     }
1265     if(nentries < minentries){
1266       minentries = nentries;
1267       placemin = idect;
1268     }
1269     //nbwe
1270     if(nentries > 0) {
1271       nbwe++;
1272       meanstats += nentries;
1273     }
1274   }//calibration groups loop
1275   
1276   if(nbwe > 0) meanstats /= nbwe;
1277   if(counter > 0) meanrelativerror /= counter;
1278
1279   AliInfo(Form("There are %f calibration groups with entries",nbwe));
1280   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1281   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1282   AliInfo(Form("The mean number of entries is %f",meanstats));
1283   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1284
1285   info[0] = nbwe;
1286   info[1] = minentries;
1287   info[2] = placemin;
1288   info[3] = maxentries;
1289   info[4] = placemax;
1290   info[5] = meanstats;
1291   info[6] = meanrelativerror;
1292
1293   if(i > 1){
1294     gStyle->SetPalette(1);
1295     gStyle->SetOptStat(1111);
1296     gStyle->SetPadBorderMode(0);
1297     gStyle->SetCanvasColor(10);
1298     gStyle->SetPadLeftMargin(0.13);
1299     gStyle->SetPadRightMargin(0.01);
1300     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1301     stat->Divide(2,1);
1302     stat->cd(1);
1303     NbEntries->Draw("");
1304     stat->cd(2);
1305     NbEntriesPerSp->SetStats(0);
1306     NbEntriesPerSp->Draw("");
1307     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1308     stat1->cd();
1309     NbEntriesPerGroup->SetStats(0);
1310     NbEntriesPerGroup->Draw("");
1311   }
1312
1313   return info;
1314
1315 }
1316 //____________________________________________________________________________
1317 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1318 {
1319   //
1320   // Return a Int_t[4] with:
1321   // 0 Mean number of entries
1322   // 1 median of number of entries
1323   // 2 rms of number of entries
1324   // 3 number of group with entries
1325   //
1326
1327   Double_t *stat      = new Double_t[4]; 
1328   stat[3]             = 0.0;
1329
1330   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
1331   Double_t *weight    = new Double_t[nbofgroups];
1332   Int_t    *nonul     = new Int_t[nbofgroups];
1333  
1334   for(Int_t k = 0; k < nbofgroups; k++){
1335     if(fEntriesCH[k] > 0) {
1336       weight[k] = 1.0;
1337       nonul[(Int_t)stat[3]] = fEntriesCH[k];
1338       stat[3]++;
1339     }
1340     else weight[k] = 0.0;
1341   }
1342   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight); 
1343   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight); 
1344   stat[2]          = TMath::RMS((Int_t)stat[3],nonul); 
1345
1346   return stat;
1347
1348 }
1349 //____________________________________________________________________________
1350 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1351 {
1352   //
1353   // Return a Int_t[4] with:
1354   // 0 Mean number of entries
1355   // 1 median of number of entries
1356   // 2 rms of number of entries
1357   // 3 number of group with entries
1358   //
1359
1360   Double_t *stat      = new Double_t[4]; 
1361   stat[3]             = 0.0;
1362
1363   Int_t    nbofgroups = 540;
1364   Double_t *weight    = new Double_t[nbofgroups];
1365   Int_t    *nonul     = new Int_t[nbofgroups]; 
1366
1367   for(Int_t k = 0; k < nbofgroups; k++){
1368     if(fEntriesLinearFitter[k] > 0) {
1369       weight[k] = 1.0;
1370       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1371       stat[3]++;     
1372     }
1373     else weight[k] = 0.0;
1374   }
1375   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight); 
1376   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight); 
1377   stat[2]          = TMath::RMS((Int_t)stat[3],nonul); 
1378
1379   return stat;
1380
1381 }
1382 //_____________________________________________________________________________
1383 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1384 {
1385   //
1386   // Should be between 0 and 6
1387   //
1388  
1389   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1390     AliInfo("The number of groups must be between 0 and 6!");
1391   } 
1392   else {
1393     fNgroupprf = numberGroupsPRF;
1394   }
1395
1396
1397 //_____________________________________________________________________________
1398 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1399 {
1400   //
1401   // Set the factor that will divide the deposited charge
1402   // to fit in the histo range [0,300]
1403   //
1404  
1405   if (RelativeScale > 0.0) {
1406     fRelativeScale = RelativeScale;
1407   } 
1408   else {
1409     AliInfo("RelativeScale must be strict positif!");
1410   }
1411
1412
1413
1414 //_____________________________________________________________________________
1415 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1416 {
1417   //
1418   // Set the mode of calibration group in the z direction for the parameter i
1419   // 
1420
1421   if ((Nz >= 0) && 
1422       (Nz <  5)) {
1423     fCalibraMode->SetNz(i, Nz); 
1424   }
1425   else { 
1426     AliInfo("You have to choose between 0 and 4");
1427   }
1428
1429 }
1430
1431 //_____________________________________________________________________________
1432 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1433 {
1434   //
1435   // Set the mode of calibration group in the rphi direction for the parameter i
1436   //
1437  
1438   if ((Nrphi >= 0) && 
1439       (Nrphi <  7)) {
1440     fCalibraMode->SetNrphi(i ,Nrphi); 
1441   }
1442   else {
1443     AliInfo("You have to choose between 0 and 6");
1444   }
1445
1446 }
1447 //____________Protected Functions______________________________________________
1448 //____________Create the 2D histo to be filled online__________________________
1449 //
1450 //_____________________________________________________________________________
1451 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1452 {
1453   //
1454   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1455   // If fNgroupprf is zero then no binning in tnp
1456   //
1457
1458   TString name("Nz");
1459   name += fCalibraMode->GetNz(2);
1460   name += "Nrphi";
1461   name += fCalibraMode->GetNrphi(2);
1462   name += "Ngp";
1463   name += fNgroupprf;
1464
1465   if(fNgroupprf != 0){
1466     
1467     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1468                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1469     fPRF2d->SetYTitle("Det/pad groups");
1470     fPRF2d->SetXTitle("Position x/W [pad width units]");
1471     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1472     fPRF2d->SetStats(0);
1473   }
1474   else{
1475     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1476                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1477     fPRF2d->SetYTitle("Det/pad groups");
1478     fPRF2d->SetXTitle("Position x/W [pad width units]");
1479     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1480     fPRF2d->SetStats(0);
1481   }
1482
1483 }
1484
1485 //_____________________________________________________________________________
1486 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1487 {
1488   //
1489   // Create the 2D histos
1490   //
1491
1492   TString name("Nz");
1493   name += fCalibraMode->GetNz(1);
1494   name += "Nrphi";
1495   name += fCalibraMode->GetNrphi(1);
1496   
1497   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1498                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1499                          ,nn,0,nn);
1500   fPH2d->SetYTitle("Det/pad groups");
1501   fPH2d->SetXTitle("time [#mus]");
1502   fPH2d->SetZTitle("<PH> [a.u.]");
1503   fPH2d->SetStats(0);
1504
1505 }
1506
1507 //_____________________________________________________________________________
1508 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1509 {
1510   //
1511   // Create the 2D histos
1512   //
1513
1514   TString name("Nz");
1515   name += fCalibraMode->GetNz(0);
1516   name += "Nrphi";
1517   name += fCalibraMode->GetNrphi(0);
1518
1519   fCH2d = new TH2I("CH2d",(const Char_t *) name
1520                    ,fNumberBinCharge,0,300,nn,0,nn);
1521   fCH2d->SetYTitle("Det/pad groups");
1522   fCH2d->SetXTitle("charge deposit [a.u]");
1523   fCH2d->SetZTitle("counts");
1524   fCH2d->SetStats(0);
1525   fCH2d->Sumw2();
1526
1527 }
1528 //____________Offine tracking in the AliTRDtracker_____________________________
1529 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1530 {
1531   //
1532   // For the offline tracking or mcm tracklets
1533   // This function will be called in the functions UpdateHistogram... 
1534   // to fill the info of a track for the relativ gain calibration
1535   //
1536         
1537   Int_t nb            =  0;   // Nombre de zones traversees
1538   Int_t fd            = -1;   // Premiere zone non nulle
1539   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1540   
1541   
1542   // See if the track goes through different zones
1543   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1544     if (fAmpTotal[k] > 0.0) {
1545       totalcharge += fAmpTotal[k];
1546       nb++;
1547       if (nb == 1) {
1548         fd = k;
1549       }
1550     }
1551   }
1552
1553  
1554   switch (nb)
1555     { 
1556     case 1:
1557       fNumberUsedCh[0]++;
1558       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1559       if (fHisto2d) {
1560         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1561         //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1562       }
1563       if (fVector2d) {
1564         fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1565       }
1566       break;
1567     case 2:
1568       if ((fAmpTotal[fd]   > 0.0) && 
1569           (fAmpTotal[fd+1] > 0.0)) {
1570         // One of the two very big
1571         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1572           if (fHisto2d) {
1573             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1574             //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1575           }
1576           if (fVector2d) {
1577             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1578           }
1579           fNumberUsedCh[1]++;
1580           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1581         }
1582         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1583           if (fHisto2d) {
1584             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1585             //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1586           }
1587           if (fVector2d) {
1588             fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1589           }
1590           fNumberUsedCh[1]++;
1591           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1592         }
1593       }
1594       if (fCalibraMode->GetNfragZ(0) > 1) {
1595         if (fAmpTotal[fd] > 0.0) {
1596           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1597             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1598               // One of the two very big
1599               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1600                 if (fHisto2d) {
1601                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1602                   //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1603                 }
1604                 if (fVector2d) {
1605                   fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1606                 }
1607                 fNumberUsedCh[1]++;
1608                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1609               }
1610               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1611                 if (fHisto2d) {
1612                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1613                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1614                 }
1615                 fNumberUsedCh[1]++;
1616                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1617                 if (fVector2d) {
1618                   fCalibraVector->UpdateVectorCH(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1619                 }
1620               }
1621             }
1622           }
1623         }
1624       }
1625       break;
1626     default: break;
1627     }
1628 }
1629 //____________Offine tracking in the AliTRDtracker_____________________________
1630 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1631 {
1632   //
1633   // For the offline tracking or mcm tracklets
1634   // This function will be called in the functions UpdateHistogram... 
1635   // to fill the info of a track for the drift velocity  calibration
1636   //
1637     
1638   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1639   Int_t fd1 = -1; // Premiere zone non nulle
1640   Int_t fd2 = -1; // Deuxieme zone non nulle
1641   Int_t k1  = -1; // Debut de la premiere zone
1642   Int_t k2  = -1; // Debut de la seconde zone
1643
1644   // See if the track goes through different zones
1645   for (Int_t k = 0; k < fTimeMax; k++) {
1646     if (fPHValue[k] > 0.0) {
1647       if (fd1 == -1) {
1648         fd1 = fPHPlace[k];
1649         k1  = k;              
1650       }
1651       if (fPHPlace[k] != fd1) {
1652         if (fd2 == -1) {
1653           k2  = k;
1654           fd2 = fPHPlace[k];
1655           nb  = 2;
1656         }
1657         if (fPHPlace[k] != fd2) {
1658           nb = 3;
1659         }
1660       }
1661     }
1662   }
1663
1664   switch(nb)
1665     {
1666     case 1:
1667       fNumberUsedPh[0]++;
1668       for (Int_t i = 0; i < fTimeMax; i++) {
1669         if (fHisto2d) {
1670           fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1671         }
1672         if (fVector2d) {
1673           fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1674         }
1675       }
1676       break;
1677     case 2:
1678       if ((fd1 == fd2+1) || 
1679           (fd2 == fd1+1)) {
1680         // One of the two fast all the think
1681         if (k2 > (k1+fDifference)) {
1682           //we choose to fill the fd1 with all the values
1683           fNumberUsedPh[1]++;
1684           for (Int_t i = 0; i < fTimeMax; i++) {
1685             if (fHisto2d) {
1686               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1687             }
1688             if (fVector2d) {
1689               fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1690             }
1691           }
1692         }
1693         if ((k2+fDifference) < fTimeMax) {
1694           //we choose to fill the fd2 with all the values
1695           fNumberUsedPh[1]++;
1696           for (Int_t i = 0; i < fTimeMax; i++) {
1697             if (fHisto2d) {
1698               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1699             }
1700           if (fVector2d) {
1701             fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1702           }
1703           }
1704         }
1705       }
1706       // Two zones voisines sinon rien!
1707       if (fCalibraMode->GetNfragZ(1) > 1) {
1708         // Case 2
1709         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1710           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1711             // One of the two fast all the think
1712             if (k2 > (k1+fDifference)) {
1713               //we choose to fill the fd1 with all the values
1714               fNumberUsedPh[1]++;
1715               for (Int_t i = 0; i < fTimeMax; i++) {
1716                 if (fHisto2d) {
1717                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1718                 }
1719                 if (fVector2d) {
1720                   fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1721                 }
1722               }
1723             }
1724             if ((k2+fDifference) < fTimeMax) {
1725               //we choose to fill the fd2 with all the values
1726               fNumberUsedPh[1]++;
1727               for (Int_t i = 0; i < fTimeMax; i++) {
1728                 if (fHisto2d) {
1729                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1730                 }
1731                 if (fVector2d) {
1732                   fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1733                 }
1734               }
1735             }
1736           }
1737         }
1738         // Two zones voisines sinon rien!
1739         // Case 3
1740         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1741           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1742             // One of the two fast all the think
1743             if (k2 > (k1 + fDifference)) {
1744               //we choose to fill the fd1 with all the values
1745               fNumberUsedPh[1]++;
1746               for (Int_t i = 0; i < fTimeMax; i++) {
1747                 if (fHisto2d) {
1748                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1749                 }
1750                 if (fVector2d) {
1751                   fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd1),i,fPHValue[i]);
1752                 }
1753               }
1754             }
1755             if ((k2+fDifference) < fTimeMax) {
1756               //we choose to fill the fd2 with all the values
1757               fNumberUsedPh[1]++;
1758               for (Int_t i = 0; i < fTimeMax; i++) {
1759                 if (fHisto2d) {
1760                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1761                 }
1762                 if (fVector2d) {
1763                   fCalibraVector->UpdateVectorPH((fCalibraMode->GetXbins(1)+fd2),i,fPHValue[i]);
1764                 }
1765               }
1766             }
1767           }
1768         }
1769       }
1770       break;
1771     default: break;
1772     } 
1773 }
1774 //____________Offine tracking in the AliTRDtracker_____________________________
1775 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
1776 {
1777   //
1778   // For the offline tracking
1779   // This function will be called in the functions UpdateHistogram... 
1780   // to fill the find the parameter P1 of a track for the drift velocity  calibration
1781   //
1782   
1783   //Number of points: if less than 3 return kFALSE
1784   Int_t Npoints = fListClusters->GetEntriesFast();
1785   if(Npoints <= 2) return kFALSE;
1786
1787   //Variables
1788   TLinearFitter linearFitterTracklet      = TLinearFitter(2,"pol1");        // TLinearFitter per tracklet
1789   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
1790   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
1791   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
1792   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
1793   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
1794   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!  
1795   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
1796   Double_t pointError                 = 0.0;                                // error after straight line fit 
1797   Int_t    detector                   = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
1798   Int_t    snpright                   = 1;                                  // if we took in the middle snp
1799   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
1800   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
1801   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
1802   Int_t     nbli                      = 0;                                  // number linear fitter points
1803   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
1804
1805   linearFitterTracklet.StoreData(kFALSE);
1806   linearFitterTracklet.ClearPoints();
1807   
1808   //if more than one row
1809   Int_t    rowp                       = -1;                              // if it crosses a pad row
1810
1811   //tiltingangle
1812   tiltingangle                        = padplane->GetTiltingAngle();
1813   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
1814
1815   //Fill with points
1816   for(Int_t k = 0; k < Npoints; k++){
1817     
1818     AliTRDcluster *cl                 = (AliTRDcluster *) fListClusters->At(k);
1819     Double_t ycluster                 = cl->GetY();
1820     Int_t time                        = cl->GetLocalTimeBin();
1821     Double_t timeis                   = time/fSf;
1822     //See if cross two pad rows
1823     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
1824     if(k==0) rowp                     = row;
1825     if(row != rowp) crossrow          = 1;
1826     //Take in the middle of the chamber
1827     if(time > (Int_t) 10) {
1828       z   = cl->GetZ();
1829       y   = cl->GetY();  
1830       snp = fPar2[time];
1831       tgl = fPar3[time];
1832     }
1833     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1834     nbli++;
1835   }
1836   if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() <= 10) snpright = 0;
1837   if(nbli <= 2) return kFALSE; 
1838   
1839   // Do the straight line fit now
1840   TVectorD pars;
1841   linearFitterTracklet.Eval();
1842   linearFitterTracklet.GetParameters(pars);
1843   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
1844   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
1845   dydt  = pars[1]; 
1846  
1847   if( TMath::Abs(snp) <  1.){
1848     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1849   } 
1850   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1851
1852   if(fDebugLevel > 0){
1853     if ( !fDebugStreamer ) {
1854       //debug stream
1855       TDirectory *backup = gDirectory;
1856       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
1857       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1858     } 
1859     
1860     (* fDebugStreamer) << "VDRIFT0"<<
1861       "Npoints="<<Npoints<<
1862       "\n"; 
1863   
1864     
1865     (* fDebugStreamer) << "VDRIFT"<<
1866       "snpright="<<snpright<<
1867       "Npoints="<<Npoints<<
1868       "nbli="<<nbli<<
1869       "detector="<<detector<<
1870       "snp="<<snp<<
1871       "tnp="<<tnp<<
1872       "tgl="<<tgl<<
1873       "tnt="<<tnt<<
1874       "y="<<y<<
1875       "z="<<z<<
1876       "dydt="<<dydt<<
1877       "dzdx="<<dzdx<<
1878       "crossrow="<<crossrow<<
1879       "errorpar="<<errorpar<<
1880       "pointError="<<pointError<<
1881       "\n";     
1882
1883   }
1884   
1885   if(Npoints < fNumberClusters) return kFALSE;
1886   if(snpright == 0) return kFALSE;
1887   if(pointError >= 0.1) return kFALSE;
1888   if(crossrow == 1) return kFALSE;
1889   
1890   if(fLinearFitterOn){
1891     //Add to the linear fitter of the detector
1892     if( TMath::Abs(snp) <  1.){
1893       Double_t x = tnp-dzdx*tnt; 
1894       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
1895       if(fLinearFitterDebugOn) {
1896         ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars[1]);
1897       }
1898       fEntriesLinearFitter[detector]++;
1899     }
1900   }
1901   //AliInfo("End of FindP1TrackPH with success!")
1902   return kTRUE;
1903
1904 }
1905 //____________Offine tracking in the AliTRDtracker_____________________________
1906 Bool_t AliTRDCalibraFillHisto::HandlePRF()
1907 {
1908   //
1909   // For the offline tracking
1910   // Fit the tracklet with a line and take the position as reference for the PRF
1911   //
1912   
1913   //Number of points
1914   Int_t Npoints  = fListClusters->GetEntriesFast();                         // number of total points
1915   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit 
1916   Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
1917  
1918
1919   // To see the difference due to the fit
1920   Double_t *padPositions;
1921   padPositions = new Double_t[Npoints];
1922   for(Int_t k = 0; k < Npoints; k++){
1923     padPositions[k] = 0.0;
1924   } 
1925
1926
1927   //Find the position by a fit
1928   TLinearFitter fitter(2,"pol1");
1929   fitter.StoreData(kFALSE);
1930   fitter.ClearPoints();
1931   for(Int_t k = 0;  k < Npoints; k++){
1932     //Take the cluster
1933     AliTRDcluster *cl  = (AliTRDcluster *) fListClusters->At(k);
1934     Short_t  *signals  = cl->GetSignals();
1935     Double_t     time  = cl->GetLocalTimeBin();
1936     //Calculate x if possible 
1937     Float_t xcenter    = 0.0;    
1938     Bool_t  echec      = kTRUE;   
1939     if((time<=7) || (time>=21)) continue; 
1940     // Center 3 balanced: position with the center of the pad
1941     if ((((Float_t) signals[3]) > 0.0) && 
1942         (((Float_t) signals[2]) > 0.0) && 
1943         (((Float_t) signals[4]) > 0.0)) {
1944       echec = kFALSE;
1945       // Security if the denomiateur is 0 
1946       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / 
1947            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1948         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1949           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) 
1950                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1951       }
1952       else {
1953         echec = kTRUE;
1954       }
1955     }
1956     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1957     if(echec) continue;
1958     //if no echec: calculate with the position of the pad
1959     // Position of the cluster
1960     Double_t       padPosition = xcenter +  cl->GetPad();
1961     padPositions[k]            = padPosition;
1962     Nb3pc++;
1963     fitter.AddPoint(&time, padPosition,1);
1964   }//clusters loop
1965
1966   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
1967   if(Nb3pc < 3) return kFALSE;
1968   fitter.Eval();
1969   TVectorD line(2);
1970   fitter.GetParameters(line);
1971   Float_t  pointError  = -1.0;
1972   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
1973   
1974
1975   // Now fill the PRF  
1976   for(Int_t k = 0;  k < Npoints; k++){
1977     //Take the cluster
1978     AliTRDcluster *cl      = (AliTRDcluster *) fListClusters->At(k);
1979     Short_t  *signals      = cl->GetSignals();              // signal
1980     Double_t     time      = cl->GetLocalTimeBin();         // time bin
1981     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1982     Float_t padPos         = cl->GetPad();                  // middle pad
1983     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit 
1984     Float_t ycenter        = 0.0;                           // relative center charge
1985     Float_t ymin           = 0.0;                           // relative left charge
1986     Float_t ymax           = 0.0;                           // relative right charge
1987     Double_t tgl           = fPar3[(Int_t)time];            // dz/dl and not dz/dx
1988     Double_t pt            = fPar4[(Int_t)time];            // pt
1989     Float_t  dzdx          = 0.0;                           // dzdx
1990
1991
1992     //Requiere simply two pads clusters at least
1993     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1994        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1995       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1996       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1997       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1998       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum; 
1999     }
2000     
2001     //calibration group
2002     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2003     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2004     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponidng group
2005     Double_t  snp          = fPar2[(Int_t)time];                        // sin angle in xy plan
2006     Float_t   xcl          = cl->GetY();                                // y cluster
2007     Float_t   qcl          = cl->GetQ();                                // charge cluster 
2008     Int_t     plane        = GetPlane(detector);                        // plane 
2009     Int_t     chamber      = GetChamber(detector);                      // chamber  
2010     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2011     Double_t  x            = dpad;                                      // reconstructed position moved
2012     Float_t   Ep           = pointError;                                // error of fit
2013     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2014     Float_t   signal3      = (Float_t)signals[3];                       // signal
2015     Float_t   signal2      = (Float_t)signals[2];                       // signal
2016     Float_t   signal4      = (Float_t)signals[4];                       // signal
2017     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2018     Float_t tnp            = 0.0;
2019     if(TMath::Abs(snp) < 1.0){
2020       tnp = snp / (TMath::Sqrt(1-snp*snp));
2021       dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2022     }
2023    
2024
2025     if(fDebugLevel > 0){
2026       if ( !fDebugStreamer ) {
2027         //debug stream
2028         TDirectory *backup = gDirectory;
2029         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibra.root");
2030         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2031       }     
2032       
2033       (* fDebugStreamer) << "PRF0"<<
2034         "caligroup="<<caligroup<<
2035         "detector="<<detector<<
2036         "plane="<<plane<<
2037         "chamber="<<chamber<<
2038         "Npoints="<<Npoints<<
2039         "Np="<<Nb3pc<<
2040         "Ep="<<Ep<<
2041         "snp="<<snp<<
2042         "tnp="<<tnp<<    
2043         "tgl="<<tgl<<  
2044         "dzdx="<<dzdx<<  
2045         "pt="<<pt<<    
2046         "padPos="<<padPos<<
2047         "padPosition="<<padPositions[k]<<
2048         "padPosTracklet="<<padPosTracklet<<
2049         "x="<<x<<
2050         "ycenter="<<ycenter<<
2051         "ymin="<<ymin<<
2052         "ymax="<<ymax<<
2053         "xcl="<<xcl<<
2054         "qcl="<<qcl<< 
2055         "signal1="<<signal1<<    
2056         "signal2="<<signal2<<
2057         "signal3="<<signal3<<
2058         "signal4="<<signal4<<
2059         "signal5="<<signal5<<
2060         "time="<<time<<
2061         "\n";     
2062       x = xdiff;
2063       Int_t type=0;
2064       Float_t y = ycenter;
2065       (* fDebugStreamer) << "PRFALL"<<
2066         "caligroup="<<caligroup<<
2067         "detector="<<detector<<
2068         "plane="<<plane<<
2069         "chamber="<<chamber<<
2070         "Npoints="<<Npoints<<
2071         "Np="<<Nb3pc<<
2072         "Ep="<<Ep<<
2073         "type="<<type<<
2074         "snp="<<snp<<
2075         "tnp="<<tnp<<
2076         "tgl="<<tgl<<  
2077         "dzdx="<<dzdx<< 
2078         "pt="<<pt<< 
2079         "padPos="<<padPos<<
2080         "padPosition="<<padPositions[k]<<
2081         "padPosTracklet="<<padPosTracklet<<
2082         "x="<<x<<
2083         "y="<<y<<           
2084         "xcl="<<xcl<<
2085         "qcl="<<qcl<<
2086         "signal1="<<signal1<<
2087         "signal2="<<signal2<<
2088         "signal3="<<signal3<<
2089         "signal4="<<signal4<<
2090         "signal5="<<signal5<<
2091         "time="<<time<<
2092         "\n";
2093       x=-(xdiff+1);
2094       y = ymin;
2095       type=-1;
2096       (* fDebugStreamer) << "PRFALL"<<
2097         "caligroup="<<caligroup<<
2098         "detector="<<detector<<
2099         "plane="<<plane<<
2100         "chamber="<<chamber<<
2101         "Npoints="<<Npoints<<
2102         "Np="<<Nb3pc<<
2103         "Ep="<<Ep<<
2104         "type="<<type<<
2105         "snp="<<snp<<
2106         "tnp="<<tnp<<
2107         "tgl="<<tgl<<  
2108         "dzdx="<<dzdx<< 
2109         "pt="<<pt<< 
2110         "padPos="<<padPos<<
2111         "padPosition="<<padPositions[k]<<
2112         "padPosTracklet="<<padPosTracklet<<
2113         "x="<<x<<
2114         "y="<<y<<
2115         "xcl="<<xcl<<
2116         "qcl="<<qcl<<
2117         "signal1="<<signal1<<
2118         "signal2="<<signal2<<
2119         "signal3="<<signal3<<
2120         "signal4="<<signal4<<
2121         "signal5="<<signal5<<
2122         "time="<<time<<
2123         "\n";
2124       x=1-xdiff;
2125       y = ymax;
2126       type=1;
2127       
2128       (* fDebugStreamer) << "PRFALL"<<
2129         "caligroup="<<caligroup<<
2130         "detector="<<detector<<
2131         "plane="<<plane<<
2132         "chamber="<<chamber<<
2133         "Npoints="<<Npoints<<
2134         "Np="<<Nb3pc<<
2135         "Ep="<<Ep<<
2136         "type="<<type<<
2137         "snp="<<snp<<
2138         "tnp="<<tnp<<   
2139         "tgl="<<tgl<<  
2140         "dzdx="<<dzdx<< 
2141         "pt="<<pt<< 
2142         "padPos="<<padPos<<
2143         "padPosition="<<padPositions[k]<<
2144         "padPosTracklet="<<padPosTracklet<<
2145         "x="<<x<<
2146         "y="<<y<<
2147         "xcl="<<xcl<<
2148         "qcl="<<qcl<<
2149         "signal1="<<signal1<<
2150         "signal2="<<signal2<<
2151         "signal3="<<signal3<<
2152         "signal4="<<signal4<<
2153         "signal5="<<signal5<<
2154         "time="<<time<<
2155         "\n";
2156       
2157     }
2158     
2159     // some cuts
2160     if(Npoints < fNumberClusters) continue;
2161     if(Nb3pc <= 5) continue;
2162     if((time >= 21) || (time < 7)) continue;
2163     if(TMath::Abs(snp) >= 1.0) continue;
2164     if(qcl < 80) continue; 
2165     
2166     Bool_t echec   = kFALSE;
2167     Double_t shift = 0.0;
2168     //Calculate the shift in x coresponding to this tnp
2169     if(fNgroupprf != 0.0){
2170       shift      = -3.0*(fNgroupprf-1)-1.5;
2171       Double_t limithigh  = -0.2*(fNgroupprf-1);
2172       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2173       else{
2174         while(tnp > limithigh){
2175           limithigh += 0.2;
2176           shift += 3.0;
2177         }
2178       }
2179     }
2180     if (fHisto2d && !echec) {
2181       fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2182       if(ymin > 0.0) {
2183         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2184         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2185       }
2186       if(ymax > 0.0) {
2187         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2188         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2189       }
2190     }
2191     //Not equivalent anymore here!
2192     if (fVector2d && !echec) {
2193       fCalibraVector->UpdateVectorPRF(caligroup,shift+dpad,ycenter);
2194       if(ymin > 0.0) {
2195         fCalibraVector->UpdateVectorPRF(caligroup,shift-(dpad+1.0),ymin);
2196         fCalibraVector->UpdateVectorPRF(caligroup,shift+(dpad+1.0),ymin);
2197       }
2198       if(ymax > 0.0) {
2199         fCalibraVector->UpdateVectorPRF(caligroup,shift+1.0-dpad,ymax);
2200         fCalibraVector->UpdateVectorPRF(caligroup,shift-1.0+dpad,ymax);
2201       }
2202     }
2203   }
2204   return kTRUE;
2205   
2206 }
2207 //
2208 //____________Some basic geometry function_____________________________________
2209 //
2210 //_____________________________________________________________________________
2211 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2212 {
2213   //
2214   // Reconstruct the plane number from the detector number
2215   //
2216
2217   return ((Int_t) (d % 6));
2218
2219 }
2220
2221 //_____________________________________________________________________________
2222 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2223 {
2224   //
2225   // Reconstruct the chamber number from the detector number
2226   //
2227   Int_t fgkNplan = 6;
2228
2229   return ((Int_t) (d % 30) / fgkNplan);
2230
2231 }
2232
2233 //_____________________________________________________________________________
2234 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2235 {
2236   //
2237   // Reconstruct the sector number from the detector number
2238   //
2239   Int_t fg = 30;
2240
2241   return ((Int_t) (d / fg));
2242
2243 }
2244 //_____________________________________________________________________
2245 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency, Bool_t force)
2246 {
2247     //
2248     // return pointer to fPH2d TProfile2D
2249     // if force is true create a new TProfile2D if it doesn't exist allready
2250     //
2251     if ( !force || fPH2d )
2252         return fPH2d;
2253
2254     // Some parameters
2255     fTimeMax = nbtimebin;
2256     fSf      = samplefrequency;
2257   
2258     /*
2259       AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
2260       ,fCalibraMode->GetNz(1)
2261       ,fCalibraMode->GetNrphi(1)));
2262       
2263       // Calcul the number of Xbins
2264       Int_t Ntotal1 = 0;
2265       fCalibraMode->ModePadCalibration(2,1);
2266       fCalibraMode->ModePadFragmentation(0,2,0,1);
2267       fCalibraMode->SetDetChamb2(1);
2268       Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
2269       fCalibraMode->ModePadCalibration(0,1);
2270       fCalibraMode->ModePadFragmentation(0,0,0,1);
2271       fCalibraMode->SetDetChamb0(1);
2272       Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
2273       AliInfo(Form("Total number of Xbins: %d",Ntotal1));
2274       
2275       CreatePH2d(Ntotal1);
2276     */  
2277
2278     CreatePH2d(540);
2279
2280     return fPH2d;
2281 }
2282 //_____________________________________________________________________
2283 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2284 {
2285     //
2286     // return pointer to TLinearFitter Calibration
2287     // if force is true create a new TLinearFitter if it doesn't exist allready
2288     //
2289     if ( !force || fLinearFitterArray.UncheckedAt(detector) )
2290         return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2291
2292     // if we are forced and TLinearFitter doesn't yet exist create it
2293
2294     // new TLinearFitter
2295     TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2296     fLinearFitterArray.AddAt(linearfitter,detector);
2297     return linearfitter;
2298 }
2299 //_____________________________________________________________________
2300 TH2F* AliTRDCalibraFillHisto::GetLinearFitterHisto(Int_t detector, Bool_t force)
2301 {
2302     //
2303     // return pointer to TH2F histo 
2304     // if force is true create a new histo if it doesn't exist allready
2305     //
2306     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
2307         return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
2308
2309     // if we are forced and TLinearFitter doesn't yes exist create it
2310
2311     // new TH2F
2312     TString name("LFDV");
2313     name += detector;
2314     
2315     TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
2316                           ,100,-1.0,1.0,100
2317                           ,-2.0,2.0);
2318     lfdv->SetXTitle("tan(phi_{track})");
2319     lfdv->SetYTitle("dy/dt");
2320     lfdv->SetZTitle("Number of clusters");
2321     lfdv->SetStats(0);
2322     lfdv->SetDirectory(0);
2323
2324     fLinearFitterHistoArray.AddAt(lfdv,detector);
2325     return lfdv;
2326 }
2327 //_____________________________________________________________________
2328 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2329 {
2330   //
2331   // FillCH2d: Marian style
2332   // 
2333   
2334   //skip simply the value out of range
2335   if((y>=300.0) || (y<0.0)) return;
2336   
2337   //Calcul the y place
2338   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2339   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2340   
2341   //Fill
2342   fCH2d->GetArray()[place]++;
2343
2344 }