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