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