Modifications in the calibration/preprocessor for new raw data version
[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 "AliTRDRawStreamV2.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(AliTRDRawStreamV2 *rawStream, Bool_t nocheck)
1208 {
1209   //
1210   // Event Processing loop - AliTRDRawStreamV2
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       //printf("detector %d, nbtimebin %d, iTimeBin %d\n",fDetectorPreviousTrack,nbtimebin,iTimeBin);
1257       
1258       Int_t fin     = TMath::Min(fTimeMax,(iTimeBin+3));
1259       Int_t n       = 0;
1260       for(Int_t itime = iTimeBin; itime < fin; itime++){
1261         // should extract baseline here!
1262         if(signal[n]>13) phvalue[itime] = signal[n];
1263         //printf("signal is %d for %d\n",signal[n],n);
1264         n++;
1265       }
1266     }
1267   
1268     // fill the last one
1269     if(fDetectorPreviousTrack != -1){
1270       if(TMath::Mean(fTimeMax,phvalue)>20.0){
1271         withInput = 2;
1272         for(Int_t k = 0; k < fTimeMax; k++){
1273           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax);
1274           phvalue[k] = 10;
1275           //row[k]     = -1;
1276           //col[k]     = -1;
1277         }
1278       }
1279     }
1280     
1281   }
1282   else{
1283
1284     while (rawStream->Next()) {
1285
1286       Int_t idetector = rawStream->GetDet();                            //  current detector
1287       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
1288         if(TMath::Mean(nbtimebin,phvalue)>20.0){
1289           withInput = 2;
1290           for(Int_t k = 0; k < nbtimebin; k++){
1291             UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1292             phvalue[k] = 10;
1293             //row[k]     = -1;
1294             //col[k]     = -1;
1295           }
1296         }
1297       }
1298       fDetectorPreviousTrack = idetector;
1299       nbtimebin         = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
1300       Int_t iTimeBin    = rawStream->GetTimeBin();                       //  current time bin
1301       //row[iTimeBin]   = rawStream->GetRow();                           //  current row
1302       //col[iTimeBin]   = rawStream->GetCol();                           //  current col     
1303       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
1304       //printf("detector %d, nbtimebin %d, iTimeBin %d\n",fDetectorPreviousTrack,nbtimebin,iTimeBin);
1305
1306       Int_t fin     = TMath::Min(nbtimebin,(iTimeBin+3));
1307       Int_t n       = 0;
1308       for(Int_t itime = iTimeBin; itime < fin; itime++){
1309         // should extract baseline here!
1310         if(signal[n]>13) phvalue[itime] = signal[n];
1311         //printf("signal is %d for %d\n",signal[n],n);
1312         n++;
1313       }
1314     }
1315     
1316     // fill the last one
1317     if(fDetectorPreviousTrack != -1){
1318       if(TMath::Mean(nbtimebin,phvalue)>20.0){
1319         withInput = 2;
1320         for(Int_t k = 0; k < nbtimebin; k++){
1321           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1322           phvalue[k] = 10;
1323           //row[k]     = -1;
1324           //col[k]     = -1;
1325         }
1326       }
1327     }
1328   }
1329   
1330   return withInput;
1331
1332 }
1333 //_____________________________________________________________________
1334 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
1335 {
1336   //
1337   //  Event processing loop - AliRawReader
1338   //
1339
1340
1341   AliTRDRawStreamV2 rawStream(rawReader);
1342
1343   rawReader->Select("TRD");
1344
1345   return ProcessEventDAQ(&rawStream, nocheck);
1346 }
1347 //_________________________________________________________________________
1348 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
1349 #ifdef ALI_DATE
1350                                                eventHeaderStruct *event,
1351                                                Bool_t nocheck
1352 #else
1353                                                eventHeaderStruct* /*event*/,
1354                                                Bool_t /*nocheck*/
1355             
1356 #endif 
1357                                    )
1358 {
1359   //
1360   //  process date event
1361   //
1362 #ifdef ALI_DATE
1363     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1364     Int_t result=ProcessEventDAQ(rawReader, nocheck);
1365     delete rawReader;
1366     return result;
1367 #else
1368     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1369     return 0;
1370 #endif
1371
1372 }
1373 //____________Online trackling in AliTRDtrigger________________________________
1374 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1375 {
1376   //
1377   // For the DAQ
1378   // Fill a simple average pulse height
1379   //
1380   
1381   // Localisation of the Xbins involved
1382   //LocalisationDetectorXbins(det);
1383
1384   // Row  and position in the pad groups
1385   //Int_t posr = 0;
1386   //if (fCalibraMode->GetNnZ(1) != 0) {
1387   //  posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1388   //}
1389  
1390   // Col of cluster and position in the pad groups
1391   //Int_t posc = 0;
1392   //if (fCalibraMode->GetNnRphi(1) != 0) {
1393   //  posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1394   //}
1395   
1396   //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);   
1397   
1398   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1399   
1400   return kTRUE;
1401   
1402 }
1403 //____________Write_____________________________________________________
1404 //_____________________________________________________________________
1405 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1406 {
1407   //
1408   //  Write infos to file
1409   //
1410   
1411   //For debugging
1412   if ( fDebugStreamer ) {
1413     delete fDebugStreamer;
1414     fDebugStreamer = 0x0;
1415   }
1416
1417   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1418                ,fNumberTrack
1419                ,fNumberUsedCh[0]
1420                ,fNumberUsedCh[1]
1421                ,fNumberUsedPh[0]
1422                ,fNumberUsedPh[1]));
1423   
1424   TDirectory *backup = gDirectory;
1425   TString option;
1426   
1427   if ( append )
1428     option = "update";
1429   else
1430     option = "recreate";
1431   
1432   TFile f(filename,option.Data());
1433   
1434   TStopwatch stopwatch;
1435   stopwatch.Start();
1436   if(fVector2d) {
1437     f.WriteTObject(fCalibraVector);
1438   }
1439
1440   if (fCH2dOn ) {
1441     if (fHisto2d) {
1442       f.WriteTObject(fCH2d);
1443     }
1444   }
1445   if (fPH2dOn ) {
1446     if (fHisto2d) {
1447       f.WriteTObject(fPH2d);
1448     }
1449   }
1450   if (fPRF2dOn) {
1451     if (fHisto2d) {
1452         f.WriteTObject(fPRF2d);
1453     }
1454   }
1455   if(fLinearFitterOn){
1456     AnalyseLinearFitter();
1457     f.WriteTObject(fLinearVdriftFit);
1458   }
1459    
1460   f.Close();
1461   
1462   if ( backup ) backup->cd();
1463   
1464   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1465                ,stopwatch.RealTime(),stopwatch.CpuTime()));
1466 }
1467 //___________________________________________probe the histos__________________________________________________
1468 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1469 {
1470   //
1471   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1472   // debug mode with 2 for TH2I and 3 for TProfile2D
1473   // It gives a pointer to a Double_t[7] with the info following...
1474   // [0] : number of calibration groups with entries
1475   // [1] : minimal number of entries found
1476   // [2] : calibration group number of the min
1477   // [3] : maximal number of entries found
1478   // [4] : calibration group number of the max
1479   // [5] : mean number of entries found
1480   // [6] : mean relative error
1481   //
1482
1483   Double_t *info = new Double_t[7];
1484    
1485   // Number of Xbins (detectors or groups of pads)
1486   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
1487   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
1488
1489   // Initialise
1490   Double_t nbwe = 0; //number of calibration groups with entries
1491   Double_t minentries = 0; //minimal number of entries found
1492   Double_t maxentries = 0; //maximal number of entries found
1493   Double_t placemin = 0; //calibration group number of the min
1494   Double_t placemax = -1; //calibration group number of the max
1495   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1496   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1497
1498   Double_t counter = 0;
1499
1500   //Debug
1501   TH1F *NbEntries = 0x0;//distribution of the number of entries
1502   TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1503   TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1504     
1505   // Beginning of the loop over the calibration groups 
1506   for (Int_t idect = 0; idect < nbins; idect++) {
1507
1508     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1509     projch->SetDirectory(0);
1510     
1511     // Number of entries for this calibration group
1512     Double_t nentries = 0.0;
1513     if((i%2) == 0){
1514       for (Int_t k = 0; k < nxbins; k++) {
1515         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1516       }
1517     }
1518     else{
1519       for (Int_t k = 0; k < nxbins; k++) {
1520         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1521         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1522           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1523           counter++;
1524         } 
1525       }
1526     }
1527
1528     //Debug
1529     if(i > 1){
1530       if((!((Bool_t)NbEntries)) && (nentries > 0)){
1531         NbEntries = new TH1F("Number of entries","Number of entries"
1532                                ,100,(Int_t)nentries/2,nentries*2);
1533         NbEntries->SetDirectory(0);
1534         NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1535                                ,nbins,0,nbins);
1536         NbEntriesPerGroup->SetDirectory(0);
1537         NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1538                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1539         NbEntriesPerSp->SetDirectory(0);
1540       }
1541       if(NbEntries){
1542         if(nentries > 0) NbEntries->Fill(nentries);
1543         NbEntriesPerGroup->Fill(idect+0.5,nentries);
1544         NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1545       }
1546     }
1547
1548     //min amd max
1549     if(nentries > maxentries){
1550       maxentries = nentries;
1551       placemax = idect;
1552     }
1553     if(idect == 0) {
1554       minentries = nentries;
1555     }
1556     if(nentries < minentries){
1557       minentries = nentries;
1558       placemin = idect;
1559     }
1560     //nbwe
1561     if(nentries > 0) {
1562       nbwe++;
1563       meanstats += nentries;
1564     }
1565   }//calibration groups loop
1566   
1567   if(nbwe > 0) meanstats /= nbwe;
1568   if(counter > 0) meanrelativerror /= counter;
1569
1570   AliInfo(Form("There are %f calibration groups with entries",nbwe));
1571   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1572   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1573   AliInfo(Form("The mean number of entries is %f",meanstats));
1574   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1575
1576   info[0] = nbwe;
1577   info[1] = minentries;
1578   info[2] = placemin;
1579   info[3] = maxentries;
1580   info[4] = placemax;
1581   info[5] = meanstats;
1582   info[6] = meanrelativerror;
1583
1584   if(i > 1){
1585     gStyle->SetPalette(1);
1586     gStyle->SetOptStat(1111);
1587     gStyle->SetPadBorderMode(0);
1588     gStyle->SetCanvasColor(10);
1589     gStyle->SetPadLeftMargin(0.13);
1590     gStyle->SetPadRightMargin(0.01);
1591     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1592     stat->Divide(2,1);
1593     stat->cd(1);
1594     NbEntries->Draw("");
1595     stat->cd(2);
1596     NbEntriesPerSp->SetStats(0);
1597     NbEntriesPerSp->Draw("");
1598     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1599     stat1->cd();
1600     NbEntriesPerGroup->SetStats(0);
1601     NbEntriesPerGroup->Draw("");
1602   }
1603
1604   return info;
1605
1606 }
1607 //____________________________________________________________________________
1608 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1609 {
1610   //
1611   // Return a Int_t[4] with:
1612   // 0 Mean number of entries
1613   // 1 median of number of entries
1614   // 2 rms of number of entries
1615   // 3 number of group with entries
1616   //
1617
1618   Double_t *stat      = new Double_t[4]; 
1619   stat[3]             = 0.0;
1620
1621   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
1622   Double_t *weight    = new Double_t[nbofgroups];
1623   Int_t    *nonul     = new Int_t[nbofgroups];
1624  
1625   for(Int_t k = 0; k < nbofgroups; k++){
1626     if(fEntriesCH[k] > 0) {
1627       weight[k] = 1.0;
1628       nonul[(Int_t)stat[3]] = fEntriesCH[k];
1629       stat[3]++;
1630     }
1631     else weight[k] = 0.0;
1632   }
1633   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight); 
1634   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight); 
1635   stat[2]          = TMath::RMS((Int_t)stat[3],nonul); 
1636
1637   return stat;
1638
1639 }
1640 //____________________________________________________________________________
1641 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1642 {
1643   //
1644   // Return a Int_t[4] with:
1645   // 0 Mean number of entries
1646   // 1 median of number of entries
1647   // 2 rms of number of entries
1648   // 3 number of group with entries
1649   //
1650
1651   Double_t *stat      = new Double_t[4]; 
1652   stat[3]             = 0.0;
1653
1654   Int_t    nbofgroups = 540;
1655   Double_t *weight    = new Double_t[nbofgroups];
1656   Int_t    *nonul     = new Int_t[nbofgroups]; 
1657
1658   for(Int_t k = 0; k < nbofgroups; k++){
1659     if(fEntriesLinearFitter[k] > 0) {
1660       weight[k] = 1.0;
1661       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1662       stat[3]++;     
1663     }
1664     else weight[k] = 0.0;
1665   }
1666   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight); 
1667   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight); 
1668   stat[2]          = TMath::RMS((Int_t)stat[3],nonul); 
1669
1670   return stat;
1671
1672 }
1673 //_____________________________________________________________________________
1674 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1675 {
1676   //
1677   // Should be between 0 and 6
1678   //
1679  
1680   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1681     AliInfo("The number of groups must be between 0 and 6!");
1682   } 
1683   else {
1684     fNgroupprf = numberGroupsPRF;
1685   }
1686
1687
1688 //_____________________________________________________________________________
1689 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1690 {
1691   //
1692   // Set the factor that will divide the deposited charge
1693   // to fit in the histo range [0,300]
1694   //
1695  
1696   if (RelativeScale > 0.0) {
1697     fRelativeScale = RelativeScale;
1698   } 
1699   else {
1700     AliInfo("RelativeScale must be strict positif!");
1701   }
1702
1703
1704
1705 //_____________________________________________________________________________
1706 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1707 {
1708   //
1709   // Set the mode of calibration group in the z direction for the parameter i
1710   // 
1711
1712   if ((Nz >= 0) && 
1713       (Nz <  5)) {
1714     fCalibraMode->SetNz(i, Nz); 
1715   }
1716   else { 
1717     AliInfo("You have to choose between 0 and 4");
1718   }
1719
1720 }
1721
1722 //_____________________________________________________________________________
1723 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1724 {
1725   //
1726   // Set the mode of calibration group in the rphi direction for the parameter i
1727   //
1728  
1729   if ((Nrphi >= 0) && 
1730       (Nrphi <  7)) {
1731     fCalibraMode->SetNrphi(i ,Nrphi); 
1732   }
1733   else {
1734     AliInfo("You have to choose between 0 and 6");
1735   }
1736
1737 }
1738 //____________Protected Functions______________________________________________
1739 //____________Create the 2D histo to be filled online__________________________
1740 //
1741 //_____________________________________________________________________________
1742 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1743 {
1744   //
1745   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1746   // If fNgroupprf is zero then no binning in tnp
1747   //
1748
1749   TString name("Nz");
1750   name += fCalibraMode->GetNz(2);
1751   name += "Nrphi";
1752   name += fCalibraMode->GetNrphi(2);
1753   name += "Ngp";
1754   name += fNgroupprf;
1755
1756   if(fNgroupprf != 0){
1757     
1758     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1759                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1760     fPRF2d->SetYTitle("Det/pad groups");
1761     fPRF2d->SetXTitle("Position x/W [pad width units]");
1762     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1763     fPRF2d->SetStats(0);
1764   }
1765   else{
1766     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1767                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1768     fPRF2d->SetYTitle("Det/pad groups");
1769     fPRF2d->SetXTitle("Position x/W [pad width units]");
1770     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1771     fPRF2d->SetStats(0);
1772   }
1773
1774 }
1775
1776 //_____________________________________________________________________________
1777 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1778 {
1779   //
1780   // Create the 2D histos
1781   //
1782
1783   TString name("Nz");
1784   name += fCalibraMode->GetNz(1);
1785   name += "Nrphi";
1786   name += fCalibraMode->GetNrphi(1);
1787   
1788   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1789                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1790                          ,nn,0,nn);
1791   fPH2d->SetYTitle("Det/pad groups");
1792   fPH2d->SetXTitle("time [#mus]");
1793   fPH2d->SetZTitle("<PH> [a.u.]");
1794   fPH2d->SetStats(0);
1795
1796 }
1797 //_____________________________________________________________________________
1798 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1799 {
1800   //
1801   // Create the 2D histos
1802   //
1803
1804   TString name("Nz");
1805   name += fCalibraMode->GetNz(0);
1806   name += "Nrphi";
1807   name += fCalibraMode->GetNrphi(0);
1808
1809   fCH2d = new TH2I("CH2d",(const Char_t *) name
1810                    ,fNumberBinCharge,0,300,nn,0,nn);
1811   fCH2d->SetYTitle("Det/pad groups");
1812   fCH2d->SetXTitle("charge deposit [a.u]");
1813   fCH2d->SetZTitle("counts");
1814   fCH2d->SetStats(0);
1815   fCH2d->Sumw2();
1816
1817 }
1818
1819 //____________Offine tracking in the AliTRDtracker_____________________________
1820 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1821 {
1822   //
1823   // For the offline tracking or mcm tracklets
1824   // This function will be called in the functions UpdateHistogram... 
1825   // to fill the info of a track for the relativ gain calibration
1826   //
1827         
1828   Int_t nb            =  0;   // Nombre de zones traversees
1829   Int_t fd            = -1;   // Premiere zone non nulle
1830   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1831
1832  
1833   
1834   
1835   // See if the track goes through different zones
1836   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1837     if (fAmpTotal[k] > 0.0) {
1838       totalcharge += fAmpTotal[k];
1839       nb++;
1840       if (nb == 1) {
1841         fd = k;
1842       }
1843     }
1844   }
1845
1846   
1847   switch (nb)
1848     { 
1849     case 1:
1850       fNumberUsedCh[0]++;
1851       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1852       if (fHisto2d) {
1853         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1854         //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1855       }
1856       if (fVector2d) {
1857         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1858       }
1859       break;
1860     case 2:
1861       if ((fAmpTotal[fd]   > 0.0) && 
1862           (fAmpTotal[fd+1] > 0.0)) {
1863         // One of the two very big
1864         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1865           if (fHisto2d) {
1866             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1867             //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1868           }
1869           if (fVector2d) {
1870             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1871           }
1872           fNumberUsedCh[1]++;
1873           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1874         }
1875         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1876           if (fHisto2d) {
1877             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1878             //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1879           }
1880           if (fVector2d) {
1881             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1882           }
1883           fNumberUsedCh[1]++;
1884           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1885         }
1886       }
1887       if (fCalibraMode->GetNfragZ(0) > 1) {
1888         if (fAmpTotal[fd] > 0.0) {
1889           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1890             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1891               // One of the two very big
1892               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1893                 if (fHisto2d) {
1894                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1895                   //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1896                 }
1897                 if (fVector2d) {
1898                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1899                 }
1900                 fNumberUsedCh[1]++;
1901                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1902               }
1903               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1904                 if (fHisto2d) {
1905                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1906                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1907                 }
1908                 fNumberUsedCh[1]++;
1909                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1910                 if (fVector2d) {
1911                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1912                 }
1913               }
1914             }
1915           }
1916         }
1917       }
1918       break;
1919     default: break;
1920     }
1921 }
1922 //____________Offine tracking in the AliTRDtracker_____________________________
1923 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1924 {
1925   //
1926   // For the offline tracking or mcm tracklets
1927   // This function will be called in the functions UpdateHistogram... 
1928   // to fill the info of a track for the drift velocity  calibration
1929   //
1930     
1931   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1932   Int_t fd1 = -1; // Premiere zone non nulle
1933   Int_t fd2 = -1; // Deuxieme zone non nulle
1934   Int_t k1  = -1; // Debut de la premiere zone
1935   Int_t k2  = -1; // Debut de la seconde zone
1936
1937   // See if the track goes through different zones
1938   for (Int_t k = 0; k < fTimeMax; k++) {
1939     if (fPHValue[k] > 0.0) {
1940       if (fd1 == -1) {
1941         fd1 = fPHPlace[k];
1942         k1  = k;              
1943       }
1944       if (fPHPlace[k] != fd1) {
1945         if (fd2 == -1) {
1946           k2  = k;
1947           fd2 = fPHPlace[k];
1948           nb  = 2;
1949         }
1950         if (fPHPlace[k] != fd2) {
1951           nb = 3;
1952         }
1953       }
1954     }
1955   }
1956
1957   
1958   switch(nb)
1959     {
1960     case 1:
1961       fNumberUsedPh[0]++;
1962       for (Int_t i = 0; i < fTimeMax; i++) {
1963         if (fHisto2d) {
1964           fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1965         }
1966         if (fVector2d) {
1967           fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1968         }
1969       }
1970       break;
1971     case 2:
1972       if ((fd1 == fd2+1) || 
1973           (fd2 == fd1+1)) {
1974         // One of the two fast all the think
1975         if (k2 > (k1+fDifference)) {
1976           //we choose to fill the fd1 with all the values
1977           fNumberUsedPh[1]++;
1978           for (Int_t i = 0; i < fTimeMax; i++) {
1979             if (fHisto2d) {
1980               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1981             }
1982             if (fVector2d) {
1983               fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1984             }
1985           }
1986         }
1987         if ((k2+fDifference) < fTimeMax) {
1988           //we choose to fill the fd2 with all the values
1989           fNumberUsedPh[1]++;
1990           for (Int_t i = 0; i < fTimeMax; i++) {
1991             if (fHisto2d) {
1992               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1993             }
1994           if (fVector2d) {
1995             fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1996           }
1997           }
1998         }
1999       }
2000       // Two zones voisines sinon rien!
2001       if (fCalibraMode->GetNfragZ(1) > 1) {
2002         // Case 2
2003         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2004           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2005             // One of the two fast all the think
2006             if (k2 > (k1+fDifference)) {
2007               //we choose to fill the fd1 with all the values
2008               fNumberUsedPh[1]++;
2009               for (Int_t i = 0; i < fTimeMax; i++) {
2010                 if (fHisto2d) {
2011                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2012                 }
2013                 if (fVector2d) {
2014                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2015                 }
2016               }
2017             }
2018             if ((k2+fDifference) < fTimeMax) {
2019               //we choose to fill the fd2 with all the values
2020               fNumberUsedPh[1]++;
2021               for (Int_t i = 0; i < fTimeMax; i++) {
2022                 if (fHisto2d) {
2023                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2024                 }
2025                 if (fVector2d) {
2026                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2027                 }
2028               }
2029             }
2030           }
2031         }
2032         // Two zones voisines sinon rien!
2033         // Case 3
2034         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2035           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2036             // One of the two fast all the think
2037             if (k2 > (k1 + fDifference)) {
2038               //we choose to fill the fd1 with all the values
2039               fNumberUsedPh[1]++;
2040               for (Int_t i = 0; i < fTimeMax; i++) {
2041                 if (fHisto2d) {
2042                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2043                 }
2044                 if (fVector2d) {
2045                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2046                 }
2047               }
2048             }
2049             if ((k2+fDifference) < fTimeMax) {
2050               //we choose to fill the fd2 with all the values
2051               fNumberUsedPh[1]++;
2052               for (Int_t i = 0; i < fTimeMax; i++) {
2053                 if (fHisto2d) {
2054                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2055                 }
2056                 if (fVector2d) {
2057                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2058                 }
2059               }
2060             }
2061           }
2062         }
2063       }
2064       break;
2065     default: break;
2066     } 
2067 }
2068 //____________Offine tracking in the AliTRDtracker_____________________________
2069 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
2070 {
2071   //
2072   // For the offline tracking
2073   // This function will be called in the functions UpdateHistogram... 
2074   // to fill the find the parameter P1 of a track for the drift velocity  calibration
2075   //
2076
2077    
2078   //Number of points: if less than 3 return kFALSE
2079   Int_t Npoints = fListClusters->GetEntriesFast();
2080   if(Npoints <= 2) return kFALSE;
2081
2082   //Variables
2083   TLinearFitter linearFitterTracklet      = TLinearFitter(2,"pol1");        // TLinearFitter per tracklet
2084   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
2085   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
2086   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
2087   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
2088   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
2089   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!  
2090   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
2091   Double_t pointError                 = 0.0;                                // error after straight line fit 
2092   Int_t    detector                   = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
2093   Int_t    snpright                   = 1;                                  // if we took in the middle snp
2094   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
2095   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
2096   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
2097   Int_t     nbli                      = 0;                                  // number linear fitter points
2098   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2099
2100   linearFitterTracklet.StoreData(kFALSE);
2101   linearFitterTracklet.ClearPoints();
2102   
2103   //if more than one row
2104   Int_t    rowp                       = -1;                              // if it crosses a pad row
2105
2106   //tiltingangle
2107   tiltingangle                        = padplane->GetTiltingAngle();
2108   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2109
2110   //Fill with points
2111   for(Int_t k = 0; k < Npoints; k++){
2112     
2113     AliTRDcluster *cl                 = (AliTRDcluster *) fListClusters->At(k);
2114     Double_t ycluster                 = cl->GetY();
2115     Int_t time                        = cl->GetLocalTimeBin();
2116     Double_t timeis                   = time/fSf;
2117     //See if cross two pad rows
2118     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
2119     if(k==0) rowp                     = row;
2120     if(row != rowp) crossrow          = 1;
2121     //Take in the middle of the chamber
2122     //FollowBack
2123     if(time > (Int_t) 10) {
2124     //Follow
2125     //if(time < (Int_t) 11) {
2126       z   = cl->GetZ();
2127       y   = cl->GetY();  
2128       snp = fPar2[time];
2129       tgl = fPar3[time];
2130     }
2131     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2132     nbli++;
2133   }
2134   //FollowBack
2135   if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
2136   //Follow
2137   //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
2138   if(nbli <= 2) return kFALSE; 
2139   
2140   // Do the straight line fit now
2141   TVectorD pars;
2142   linearFitterTracklet.Eval();
2143   linearFitterTracklet.GetParameters(pars);
2144   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2145   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
2146   dydt  = pars[1]; 
2147  
2148   if( TMath::Abs(snp) <  1.){
2149     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2150   } 
2151   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2152
2153   if(fDebugLevel > 0){
2154     if ( !fDebugStreamer ) {
2155       //debug stream
2156       TDirectory *backup = gDirectory;
2157       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2158       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2159     } 
2160     
2161     (* fDebugStreamer) << "VDRIFT0"<<
2162       "Npoints="<<Npoints<<
2163       "\n"; 
2164   
2165     
2166     (* fDebugStreamer) << "VDRIFT"<<
2167       "snpright="<<snpright<<
2168       "Npoints="<<Npoints<<
2169       "nbli="<<nbli<<
2170       "detector="<<detector<<
2171       "snp="<<snp<<
2172       "tnp="<<tnp<<
2173       "tgl="<<tgl<<
2174       "tnt="<<tnt<<
2175       "y="<<y<<
2176       "z="<<z<<
2177       "dydt="<<dydt<<
2178       "dzdx="<<dzdx<<
2179       "crossrow="<<crossrow<<
2180       "errorpar="<<errorpar<<
2181       "pointError="<<pointError<<
2182       "\n";     
2183
2184   }
2185   
2186   if(Npoints < fNumberClusters) return kFALSE;
2187   if(snpright == 0) return kFALSE;
2188   if(pointError >= 0.1) return kFALSE;
2189   if(crossrow == 1) return kFALSE;
2190   
2191   if(fLinearFitterOn){
2192     //Add to the linear fitter of the detector
2193     if( TMath::Abs(snp) <  1.){
2194       Double_t x = tnp-dzdx*tnt; 
2195       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2196       if(fLinearFitterDebugOn) {
2197         fLinearVdriftFit->Update(detector,x,pars[1]);
2198       }
2199       fEntriesLinearFitter[detector]++;
2200     }
2201   }
2202   //AliInfo("End of FindP1TrackPH with success!")
2203   return kTRUE;
2204
2205 }
2206 //____________Offine tracking in the AliTRDtracker_____________________________
2207 Bool_t AliTRDCalibraFillHisto::HandlePRF()
2208 {
2209   //
2210   // For the offline tracking
2211   // Fit the tracklet with a line and take the position as reference for the PRF
2212   //
2213
2214   //Number of points
2215   Int_t Npoints  = fListClusters->GetEntriesFast();                         // number of total points
2216   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit 
2217   Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
2218  
2219
2220   // To see the difference due to the fit
2221   Double_t *padPositions;
2222   padPositions = new Double_t[Npoints];
2223   for(Int_t k = 0; k < Npoints; k++){
2224     padPositions[k] = 0.0;
2225   } 
2226
2227
2228   //Find the position by a fit
2229   TLinearFitter fitter(2,"pol1");
2230   fitter.StoreData(kFALSE);
2231   fitter.ClearPoints();
2232   for(Int_t k = 0;  k < Npoints; k++){
2233     //Take the cluster
2234     AliTRDcluster *cl  = (AliTRDcluster *) fListClusters->At(k);
2235     Short_t  *signals  = cl->GetSignals();
2236     Double_t     time  = cl->GetLocalTimeBin();
2237     //Calculate x if possible 
2238     Float_t xcenter    = 0.0;    
2239     Bool_t  echec      = kTRUE;   
2240     if((time<=7) || (time>=21)) continue; 
2241     // Center 3 balanced: position with the center of the pad
2242     if ((((Float_t) signals[3]) > 0.0) && 
2243         (((Float_t) signals[2]) > 0.0) && 
2244         (((Float_t) signals[4]) > 0.0)) {
2245       echec = kFALSE;
2246       // Security if the denomiateur is 0 
2247       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / 
2248            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2249         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2250           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) 
2251                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2252       }
2253       else {
2254         echec = kTRUE;
2255       }
2256     }
2257     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2258     if(echec) continue;
2259     //if no echec: calculate with the position of the pad
2260     // Position of the cluster
2261     Double_t       padPosition = xcenter +  cl->GetPad();
2262     padPositions[k]            = padPosition;
2263     Nb3pc++;
2264     fitter.AddPoint(&time, padPosition,1);
2265   }//clusters loop
2266
2267   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
2268   if(Nb3pc < 3) return kFALSE;
2269   fitter.Eval();
2270   TVectorD line(2);
2271   fitter.GetParameters(line);
2272   Float_t  pointError  = -1.0;
2273   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
2274   
2275
2276   // Now fill the PRF  
2277   for(Int_t k = 0;  k < Npoints; k++){
2278     //Take the cluster
2279     AliTRDcluster *cl      = (AliTRDcluster *) fListClusters->At(k);
2280     Short_t  *signals      = cl->GetSignals();              // signal
2281     Double_t     time      = cl->GetLocalTimeBin();         // time bin
2282     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
2283     Float_t padPos         = cl->GetPad();                  // middle pad
2284     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit 
2285     Float_t ycenter        = 0.0;                           // relative center charge
2286     Float_t ymin           = 0.0;                           // relative left charge
2287     Float_t ymax           = 0.0;                           // relative right charge
2288     Double_t tgl           = fPar3[(Int_t)time];            // dz/dl and not dz/dx
2289     Double_t pt            = fPar4[(Int_t)time];            // pt
2290     Float_t  dzdx          = 0.0;                           // dzdx
2291
2292
2293     //Requiere simply two pads clusters at least
2294     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2295        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2296       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2297       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2298       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
2299       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum; 
2300     }
2301     
2302     //calibration group
2303     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2304     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2305     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
2306     Double_t  snp          = fPar2[(Int_t)time];                        // sin angle in xy plan
2307     Float_t   xcl          = cl->GetY();                                // y cluster
2308     Float_t   qcl          = cl->GetQ();                                // charge cluster 
2309     Int_t     plane        = GetPlane(detector);                        // plane 
2310     Int_t     chamber      = GetChamber(detector);                      // chamber  
2311     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2312     Double_t  x            = dpad;                                      // reconstructed position moved
2313     Float_t   Ep           = pointError;                                // error of fit
2314     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2315     Float_t   signal3      = (Float_t)signals[3];                       // signal
2316     Float_t   signal2      = (Float_t)signals[2];                       // signal
2317     Float_t   signal4      = (Float_t)signals[4];                       // signal
2318     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2319     Float_t tnp            = 0.0;
2320     if(TMath::Abs(snp) < 1.0){
2321       tnp = snp / (TMath::Sqrt(1-snp*snp));
2322       dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2323     }
2324    
2325
2326     if(fDebugLevel > 0){
2327       if ( !fDebugStreamer ) {
2328         //debug stream
2329         TDirectory *backup = gDirectory;
2330         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2331         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2332       }     
2333       
2334       (* fDebugStreamer) << "PRF0"<<
2335         "caligroup="<<caligroup<<
2336         "detector="<<detector<<
2337         "plane="<<plane<<
2338         "chamber="<<chamber<<
2339         "Npoints="<<Npoints<<
2340         "Np="<<Nb3pc<<
2341         "Ep="<<Ep<<
2342         "snp="<<snp<<
2343         "tnp="<<tnp<<    
2344         "tgl="<<tgl<<  
2345         "dzdx="<<dzdx<<  
2346         "pt="<<pt<<    
2347         "padPos="<<padPos<<
2348         "padPosition="<<padPositions[k]<<
2349         "padPosTracklet="<<padPosTracklet<<
2350         "x="<<x<<
2351         "ycenter="<<ycenter<<
2352         "ymin="<<ymin<<
2353         "ymax="<<ymax<<
2354         "xcl="<<xcl<<
2355         "qcl="<<qcl<< 
2356         "signal1="<<signal1<<    
2357         "signal2="<<signal2<<
2358         "signal3="<<signal3<<
2359         "signal4="<<signal4<<
2360         "signal5="<<signal5<<
2361         "time="<<time<<
2362         "\n";     
2363       x = xdiff;
2364       Int_t type=0;
2365       Float_t y = ycenter;
2366       (* fDebugStreamer) << "PRFALL"<<
2367         "caligroup="<<caligroup<<
2368         "detector="<<detector<<
2369         "plane="<<plane<<
2370         "chamber="<<chamber<<
2371         "Npoints="<<Npoints<<
2372         "Np="<<Nb3pc<<
2373         "Ep="<<Ep<<
2374         "type="<<type<<
2375         "snp="<<snp<<
2376         "tnp="<<tnp<<
2377         "tgl="<<tgl<<  
2378         "dzdx="<<dzdx<< 
2379         "pt="<<pt<< 
2380         "padPos="<<padPos<<
2381         "padPosition="<<padPositions[k]<<
2382         "padPosTracklet="<<padPosTracklet<<
2383         "x="<<x<<
2384         "y="<<y<<           
2385         "xcl="<<xcl<<
2386         "qcl="<<qcl<<
2387         "signal1="<<signal1<<
2388         "signal2="<<signal2<<
2389         "signal3="<<signal3<<
2390         "signal4="<<signal4<<
2391         "signal5="<<signal5<<
2392         "time="<<time<<
2393         "\n";
2394       x=-(xdiff+1);
2395       y = ymin;
2396       type=-1;
2397       (* fDebugStreamer) << "PRFALL"<<
2398         "caligroup="<<caligroup<<
2399         "detector="<<detector<<
2400         "plane="<<plane<<
2401         "chamber="<<chamber<<
2402         "Npoints="<<Npoints<<
2403         "Np="<<Nb3pc<<
2404         "Ep="<<Ep<<
2405         "type="<<type<<
2406         "snp="<<snp<<
2407         "tnp="<<tnp<<
2408         "tgl="<<tgl<<  
2409         "dzdx="<<dzdx<< 
2410         "pt="<<pt<< 
2411         "padPos="<<padPos<<
2412         "padPosition="<<padPositions[k]<<
2413         "padPosTracklet="<<padPosTracklet<<
2414         "x="<<x<<
2415         "y="<<y<<
2416         "xcl="<<xcl<<
2417         "qcl="<<qcl<<
2418         "signal1="<<signal1<<
2419         "signal2="<<signal2<<
2420         "signal3="<<signal3<<
2421         "signal4="<<signal4<<
2422         "signal5="<<signal5<<
2423         "time="<<time<<
2424         "\n";
2425       x=1-xdiff;
2426       y = ymax;
2427       type=1;
2428       
2429       (* fDebugStreamer) << "PRFALL"<<
2430         "caligroup="<<caligroup<<
2431         "detector="<<detector<<
2432         "plane="<<plane<<
2433         "chamber="<<chamber<<
2434         "Npoints="<<Npoints<<
2435         "Np="<<Nb3pc<<
2436         "Ep="<<Ep<<
2437         "type="<<type<<
2438         "snp="<<snp<<
2439         "tnp="<<tnp<<   
2440         "tgl="<<tgl<<  
2441         "dzdx="<<dzdx<< 
2442         "pt="<<pt<< 
2443         "padPos="<<padPos<<
2444         "padPosition="<<padPositions[k]<<
2445         "padPosTracklet="<<padPosTracklet<<
2446         "x="<<x<<
2447         "y="<<y<<
2448         "xcl="<<xcl<<
2449         "qcl="<<qcl<<
2450         "signal1="<<signal1<<
2451         "signal2="<<signal2<<
2452         "signal3="<<signal3<<
2453         "signal4="<<signal4<<
2454         "signal5="<<signal5<<
2455         "time="<<time<<
2456         "\n";
2457       
2458     }
2459     
2460     // some cuts
2461     if(Npoints < fNumberClusters) continue;
2462     if(Nb3pc <= 5) continue;
2463     if((time >= 21) || (time < 7)) continue;
2464     if(TMath::Abs(snp) >= 1.0) continue;
2465     if(qcl < 80) continue; 
2466     
2467     Bool_t echec   = kFALSE;
2468     Double_t shift = 0.0;
2469     //Calculate the shift in x coresponding to this tnp
2470     if(fNgroupprf != 0.0){
2471       shift      = -3.0*(fNgroupprf-1)-1.5;
2472       Double_t limithigh  = -0.2*(fNgroupprf-1);
2473       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2474       else{
2475         while(tnp > limithigh){
2476           limithigh += 0.2;
2477           shift += 3.0;
2478         }
2479       }
2480     }
2481     if (fHisto2d && !echec) {
2482       if(TMath::Abs(dpad) < 1.5) {
2483         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2484         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2485       }
2486       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2487         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2488         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2489       }
2490       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2491         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2492         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2493       }
2494     }
2495     //Not equivalent anymore here!
2496     if (fVector2d && !echec) {
2497       if(TMath::Abs(dpad) < 1.5) {
2498         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2499         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2500       }
2501       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2502         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2503         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2504       }
2505       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
2506         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2507         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2508       }
2509     }
2510   }
2511   return kTRUE;
2512   
2513 }
2514 //____________Offine tracking in the AliTRDtracker_____________________________
2515 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2516 {
2517   //
2518   // For the offline tracking
2519   // This function will be called in the functions UpdateHistogram... 
2520   // to fill the find the parameter P1 of a track for the drift velocity  calibration
2521   //
2522
2523     
2524   //Number of points: if less than 3 return kFALSE
2525   Int_t Npoints = index1-index0;
2526   if(Npoints <= 2) return kFALSE;
2527
2528   //Variables
2529   TLinearFitter linearFitterTracklet  = TLinearFitter(2,"pol1");            // TLinearFitter per tracklet
2530   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
2531   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
2532   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
2533   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
2534   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
2535   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!  
2536   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
2537   Double_t pointError                 = 0.0;                                // error after straight line fit 
2538   Int_t    detector                   = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
2539   //Int_t    snpright                   = 1;                                  // if we took in the middle snp
2540   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
2541   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
2542   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
2543   Int_t     nbli                      = 0;                                  // number linear fitter points
2544   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2545
2546   linearFitterTracklet.StoreData(kFALSE);
2547   linearFitterTracklet.ClearPoints();
2548   
2549   //if more than one row
2550   Int_t    rowp                       = -1;                              // if it crosses a pad row
2551
2552   //tiltingangle
2553   tiltingangle                        = padplane->GetTiltingAngle();
2554   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2555
2556   //Fill with points
2557   for(Int_t k = 0; k < Npoints; k++){
2558     
2559     AliTRDcluster *cl                 = (AliTRDcluster *) t->GetCluster(k+index0);
2560     Double_t ycluster                 = cl->GetY();
2561     Int_t time                        = cl->GetLocalTimeBin();
2562     Double_t timeis                   = time/fSf;
2563     //See if cross two pad rows
2564     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
2565     if(k==0) rowp                     = row;
2566     if(row != rowp) crossrow          = 1;
2567     //Take in the middle of the chamber
2568     //FollowBack
2569     //if(time > (Int_t) 10) {
2570     //Follow
2571     if(time < (Int_t) 11) {
2572       z   = cl->GetZ();
2573       y   = cl->GetY();  
2574     }
2575     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2576     nbli++;
2577   }
2578
2579   // take now the snp, tnp and tgl from the track
2580   snp = t->GetSnpPlane(GetPlane(detector));
2581   tgl = t->GetTglPlane(GetPlane(detector));
2582
2583   //FollowBack
2584   //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0;
2585   //Follow
2586   //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0;
2587   if(nbli <= 2) return kFALSE; 
2588   
2589   // Do the straight line fit now
2590   TVectorD pars;
2591   linearFitterTracklet.Eval();
2592   linearFitterTracklet.GetParameters(pars);
2593   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2594   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
2595   dydt  = pars[1]; 
2596  
2597   if( TMath::Abs(snp) <  1.){
2598     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2599   } 
2600   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2601
2602   if(fDebugLevel > 0){
2603     if ( !fDebugStreamer ) {
2604       //debug stream
2605       TDirectory *backup = gDirectory;
2606       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2607       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2608     } 
2609     
2610         
2611     (* fDebugStreamer) << "VDRIFT"<<
2612       //"snpright="<<snpright<<
2613       "Npoints="<<Npoints<<
2614       "nbli="<<nbli<<
2615       "detector="<<detector<<
2616       "snp="<<snp<<
2617       "tnp="<<tnp<<
2618       "tgl="<<tgl<<
2619       "tnt="<<tnt<<
2620       "y="<<y<<
2621       "z="<<z<<
2622       "dydt="<<dydt<<
2623       "dzdx="<<dzdx<<
2624       "crossrow="<<crossrow<<
2625       "errorpar="<<errorpar<<
2626       "pointError="<<pointError<<
2627       "\n";     
2628
2629   }
2630   
2631   if(Npoints < fNumberClusters) return kFALSE;
2632   //if(snpright == 0) return kFALSE;
2633   if(pointError >= 0.1) return kFALSE;
2634   if(crossrow == 1) return kFALSE;
2635   
2636   if(fLinearFitterOn){
2637     //Add to the linear fitter of the detector
2638     if( TMath::Abs(snp) <  1.){
2639       Double_t x = tnp-dzdx*tnt; 
2640       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2641       if(fLinearFitterDebugOn) {
2642         fLinearVdriftFit->Update(detector,x,pars[1]);
2643       }
2644       fEntriesLinearFitter[detector]++;
2645     }
2646   }
2647   //AliInfo("End of FindP1TrackPH with success!")
2648   return kTRUE;
2649
2650 }
2651 //____________Offine tracking in the AliTRDtracker_____________________________
2652 Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2653 {
2654   //
2655   // For the offline tracking
2656   // Fit the tracklet with a line and take the position as reference for the PRF
2657   //
2658
2659   //Number of points
2660   Int_t Npoints  = index1-index0;                                           // number of total points
2661   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit 
2662   Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
2663  
2664
2665   // To see the difference due to the fit
2666   Double_t *padPositions;
2667   padPositions = new Double_t[Npoints];
2668   for(Int_t k = 0; k < Npoints; k++){
2669     padPositions[k] = 0.0;
2670   } 
2671
2672
2673   //Find the position by a fit
2674   TLinearFitter fitter(2,"pol1");
2675   fitter.StoreData(kFALSE);
2676   fitter.ClearPoints();
2677   for(Int_t k = 0;  k < Npoints; k++){
2678     //Take the cluster
2679     AliTRDcluster *cl  = (AliTRDcluster *) t->GetCluster(k+index0);
2680     Short_t  *signals  = cl->GetSignals();
2681     Double_t     time  = cl->GetLocalTimeBin();
2682     //Calculate x if possible 
2683     Float_t xcenter    = 0.0;    
2684     Bool_t  echec      = kTRUE;   
2685     if((time<=7) || (time>=21)) continue; 
2686     // Center 3 balanced: position with the center of the pad
2687     if ((((Float_t) signals[3]) > 0.0) && 
2688         (((Float_t) signals[2]) > 0.0) && 
2689         (((Float_t) signals[4]) > 0.0)) {
2690       echec = kFALSE;
2691       // Security if the denomiateur is 0 
2692       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) / 
2693            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2694         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2695           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) 
2696                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2697       }
2698       else {
2699         echec = kTRUE;
2700       }
2701     }
2702     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2703     if(echec) continue;
2704     //if no echec: calculate with the position of the pad
2705     // Position of the cluster
2706     Double_t       padPosition = xcenter +  cl->GetPad();
2707     padPositions[k]            = padPosition;
2708     Nb3pc++;
2709     fitter.AddPoint(&time, padPosition,1);
2710   }//clusters loop
2711
2712   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
2713   if(Nb3pc < 3) return kFALSE;
2714   fitter.Eval();
2715   TVectorD line(2);
2716   fitter.GetParameters(line);
2717   Float_t  pointError  = -1.0;
2718   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
2719   
2720   // Take the tgl and snp with the track t now
2721   Double_t  tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
2722   Double_t  snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
2723   Float_t  dzdx = 0.0;                                // dzdx
2724   Float_t  tnp  = 0.0;
2725   if(TMath::Abs(snp) < 1.0){
2726     tnp = snp / (TMath::Sqrt(1-snp*snp));
2727     dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2728   }
2729
2730
2731   // Now fill the PRF  
2732   for(Int_t k = 0;  k < Npoints; k++){
2733     //Take the cluster
2734     AliTRDcluster *cl      = (AliTRDcluster *) t->GetCluster(k+index0);
2735     Short_t  *signals      = cl->GetSignals();              // signal
2736     Double_t     time      = cl->GetLocalTimeBin();         // time bin
2737     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
2738     Float_t padPos         = cl->GetPad();                  // middle pad
2739     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit 
2740     Float_t ycenter        = 0.0;                           // relative center charge
2741     Float_t ymin           = 0.0;                           // relative left charge
2742     Float_t ymax           = 0.0;                           // relative right charge
2743   
2744
2745
2746     //Requiere simply two pads clusters at least
2747     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2748        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2749       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2750       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2751       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
2752       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum; 
2753     }
2754     
2755     //calibration group
2756     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2757     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2758     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
2759     Float_t   xcl          = cl->GetY();                                // y cluster
2760     Float_t   qcl          = cl->GetQ();                                // charge cluster 
2761     Int_t     plane        = GetPlane(detector);                        // plane 
2762     Int_t     chamber      = GetChamber(detector);                      // chamber  
2763     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2764     Double_t  x            = dpad;                                      // reconstructed position moved
2765     Float_t   Ep           = pointError;                                // error of fit
2766     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2767     Float_t   signal3      = (Float_t)signals[3];                       // signal
2768     Float_t   signal2      = (Float_t)signals[2];                       // signal
2769     Float_t   signal4      = (Float_t)signals[4];                       // signal
2770     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2771    
2772    
2773
2774     if(fDebugLevel > 0){
2775       if ( !fDebugStreamer ) {
2776         //debug stream
2777         TDirectory *backup = gDirectory;
2778         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2779         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2780       }     
2781       
2782       (* fDebugStreamer) << "PRF0"<<
2783         "caligroup="<<caligroup<<
2784         "detector="<<detector<<
2785         "plane="<<plane<<
2786         "chamber="<<chamber<<
2787         "Npoints="<<Npoints<<
2788         "Np="<<Nb3pc<<
2789         "Ep="<<Ep<<
2790         "snp="<<snp<<
2791         "tnp="<<tnp<<    
2792         "tgl="<<tgl<<  
2793         "dzdx="<<dzdx<<  
2794         "padPos="<<padPos<<
2795         "padPosition="<<padPositions[k]<<
2796         "padPosTracklet="<<padPosTracklet<<
2797         "x="<<x<<
2798         "ycenter="<<ycenter<<
2799         "ymin="<<ymin<<
2800         "ymax="<<ymax<<
2801         "xcl="<<xcl<<
2802         "qcl="<<qcl<< 
2803         "signal1="<<signal1<<    
2804         "signal2="<<signal2<<
2805         "signal3="<<signal3<<
2806         "signal4="<<signal4<<
2807         "signal5="<<signal5<<
2808         "time="<<time<<
2809         "\n";     
2810       x = xdiff;
2811       Int_t type=0;
2812       Float_t y = ycenter;
2813       (* fDebugStreamer) << "PRFALL"<<
2814         "caligroup="<<caligroup<<
2815         "detector="<<detector<<
2816         "plane="<<plane<<
2817         "chamber="<<chamber<<
2818         "Npoints="<<Npoints<<
2819         "Np="<<Nb3pc<<
2820         "Ep="<<Ep<<
2821         "type="<<type<<
2822         "snp="<<snp<<
2823         "tnp="<<tnp<<
2824         "tgl="<<tgl<<  
2825         "dzdx="<<dzdx<< 
2826         "padPos="<<padPos<<
2827         "padPosition="<<padPositions[k]<<
2828         "padPosTracklet="<<padPosTracklet<<
2829         "x="<<x<<
2830         "y="<<y<<           
2831         "xcl="<<xcl<<
2832         "qcl="<<qcl<<
2833         "signal1="<<signal1<<
2834         "signal2="<<signal2<<
2835         "signal3="<<signal3<<
2836         "signal4="<<signal4<<
2837         "signal5="<<signal5<<
2838         "time="<<time<<
2839         "\n";
2840       x=-(xdiff+1);
2841       y = ymin;
2842       type=-1;
2843       (* fDebugStreamer) << "PRFALL"<<
2844         "caligroup="<<caligroup<<
2845         "detector="<<detector<<
2846         "plane="<<plane<<
2847         "chamber="<<chamber<<
2848         "Npoints="<<Npoints<<
2849         "Np="<<Nb3pc<<
2850         "Ep="<<Ep<<
2851         "type="<<type<<
2852         "snp="<<snp<<
2853         "tnp="<<tnp<<
2854         "tgl="<<tgl<<  
2855         "dzdx="<<dzdx<< 
2856         "padPos="<<padPos<<
2857         "padPosition="<<padPositions[k]<<
2858         "padPosTracklet="<<padPosTracklet<<
2859         "x="<<x<<
2860         "y="<<y<<
2861         "xcl="<<xcl<<
2862         "qcl="<<qcl<<
2863         "signal1="<<signal1<<
2864         "signal2="<<signal2<<
2865         "signal3="<<signal3<<
2866         "signal4="<<signal4<<
2867         "signal5="<<signal5<<
2868         "time="<<time<<
2869         "\n";
2870       x=1-xdiff;
2871       y = ymax;
2872       type=1;
2873       
2874       (* fDebugStreamer) << "PRFALL"<<
2875         "caligroup="<<caligroup<<
2876         "detector="<<detector<<
2877         "plane="<<plane<<
2878         "chamber="<<chamber<<
2879         "Npoints="<<Npoints<<
2880         "Np="<<Nb3pc<<
2881         "Ep="<<Ep<<
2882         "type="<<type<<
2883         "snp="<<snp<<
2884         "tnp="<<tnp<<   
2885         "tgl="<<tgl<<  
2886         "dzdx="<<dzdx<< 
2887         "padPos="<<padPos<<
2888         "padPosition="<<padPositions[k]<<
2889         "padPosTracklet="<<padPosTracklet<<
2890         "x="<<x<<
2891         "y="<<y<<
2892         "xcl="<<xcl<<
2893         "qcl="<<qcl<<
2894         "signal1="<<signal1<<
2895         "signal2="<<signal2<<
2896         "signal3="<<signal3<<
2897         "signal4="<<signal4<<
2898         "signal5="<<signal5<<
2899         "time="<<time<<
2900         "\n";
2901       
2902     }
2903     
2904     // some cuts
2905     if(Npoints < fNumberClusters) continue;
2906     if(Nb3pc <= 5) continue;
2907     if((time >= 21) || (time < 7)) continue;
2908     if(TMath::Abs(snp) >= 1.0) continue;
2909     if(qcl < 80) continue; 
2910     
2911     Bool_t echec   = kFALSE;
2912     Double_t shift = 0.0;
2913     //Calculate the shift in x coresponding to this tnp
2914     if(fNgroupprf != 0.0){
2915       shift      = -3.0*(fNgroupprf-1)-1.5;
2916       Double_t limithigh  = -0.2*(fNgroupprf-1);
2917       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2918       else{
2919         while(tnp > limithigh){
2920           limithigh += 0.2;
2921           shift += 3.0;
2922         }
2923       }
2924     }
2925     if (fHisto2d && !echec) {
2926       if(TMath::Abs(dpad) < 1.5) {
2927         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2928         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2929       }
2930       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2931         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2932         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2933       }
2934       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2935         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2936         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2937       }
2938     }
2939     //Not equivalent anymore here!
2940     if (fVector2d && !echec) {
2941       if(TMath::Abs(dpad) < 1.5) {
2942         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2943         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2944       }
2945       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2946         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2947         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2948       }
2949       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
2950         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2951         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2952       }
2953     }
2954   }
2955   return kTRUE;
2956   
2957 }
2958 //
2959 //____________Some basic geometry function_____________________________________
2960 //
2961 //_____________________________________________________________________________
2962 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2963 {
2964   //
2965   // Reconstruct the plane number from the detector number
2966   //
2967
2968   return ((Int_t) (d % 6));
2969
2970 }
2971
2972 //_____________________________________________________________________________
2973 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2974 {
2975   //
2976   // Reconstruct the chamber number from the detector number
2977   //
2978   Int_t fgkNplan = 6;
2979
2980   return ((Int_t) (d % 30) / fgkNplan);
2981
2982 }
2983
2984 //_____________________________________________________________________________
2985 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2986 {
2987   //
2988   // Reconstruct the sector number from the detector number
2989   //
2990   Int_t fg = 30;
2991
2992   return ((Int_t) (d / fg));
2993
2994 }
2995 //_____________________________________________________________________
2996 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2997 {
2998     //
2999     // return pointer to fPH2d TProfile2D
3000     // create a new TProfile2D if it doesn't exist allready
3001     //
3002     if ( fPH2d )
3003         return fPH2d;
3004
3005     // Some parameters
3006     fTimeMax = nbtimebin;
3007     fSf      = samplefrequency;
3008   
3009     /*
3010       AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
3011       ,fCalibraMode->GetNz(1)
3012       ,fCalibraMode->GetNrphi(1)));
3013       
3014       // Calcul the number of Xbins
3015       Int_t Ntotal1 = 0;
3016       fCalibraMode->ModePadCalibration(2,1);
3017       fCalibraMode->ModePadFragmentation(0,2,0,1);
3018       fCalibraMode->SetDetChamb2(1);
3019       Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
3020       fCalibraMode->ModePadCalibration(0,1);
3021       fCalibraMode->ModePadFragmentation(0,0,0,1);
3022       fCalibraMode->SetDetChamb0(1);
3023       Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
3024       AliInfo(Form("Total number of Xbins: %d",Ntotal1));
3025       
3026       CreatePH2d(Ntotal1);
3027     */  
3028
3029     CreatePH2d(540);
3030
3031     return fPH2d;
3032 }
3033 //_____________________________________________________________________
3034 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3035 {
3036     //
3037     // return pointer to TLinearFitter Calibration
3038     // if force is true create a new TLinearFitter if it doesn't exist allready
3039     //
3040
3041   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3042     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3043   }
3044
3045   // if we are forced and TLinearFitter doesn't yet exist create it
3046
3047   // new TLinearFitter
3048   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3049   fLinearFitterArray.AddAt(linearfitter,detector);
3050   return linearfitter;
3051 }
3052
3053 //_____________________________________________________________________
3054 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3055 {
3056   //
3057   // FillCH2d: Marian style
3058   // 
3059   
3060   //skip simply the value out of range
3061   if((y>=300.0) || (y<0.0)) return;
3062   
3063   //Calcul the y place
3064   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3065   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3066   
3067   //Fill
3068   fCH2d->GetArray()[place]++;
3069
3070 }
3071
3072 //____________________________________________________________________________
3073 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3074 {
3075   //
3076   // Analyse array of linear fitter because can not be written
3077   // Store two arrays: one with the param the other one with the error param + number of entries
3078   //
3079
3080   for(Int_t k = 0; k < 540; k++){
3081     TLinearFitter *linearfitter = GetLinearFitter(k);
3082     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3083       TVectorD  *par  = new TVectorD(2);
3084       TVectorD   pare = TVectorD(2);
3085       TVectorD  *parE = new TVectorD(3);
3086       linearfitter->Eval();
3087       linearfitter->GetParameters(*par);
3088       linearfitter->GetErrors(pare);
3089       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3090       (*parE)[0] = pare[0]*ppointError;
3091       (*parE)[1] = pare[1]*ppointError;
3092       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3093       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3094       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3095     }
3096   }
3097 }
3098 //________________________________________________________________________________
3099 Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const
3100 {
3101    // Partie entiere of the (x+0.5)
3102    
3103   int i;
3104   if (x >= (-0.5)) {
3105     i = int(x + 0.5);
3106     //if (x + 0.5 == Float_t(i) && i & 1) i--;
3107   } else {
3108     i = int(x - 0.5);
3109     //if (x - 0.5 == Float_t(i) && i & 1) i++;
3110     if((x-0.5)==i) i++;
3111
3112   }
3113   return i;
3114 }