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