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