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