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