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