changes from fzhou
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibTask.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /////////////////////////////////////////////////////////////////////////////////
18 //                                                                             
19 // AliTRDCalibTask                                                               
20 //                                                                             
21 // Offline TRD calibration task                                
22 //                        
23 // Author:
24 //   R. Bailhache (R.Bailhache@gsi.de)
25 //                            
26 //////////////////////////////////////////////////////////////////////////////////////
27
28 #include <iostream>
29 using namespace std;
30 #include "Riostream.h"
31 #include "TChain.h"
32 #include "TTree.h"
33 #include "TFile.h"
34 #include "TProfile2D.h"
35 #include "TH2I.h"
36 #include "TH1F.h"
37 #include "TH2F.h"
38 #include "TH2S.h"
39 #include "TList.h"
40 #include "TMath.h"
41 #include "TObject.h"
42 #include "TObjArray.h"
43 #include "TString.h"
44 #include "TCanvas.h"
45 #include "TLegend.h"
46 #include "TStyle.h"
47 #include "TLine.h"
48 #include "TIterator.h"
49 #include "TLinearFitter.h"
50 #include "TVectorD.h"
51
52 #include "AliAnalysisTask.h"
53 #include "AliAnalysisManager.h"
54
55 #include "AliExternalTrackParam.h"
56 #include "AliESDVertex.h"
57 #include "AliESDEvent.h"
58 #include "AliESDfriend.h"
59 #include "AliCentrality.h"
60 #include "AliESDInputHandler.h"
61 #include "AliESDtrack.h"
62 #include "AliESDfriendTrack.h"
63 #include "AliTRDtrackV1.h"
64 #include "AliTRDseedV1.h"
65 #include "AliTRDcluster.h"
66 #include "AliTRDgeometry.h"
67 #include "AliESDtrackCuts.h"
68 #include "AliESDVertex.h"
69 #include "AliTRDCalDet.h"
70
71 #include "AliTRDCalibraVector.h"
72 #include "AliTRDCalibraFillHisto.h"
73 #include "AliTRDCalibraVdriftLinearFit.h" 
74 #include "AliTRDCalibraExbAltFit.h" 
75
76 #include "AliTRDcalibDB.h"
77 #include "AliCDBId.h"
78 #include "AliLog.h"
79
80 #include "AliTRDCalibChamberStatus.h"
81 #include "AliTRDdEdxBaseUtils.h"
82 #include "AliTRDdEdxCalibHistArray.h"
83 #include "AliTRDdEdxCalibUtils.h"
84
85 #include "AliTRDCalibTask.h"
86
87
88 ClassImp(AliTRDCalibTask)
89
90 //________________________________________________________________________
91   AliTRDCalibTask::AliTRDCalibTask(const char *name) 
92     : AliAnalysisTaskSE(name), fESD(0),
93       fkEsdTrack(0),
94       fFriendTrack(0),
95       fCalibObject(0),
96       fTrdTrack(0),
97       fCl(0),
98       fListHist(0),
99       fTRDCalibraFillHisto(0),
100       fTRDChamberStatus(0),
101       fNEvents(0),
102       fNEventsInput(0),
103       fNbTRDTrack(0),
104       fNbTRDTrackOffline(0),
105       fNbTRDTrackStandalone(0),
106       fNbTPCTRDtrack(0),
107       fNbGoodTracks(0),
108       fNbTimeBin(0),
109       fNbTimeBinOffline(0),
110       fNbTimeBinStandalone(0),
111       fNbClusters(0),
112       fNbClustersOffline(0),
113       fNbClustersStandalone(0),
114       fNbTracklets(0),
115       fNbTrackletsOffline(0),
116       fNbTrackletsStandalone(0),
117       fAbsoluteGain(0),
118       fCH2dSum(0),
119       fPH2dSum(0),
120       fCH2dSM(0),
121       fPH2dSM(0),
122       fCH2dTest(0),
123       fPH2dTest(0),
124       fLinearVdriftTest(0),
125       fOnInstance(kTRUE),
126       fHisto2d(kTRUE),
127       fVector2d(kFALSE),
128       fVdriftLinear(kTRUE),
129       fExbAlt(kFALSE),
130       fDebugLevelTRDCalibraFillHisto(0),
131       fNbTimeBins(0),
132       fNumberBinCharge(50),
133       fSelectTrigger(kTRUE),
134       fSelectedTrigger(new TObjArray()),
135       fRejected(kTRUE),
136       fEsdTrackCuts(0),
137       fRequirePrimaryVertex(kFALSE),
138       fVtxTPC(kFALSE),
139       fVtxSPD(kFALSE),
140       fMinNbContributors(0),
141       fRangePrimaryVertexZ(9999999.0),
142       fMinNbTracks(9),
143       fMaxNbTracks(999999999),
144       fCutWithVdriftCalib(kFALSE),
145       fMinNbTRDtracklets(0),
146       fMinTRDMomentum(0),
147       fScaleGainWithTPCSignal(kFALSE),
148       fLow(0),
149       fHigh(30),
150       fFillZero(kFALSE),
151       fNormalizeNbOfCluster(kFALSE),
152       fRelativeScale(0.0),
153       fMaxCluster(100.0),
154       fNbMaxCluster(2),
155       fOfflineTracks(kFALSE),
156       fStandaloneTracks(kFALSE),
157       fFirstRunGain(-1),
158       fVersionGainUsed(-1),
159       fSubVersionGainUsed(-1),
160       fFirstRunGainLocal(-1),
161       fVersionGainLocalUsed(-1),
162       fSubVersionGainLocalUsed(-1),
163       fFirstRunVdrift(-1),
164       fVersionVdriftUsed(-1), 
165       fSubVersionVdriftUsed(-1),
166       fFirstRunExB(-1),
167       fVersionExBUsed(-1), 
168       fSubVersionExBUsed(-1),
169       fCalDetGain(0x0),
170       fMaxEvent(0),
171       fCounter(0),
172       fPHQon(kTRUE)
173 {
174   //
175   // Default constructor
176   //
177
178   fNz[0] = 0;
179   fNz[1] = 0;
180   fNz[2] = 0;
181   
182   fNrphi[0] = 0;
183   fNrphi[1] = 0;
184   fNrphi[2] = 0;
185
186   
187   // Define input and output slots here
188   // Input slot #0 works with a TChain
189   DefineInput(0, TChain::Class());
190         
191   // Output slot #0 writes into a TList container
192   DefineOutput(1, TList::Class());
193   
194    
195 }
196 //____________________________________________________________________________________
197 AliTRDCalibTask::~AliTRDCalibTask()
198 {
199   //
200   // AliTRDCalibTask destructor
201   //
202
203   // Pointeur
204   if(fNEvents) delete fNEvents;
205   if(fNEventsInput) delete fNEventsInput;
206   if(fNbTRDTrack) delete fNbTRDTrack;
207   if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;
208   if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;
209   if(fNbTPCTRDtrack) delete fNbTPCTRDtrack;
210   if(fNbGoodTracks) delete fNbGoodTracks;
211   if(fNbTimeBin) delete fNbTimeBin;
212   if(fNbTimeBinOffline) delete fNbTimeBinOffline;
213   if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;
214   if(fNbClusters) delete fNbClusters;
215   if(fNbClustersOffline) delete fNbClustersOffline;
216   if(fNbClustersStandalone) delete fNbClustersStandalone;
217   if(fNbTracklets) delete fNbTracklets;
218   if(fNbTrackletsOffline) delete fNbTrackletsOffline;
219   if(fNbTrackletsStandalone) delete fNbTrackletsStandalone;
220   if(fAbsoluteGain) delete fAbsoluteGain;
221   if(fCH2dSum) delete fCH2dSum;
222   if(fPH2dSum) delete fPH2dSum;
223   if(fCH2dSM) delete fCH2dSM;
224   if(fPH2dSM) delete fPH2dSM;
225   if(fCH2dTest) delete fCH2dTest;
226   if(fPH2dTest) delete fPH2dTest;
227   if(fLinearVdriftTest) delete fLinearVdriftTest;
228   if(IsPHQon()){
229     AliTRDdEdxCalibUtils::DeleteHistArray();
230   }
231
232   if(fCalDetGain) delete fCalDetGain;
233   
234   if(fSelectedTrigger) {
235     fSelectedTrigger->Delete();
236     delete fSelectedTrigger;
237   }
238   if(fEsdTrackCuts) {
239     delete fEsdTrackCuts;
240   }
241   
242   if(fTRDChamberStatus) delete fTRDChamberStatus;
243   
244 }
245
246 //________________________________________________________________________
247 void AliTRDCalibTask::UserCreateOutputObjects() 
248 {
249   //
250   // Create the histos
251   //
252   //cout << "AliTRDCalibTask::CreateOutputObjects() IN" << endl;
253
254   // Number of time bins
255   if(fNbTimeBins==0) {
256     AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
257     fNbTimeBins = cal->GetNumberOfTimeBinsDCS();
258     if(fNbTimeBins <= 0){ 
259       AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));
260       fNbTimeBins = 30;
261     }
262   }
263
264   // output list
265   fListHist = new TList();
266   fListHist->SetOwner();
267   
268   // init chamber status
269   fTRDChamberStatus = new AliTRDCalibChamberStatus();
270   fTRDChamberStatus->Init();
271   fListHist->Add(fTRDChamberStatus->GetSparseI());
272   
273   // instance calibration 
274   fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();
275   if(fOnInstance) {
276     fTRDCalibraFillHisto->SetNumberBinCharge(fNumberBinCharge); // set number of bin of the charge distribution
277     fTRDCalibraFillHisto->SetCutWithVdriftCalib(fCutWithVdriftCalib); // cut vdrift
278     fTRDCalibraFillHisto->SetMinNbTRDtracklets(fMinNbTRDtracklets); // min number of TRD tracklets
279     fTRDCalibraFillHisto->SetMinTRDMomentum(fMinTRDMomentum); // min TRD momentum
280     fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms
281     fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors
282     fTRDCalibraFillHisto->SetCH2dOn();  // choose to calibrate the gain
283     fTRDCalibraFillHisto->SetPH2dOn();  // choose to calibrate the drift velocity
284     fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF
285     fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT
286     fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift
287     fTRDCalibraFillHisto->SetExbAltFitOn(fExbAlt); // Alternative method for exb
288     fTRDCalibraFillHisto->SetScaleWithTPCSignal(fScaleGainWithTPCSignal); // Scale Gain with TPC signal
289     for(Int_t k = 0; k < 3; k++){
290       if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {
291         fTRDCalibraFillHisto->SetNz(k,fNz[k]);                                    // Mode calibration
292         fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]);                             // Mode calibration
293       }
294       else {
295         if((fNz[k] == 100) && (fNrphi[k] == 100))  {
296           if(fVector2d) AliInfo("The mode all together is not supported by the vector method");
297           fTRDCalibraFillHisto->SetAllTogether(k);
298         }
299         if((fNz[k] == 10) && (fNrphi[k] == 10))  {
300           if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");
301           fTRDCalibraFillHisto->SetPerSuperModule(k);
302         }
303       }
304     }
305     // Variables for how to fill
306     fTRDCalibraFillHisto->SetFillWithZero(fFillZero);
307     fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fNormalizeNbOfCluster); 
308     fTRDCalibraFillHisto->SetMaxCluster(fMaxCluster);
309     fTRDCalibraFillHisto->SetNbMaxCluster(fNbMaxCluster);
310   
311     // Init with 30 timebins
312     fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos
313     fTRDCalibraFillHisto->SetNumberClusters(fLow); // At least 11 clusters
314     fTRDCalibraFillHisto->SetNumberClustersf(fHigh); // At least 11 clusters
315     
316     // For testing only
317     fTRDCalibraFillHisto->SetDebugLevel(fDebugLevelTRDCalibraFillHisto); //debug stuff
318     
319     if(fHisto2d) {  
320       fListHist->Add(fTRDCalibraFillHisto->GetCH2d());
321       fListHist->Add(fTRDCalibraFillHisto->GetPH2d()); 
322       fListHist->Add(fTRDCalibraFillHisto->GetPRF2d());
323     } 
324     if(fVdriftLinear) fListHist->Add((TObject *)fTRDCalibraFillHisto->GetVdriftLinearFit());
325     if(fVector2d) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector
326     if(fExbAlt) fListHist->Add((TObject *)fTRDCalibraFillHisto->GetExbAltFit());
327   }
328   fRelativeScale = fTRDCalibraFillHisto->GetRelativeScale(); // Get the relative scale for the gain
329   
330   fNEvents = new TH1I(Form("NEvents_%s",(const char*)fName),"NEvents", 2, 0, 2);
331   fListHist->Add(fNEvents);
332   fNEventsInput = new TH1I(Form("NEventsInput_%s",(const char*)fName),"NEventsInput", 2, 0, 2);
333   fListHist->Add(fNEventsInput);
334   
335   // absolute gain calibration even without AliESDfriend
336   Int_t nBinsPt = 25;
337   Double_t minPt = 0.001;
338   Double_t maxPt = 10.0;
339   
340   Double_t *binLimLogPt = new Double_t[nBinsPt+1];
341   Double_t *binLimPt    = new Double_t[nBinsPt+1];
342   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
343   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
344   
345   fAbsoluteGain = new TH2F(Form("AbsoluteGain_%s",(const char*)fName),"AbsoluteGain", 200, 0.0, 700.0, nBinsPt, binLimPt);
346   fAbsoluteGain->SetYTitle("Momentum at TRD");
347   fAbsoluteGain->SetXTitle("charge deposit [a.u]");
348   fAbsoluteGain->SetZTitle("counts");
349   fAbsoluteGain->SetStats(0);
350   fAbsoluteGain->Sumw2();
351   fListHist->Add(fAbsoluteGain);
352   
353   if(IsPHQon()){
354     printf("\n        AliTRDCalibTask PHQ is on!!     \n\n");
355     AliTRDdEdxBaseUtils::PrintControl();
356     AliTRDdEdxCalibUtils::IniHistArray(fListHist, kTRUE);
357   }
358   else{
359     printf("\n        AliTRDCalibTask PHQ is off!!     \n\n");
360   }
361
362   /////////////////////////////////////////
363   // First debug level
364   ///////////////////////////////////////
365   if(fDebug > 0) {
366
367     fLinearVdriftTest = new TH2S(Form("LFDV0testversion_%s",(const char*)fName),"LFDV0testversion",36,-0.9,0.9,48,-1.2,1.2);
368     fLinearVdriftTest->SetXTitle("tan(phi_{track})");
369     fLinearVdriftTest->SetYTitle("dy/dt");
370     fLinearVdriftTest->SetZTitle("Number of tracklets");
371     fLinearVdriftTest->SetStats(0);
372     fLinearVdriftTest->SetDirectory(0);
373     
374     // Standart with AliESDfriend
375     fPH2dTest = new TProfile2D(Form("PH2dTest_%s",(const char*)fName),"Nz0Nrphi0"
376                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
377                            ,540,0,540);
378     fPH2dTest->SetYTitle("Det/pad groups");
379     fPH2dTest->SetXTitle("time [#mus]");
380     fPH2dTest->SetZTitle("<PH> [a.u.]");
381     fPH2dTest->SetStats(0);
382     //
383     fCH2dTest = new TH2I(Form("CH2dTest_%s",(const char*)fName),"Nz0Nrphi0",50,0,300,540,0,540);
384     fCH2dTest->SetYTitle("Det/pad groups");
385     fCH2dTest->SetXTitle("charge deposit [a.u]");
386     fCH2dTest->SetZTitle("counts");
387     fCH2dTest->SetStats(0);
388     fCH2dTest->Sumw2();
389
390     //
391     fPH2dSM = new TProfile2D(Form("PH2dSM_%s",(const char*)fName),"Nz10Nrphi10"
392                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
393                            ,18,0,18);
394     fPH2dSM->SetYTitle("Det/pad groups");
395     fPH2dSM->SetXTitle("time [#mus]");
396     fPH2dSM->SetZTitle("<PH> [a.u.]");
397     fPH2dSM->SetStats(0);
398     //
399     fCH2dSM = new TH2I(Form("CH2dSM_%s",(const char*)fName),"Nz10Nrphi10",50,0,300,18,0,18);
400     fCH2dSM->SetYTitle("Det/pad groups");
401     fCH2dSM->SetXTitle("charge deposit [a.u]");
402     fCH2dSM->SetZTitle("counts");
403     fCH2dSM->SetStats(0);
404     fCH2dSM->Sumw2();
405     //
406     fPH2dSum = new TProfile2D(Form("PH2dSum_%s",(const char*)fName),"Nz100Nrphi100"
407                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
408                             ,1,0,1);
409     fPH2dSum->SetYTitle("Det/pad groups");
410     fPH2dSum->SetXTitle("time [#mus]");
411     fPH2dSum->SetZTitle("<PH> [a.u.]");
412     fPH2dSum->SetStats(0);
413     //
414     fCH2dSum = new TH2I(Form("CH2dSum_%s",(const char*)fName),"Nz100Nrphi100",50,0,300,1,0,1);
415     fCH2dSum->SetYTitle("Det/pad groups");
416     fCH2dSum->SetXTitle("charge deposit [a.u]");
417     fCH2dSum->SetZTitle("counts");
418     fCH2dSum->SetStats(0);
419     fCH2dSum->Sumw2();
420     
421     
422     // Add them
423     fListHist->Add(fLinearVdriftTest);
424     fListHist->Add(fPH2dTest);
425     fListHist->Add(fCH2dTest);
426     fListHist->Add(fPH2dSM);
427     fListHist->Add(fCH2dSM);
428     fListHist->Add(fPH2dSum);
429     fListHist->Add(fCH2dSum);
430
431   }
432
433   /////////////////////////////////////////
434   // Second debug level
435   ///////////////////////////////////////
436   if(fDebug > 1) {
437
438     fNbGoodTracks = new TH2F(Form("NbGoodTracks_%s",(const char*)fName),"NbGoodTracks",500,0.0,2500.0,200,0.0,100.0);
439     fNbGoodTracks->SetXTitle("Nb of good tracks");
440     fNbGoodTracks->SetYTitle("Centrality");
441     fNbGoodTracks->SetStats(0);
442
443     fNbTRDTrack = new TH1F(Form("TRDTrack_%s",(const char*)fName),"TRDTrack",50,0,50);
444     fNbTRDTrack->Sumw2();
445     fNbTRDTrackOffline = new TH1F(Form("TRDTrackOffline_%s",(const char*)fName),"TRDTrackOffline",50,0,50);
446     fNbTRDTrackOffline->Sumw2();
447     fNbTRDTrackStandalone = new TH1F(Form("TRDTrackStandalone_%s",(const char*)fName),"TRDTrackStandalone",50,0,50);
448     fNbTRDTrackStandalone->Sumw2();
449     fNbTPCTRDtrack = new TH2F(Form("NbTPCTRDtrack_%s",(const char*)fName),"NbTPCTRDtrack",100,0,100,100,0,100);
450     fNbTPCTRDtrack->Sumw2();
451     //
452     fNbTimeBin = new TH1F(Form("NbTimeBin_%s",(const char*)fName),"NbTimeBin",35,0,35);
453     fNbTimeBin->Sumw2();
454     fNbTimeBinOffline = new TH1F(Form("NbTimeBinOffline_%s",(const char*)fName),"NbTimeBinOffline",35,0,35);
455     fNbTimeBinOffline->Sumw2();
456     fNbTimeBinStandalone = new TH1F(Form("NbTimeBinStandalone_%s",(const char*)fName),"NbTimeBinStandalone",35,0,35);
457     fNbTimeBinStandalone->Sumw2();
458     //
459     fNbClusters = new TH1F(Form("NbClusters_%s",(const char*)fName),"",35,0,35);
460     fNbClusters->Sumw2();
461     fNbClustersOffline = new TH1F(Form("NbClustersOffline_%s",(const char*)fName),"",35,0,35);
462     fNbClustersOffline->Sumw2();
463     fNbClustersStandalone = new TH1F(Form("NbClustersStandalone_%s",(const char*)fName),"",35,0,35);
464     fNbClustersStandalone->Sumw2();
465     //
466     fNbTracklets = new TH1F(Form("NbTracklets_%s",(const char*)fName),"NbTracklets",540,0.,540.);
467     fNbTracklets->Sumw2();
468     fNbTrackletsOffline = new TH1F(Form("NbTrackletsOffline_%s",(const char*)fName),"NbTrackletsOffline",540,0.,540.);
469     fNbTrackletsOffline->Sumw2();
470     fNbTrackletsStandalone = new TH1F(Form("NbTrackletsStandalone_%s",(const char*)fName),"NbTrackletsStandalone",540,0.,540.);
471     fNbTrackletsStandalone->Sumw2();
472    
473     fListHist->Add(fNbGoodTracks);
474    
475     fListHist->Add(fNbTRDTrack);
476     fListHist->Add(fNbTRDTrackOffline);
477     fListHist->Add(fNbTRDTrackStandalone);
478     fListHist->Add(fNbTPCTRDtrack);
479     
480     fListHist->Add(fNbTimeBin);
481     fListHist->Add(fNbTimeBinOffline);
482     fListHist->Add(fNbTimeBinStandalone);
483     fListHist->Add(fNbClusters);
484     fListHist->Add(fNbClustersOffline);
485     fListHist->Add(fNbClustersStandalone);
486     fListHist->Add(fNbTracklets);
487     fListHist->Add(fNbTrackletsOffline);
488     fListHist->Add(fNbTrackletsStandalone);
489     
490   }
491
492   delete [] binLimLogPt;
493   delete [] binLimPt;
494
495   PostData(1,fListHist);
496
497   //cout << "AliTRDCalibTask::UserCreateOutputObjects() OUT" << endl;
498
499 }
500
501 //________________________________________________________________________
502 void AliTRDCalibTask::UserExec(Option_t *) 
503 {
504   //
505   // Filling of the histos
506   //
507   //cout << "AliTRDCalibTask::Exec() IN" << endl;
508   
509   // Init Versions and subversions used
510   if((fFirstRunGain==-1) || (fVersionGainUsed==-1) || (fSubVersionGainUsed==-1) || (fFirstRunGainLocal==-1) || (fVersionGainLocalUsed==-1) || (fSubVersionGainLocalUsed==-1) || (fFirstRunVdrift==-1) || (fVersionVdriftUsed==-1) || (fSubVersionVdriftUsed==-1)) {
511     if(!SetVersionSubversion()) {
512       PostData(1, fListHist);
513       return;
514     }
515   }
516   if(fCounter==0) {
517     if(fOnInstance) {
518       fTRDCalibraFillHisto->SetFirstRunGain(fFirstRunGain); // Gain Used
519       fTRDCalibraFillHisto->SetVersionGainUsed(fVersionGainUsed); // Gain Used
520       fTRDCalibraFillHisto->SetSubVersionGainUsed(fSubVersionGainUsed); // Gain Used
521       fTRDCalibraFillHisto->SetFirstRunGainLocal(fFirstRunGainLocal); // Gain Used
522       fTRDCalibraFillHisto->SetVersionGainLocalUsed(fVersionGainLocalUsed); // Gain Used
523       fTRDCalibraFillHisto->SetSubVersionGainLocalUsed(fSubVersionGainLocalUsed); // Gain Used
524       fTRDCalibraFillHisto->SetFirstRunVdrift(fFirstRunVdrift); // Vdrift Used
525       fTRDCalibraFillHisto->SetVersionVdriftUsed(fVersionVdriftUsed); // Vdrift Used
526       fTRDCalibraFillHisto->SetSubVersionVdriftUsed(fSubVersionVdriftUsed); // Vdrift Used
527       if((fFirstRunExB != -1) && (fVersionExBUsed != -1) && (fSubVersionExBUsed != -1)){
528         fTRDCalibraFillHisto->SetFirstRunExB(fFirstRunExB); // ExB Used
529         fTRDCalibraFillHisto->SetVersionExBUsed(fVersionExBUsed); // ExB Used
530         fTRDCalibraFillHisto->SetSubVersionExBUsed(fSubVersionExBUsed); // ExB Used
531       }
532       fTRDCalibraFillHisto->InitCalDet();
533     }
534     if(fDebug > 1){
535       // title CH2dTest
536       TString name("Ver");
537       name += fVersionGainUsed;
538       name += "Subver";
539       name += fSubVersionGainUsed;
540       name += "FirstRun";
541       name += fFirstRunGain;
542       name += "Nz0Nrphi0";
543       fCH2dTest->SetTitle(name);  
544       // title PH2dTest
545       TString namee("Ver");
546       namee += fVersionVdriftUsed;
547       namee += "Subver";
548       namee += fSubVersionVdriftUsed;
549       namee += "FirstRun";
550       namee += fFirstRunVdrift;
551       namee += "Nz0Nrphi0";
552       fPH2dTest->SetTitle(namee); 
553     }
554   }
555   
556   //  AliLog::SetGlobalLogLevel(AliLog::kError);
557   //  cout << "AliTRDCalibTask::Exec() 1" << endl;
558   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
559   if(!fESD){
560     AliError("ESD Event missing");
561     PostData(1, fListHist);
562     return;
563   }
564
565   const char* type = fESD->GetBeamType();
566   
567   
568   //printf("Counter %d\n",fCounter);
569   
570   fCounter++;
571   fNEventsInput->Fill(1);
572
573   //cout << "maxEvent = " << fMaxEvent << endl;
574   //if(fCounter%100==0) cout << "fCounter = " << fCounter << endl;
575   if((fMaxEvent != 0) && (fMaxEvent < fCounter)) {
576     PostData(1, fListHist);
577     return;
578   }
579   //if(fCounter%100==0) cout << "fCounter1 = " << fCounter << endl;
580   //cout << "event = " << fCounter << endl;
581   
582   //printf("Counter %d\n",fCounter);
583
584   ///////////////////
585   // Check trigger
586   ///////////////////
587   Bool_t pass = kTRUE;
588
589   if (fSelectTrigger) {
590    
591     //printf("Will check the triggers\n");
592
593     Int_t numberOfTriggerSelected = fSelectedTrigger->GetEntriesFast();
594     //printf("numberofTriggerSelected %d\n",numberOfTriggerSelected);
595     if(fRejected) {
596       pass = kTRUE;
597       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
598         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
599         const TString tString=obString->GetString();
600         if(fESD->IsTriggerClassFired((const char*)tString)) {
601           pass = kFALSE;
602         }
603       }
604     }
605     else {
606       pass = kFALSE;
607       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
608         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
609         const TString tString=obString->GetString();
610         if(fESD->IsTriggerClassFired((const char*)tString)) {
611           pass = kTRUE;
612         }
613       }
614     }
615     if(!pass) {
616       PostData(1, fListHist);
617       return;
618     }   
619
620   }
621     
622   //printf("Class Fired %s\n",(const char*)fESD->GetFiredTriggerClasses());
623   //printf("Trigger passed\n");
624   
625   ///////////////////////////////
626   // Require a primary vertex
627   //////////////////////////////
628   if(fRequirePrimaryVertex) {
629     const AliESDVertex* vtxESD = 0x0;
630     if      (fVtxTPC) vtxESD = fESD->GetPrimaryVertexTPC() ;
631     else if (fVtxSPD) vtxESD = fESD->GetPrimaryVertexSPD() ;
632     else              vtxESD = fESD->GetPrimaryVertexTracks() ;
633     if(!vtxESD){
634       PostData(1, fListHist);
635       return;
636     }
637     Int_t nCtrb = vtxESD->GetNContributors();
638     if(nCtrb < fMinNbContributors) {
639       PostData(1, fListHist);     
640       return;
641     }
642     Double_t zPosition = vtxESD->GetZ();
643     if(TMath::Abs(zPosition) > fRangePrimaryVertexZ) {
644       PostData(1, fListHist);
645       return;
646     }     
647     
648   }
649   
650   //printf("Primary vertex passed\n");
651   
652   //////////////////////////////////////
653   // Requirement on number of good tracks
654   //////////////////////////////////////
655   Int_t nGoodParticles = 0;
656   Double_t nbTracks = fESD->GetNumberOfTracks();
657   for(Int_t itrack = 0; itrack < nbTracks; itrack++) {
658     if(ParticleGood(itrack)) nGoodParticles++;  
659   }
660   if(fDebug > 1)  {
661     // Centrality
662     AliCentrality *esdCentrality = fESD->GetCentrality();
663     Float_t centrality = esdCentrality->GetCentralityPercentile("V0M");
664     //Float_t centralityb = esdCentrality->GetCentralityPercentile("CL1");
665     fNbGoodTracks->Fill(nGoodParticles,centrality);
666     //printf("centrality %f, centralityb %f\n",centrality,centralityb);
667   }
668   
669   //printf("Beam type %s\n",(const char*)type);
670   if (!strstr(type,"p-p")) {
671     //if (strstr(type,"A-A")) {
672     //printf("Will check the number of good tracks %f %f %f\n",nGoodParticles,fMinNbTracks,fMaxNbTracks);
673     if((nGoodParticles < fMinNbTracks) || (nGoodParticles > fMaxNbTracks)) {
674       PostData(1, fListHist);
675       return;
676     }
677   }
678   
679   fNEvents->Fill(1);
680   
681   // In total
682   Int_t nbTrdTracks = 0;
683   // standalone
684   Int_t nbTrdTracksStandalone = 0;
685   // offline
686   Int_t nbTrdTracksOffline = 0;
687   // TPC
688   Int_t nbtrackTPC = 0;
689   
690
691
692   if (nbTracks <= 0.0) {
693     
694     if(fDebug > 1) {
695       fNbTRDTrack->Fill(nbTrdTracks);
696       fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
697       fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
698     }
699     PostData(1, fListHist);
700     return;
701   }
702
703
704   fESDfriend = dynamic_cast<AliESDfriend*> (fESD->FindListObject("AliESDfriend"));
705   if(!fESDfriend){
706     AliError("fESDfriend not available");
707     PostData(1, fListHist);
708     return;
709   }
710
711   if(fESDfriend->TestSkipBit()) {
712     PostData(1, fListHist);
713     return;
714   }
715
716   //printf("has friends\n");
717
718   /////////////////////////////////////
719   // Loop on AliESDtrack
720   ////////////////////////////////////
721   //printf("Nb of tracks %f\n",nbTracks);      
722   for(int itrk=0; itrk < nbTracks; ++itrk){
723     
724     // Get ESD track
725     fkEsdTrack = fESD->GetTrack(itrk);
726     if(!fkEsdTrack) continue;
727     ULong_t status = fkEsdTrack->GetStatus(); 
728     if(status&(AliESDtrack::kTPCout)) ++nbtrackTPC;
729     
730     fFriendTrack = fESDfriend->GetTrack(itrk);
731     if(!fFriendTrack)  {
732       //printf("No friend track %d\n",itrk);
733       continue;
734     }
735
736     // Other cuts
737     fTrdTrack = 0x0;
738     Bool_t good = kTRUE;
739     Bool_t standalonetrack = kFALSE;
740     Bool_t offlinetrack = kFALSE;
741     //ULong_t status = fkEsdTrack->GetStatus();
742     
743     //////////////////////////////////////
744     // Loop on calibration objects
745     //////////////////////////////////////
746     Int_t icalib=0;
747     Int_t nTRDtrackV1=0;
748     while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
749       //printf("Name %s\n",fCalibObject->IsA()->GetName());
750       if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
751       //printf("Find the calibration object\n");
752       ++nTRDtrackV1;
753       
754       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
755         standalonetrack = kTRUE;
756       }
757       if((status&(AliESDtrack::kTRDin))) {
758         offlinetrack = kTRUE;
759       }
760       if(fOfflineTracks){
761         if(!offlinetrack){
762           good = kFALSE;
763         }
764       }
765       else if(fStandaloneTracks){
766         if(!standalonetrack){
767           good = kFALSE;
768         }
769       }
770       
771       fTrdTrack = (AliTRDtrackV1 *)fCalibObject;
772       // process chamberstatus
773       fTRDChamberStatus->ProcessTrack(fTrdTrack);
774     }
775
776     // Quality cuts on the AliESDtrack
777     if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {
778       //printf("Not a good track\n");
779       continue;
780     }
781
782     // First Absolute gain calibration
783     Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
784     Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID(); 
785     //printf("Number of trd tracklets %d and PID trd tracklets %d\n",trdNTracklets,trdNTrackletsPID);
786     if((trdNTracklets > 0) && (trdNTrackletsPID > 0)) {
787       for(Int_t iPlane = 0; iPlane < 6; ++iPlane){
788         //Double_t slide = fkEsdTrack->GetTRDslice(iPlane);
789         //printf("Number of slide %d\n",fkEsdTrack->GetNumberOfTRDslices());
790         //Double_t momentum = fkEsdTrack->GetTRDmomentum(iPlane);
791         //printf("momentum %f, slide %f\n",momentum,slide);
792         if(fkEsdTrack->GetTRDslice(iPlane) > 0.0) 
793           fAbsoluteGain->Fill(fkEsdTrack->GetTRDslice(iPlane)*8.0/100.0,
794                               fkEsdTrack->GetTRDmomentum(iPlane)); 
795       }
796     }     
797     
798     
799     if(!fTrdTrack) continue; 
800
801     if(good && fOnInstance) {
802       //cout << "good" << endl;
803       fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack,fkEsdTrack);
804       //printf("Fill fTRDCalibraFillHisto\n");
805     }
806
807     if(IsPHQon()){
808       const Double_t mag = AliTRDdEdxBaseUtils::IsExBOn() ? fESD->GetMagneticField() : -1;
809       const Int_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? fkEsdTrack->Charge() : -1;
810       const Double_t toTPCscale = AliTRDdEdxCalibUtils::GetCalibTPCscale(fkEsdTrack->GetTPCncls(), fkEsdTrack->GetTPCsignal());
811       if(toTPCscale>0){
812         AliTRDdEdxCalibUtils::FillHist(fTrdTrack, 0, mag, charge, toTPCscale);
813       }
814     }
815
816     //////////////////////////////////
817     // Debug 
818     ////////////////////////////////
819     
820     if(fDebug > 0) {
821       
822       //printf("Enter debug\n");
823       
824       Int_t nbtracklets = 0;
825       
826       // Check some stuff
827       Bool_t standalonetracklet = kFALSE;  
828       const AliTRDseedV1 *tracklet = 0x0;
829       //////////////////////////////////////
830       // Loop tracklets
831       ///////////////////////////////////// 
832       Int_t nbclusters=0;
833       Double_t phtb[AliTRDseedV1::kNtb];
834       memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
835       Double_t sum = 0.0;
836       Float_t normalisation = 1.13;
837       Int_t detector = 0;
838       Int_t sector = 0;
839       for(Int_t itr = 0; itr < 6; ++itr){
840         
841         if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;
842         if(!tracklet->IsOK()) continue;
843         ++nbtracklets;
844         standalonetracklet = kFALSE; 
845         if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
846         
847         nbclusters = 0;
848         memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
849         sum = 0.0;
850         detector = 0;
851         sector = 0;
852         //Int_t crossrow = 0;
853         
854           // Check no shared clusters
855           //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
856           //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
857           // }
858         
859           // Loop on clusters
860         Int_t time = 0;
861         Float_t ch = 0;
862         Float_t qcl = 0;
863         for(int ic=0; ic<AliTRDseedV1::kNtb; ++ic){
864           
865           if(!(fCl = tracklet->GetClusters(ic))) continue;
866           ++nbclusters;
867           time = fCl->GetPadTime();
868           //ch =  tracklet->GetdQdl(ic);
869           ch =  tracklet->GetQperTB(ic);
870           qcl = TMath::Abs(fCl->GetQ());
871           detector = fCl->GetDetector();          
872           // Add the charge if shared cluster
873           if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
874             if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {
875               qcl += TMath::Abs(fCl->GetQ());
876               //printf("Add the cluster charge\n");
877             }
878           }
879           if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
880           if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation; 
881           else sum += ch/normalisation;
882           
883           if(fDebug > 1) {
884             fNbTimeBin->Fill(time);
885             if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
886             else fNbTimeBinOffline->Fill(time);
887           }
888         }
889         sector = AliTRDgeometry::GetSector(detector);
890         
891         if(fDebug > 1) {
892           fNbTracklets->Fill(detector);
893           if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);
894           else fNbTrackletsOffline->Fill(detector);
895           
896           fNbClusters->Fill(nbclusters);
897           if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
898           else  fNbClustersOffline->Fill(nbclusters);
899         }          
900         
901         if((nbclusters > fLow) && (nbclusters < fHigh)){
902           if(fRelativeScale > 0.0) sum = sum/fRelativeScale;
903           fCH2dTest->Fill(sum,detector+0.5);           
904           fCH2dSM->Fill(sum,sector+0.5);
905           fCH2dSum->Fill(sum,0.5);
906           Bool_t checknoise = kTRUE;
907           if(fMaxCluster > 0) {
908             if(phtb[0] > fMaxCluster) checknoise = kFALSE;
909             if(fNbTimeBins > fNbMaxCluster) {
910               for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){
911                 if(phtb[k] > fMaxCluster) checknoise = kFALSE;
912               }
913             }
914           }
915           if(checknoise) {             
916             for(int ic=0; ic<fNbTimeBins; ic++){
917               if(fFillZero) {
918                 fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
919                 fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
920                 fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
921               }
922               else {
923                 if(phtb[ic] > 0.0) {
924                   fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
925                   fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
926                   fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
927                 }
928               }
929             }
930           }
931         }
932         if(detector == 0) FindP1TrackPHtrackletV1Test(tracklet,nbclusters);
933         
934       } // loop on tracklets
935       
936       
937     } // debug
938     
939     if(nTRDtrackV1 > 0) {
940       ++nbTrdTracks;      
941       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
942         ++nbTrdTracksStandalone;
943       }
944       if((status&(AliESDtrack::kTRDin))) {
945         ++nbTrdTracksOffline;
946       }
947     }
948     //delete fFriendTrack;
949   } // loop ESD track
950   
951   if(fDebug > 1) {
952     fNbTRDTrack->Fill(nbTrdTracks);
953     fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
954     fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
955     fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
956   }
957   
958   // Post output data
959   PostData(1, fListHist);
960   //cout << "AliTRDCalibTask::Exec() OUT" << endl;
961 }
962      
963 //________________________________________________________________________
964 void AliTRDCalibTask::Terminate(Option_t *) 
965 {
966   //
967   // Terminate
968   //
969   
970   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
971
972  
973 }
974 //_______________________________________________________
975 Bool_t AliTRDCalibTask::Load(const Char_t *filename)
976 {
977   //
978   // Generic container loader
979   //
980
981   if(!TFile::Open(filename)){
982     //AliWarning(Form("Couldn't open file %s.", filename));
983     return kFALSE;
984   }
985   TList *o = 0x0;
986   if(!(o = (TList*)gFile->Get(GetName()))){
987     //AliWarning("Missing histogram container.");
988     return kFALSE;
989   }
990   fListHist = (TList*)o->Clone(GetName());
991   gFile->Close();
992   return kTRUE;
993 }
994 //_______________________________________________________
995 Bool_t AliTRDCalibTask::Load(TList *lister)
996 {
997   //
998   // Generic container loader
999   //
1000
1001   fListHist = (TList*)lister->Clone(GetName());
1002   return kTRUE;
1003 }
1004 //_______________________________________________________________________________________
1005 void  AliTRDCalibTask::AddTask(const AliTRDCalibTask * calibTask) {
1006
1007   //
1008   // Add stats
1009   //
1010
1011   TList *listcalibTask = calibTask->GetList();
1012   if(!listcalibTask) return;
1013
1014   THnSparseI *histoEntries = (THnSparseI *) listcalibTask->FindObject("NumberOfEntries");
1015
1016   TH1I *nEvents  = (TH1I *) listcalibTask->FindObject(Form("NEvents_%s",(const char*)calibTask->GetName()));
1017   TH1I *nEventsInput  = (TH1I *) listcalibTask->FindObject(Form("NEventsInput_%s",(const char*)calibTask->GetName()));
1018   TH2F *absoluteGain  = (TH2F *) listcalibTask->FindObject(Form("AbsoluteGain_%s",(const char*)calibTask->GetName()));
1019
1020   TH1F *trdTrack = (TH1F *) listcalibTask->FindObject(Form("TRDTrack_%s",(const char*)calibTask->GetName()));
1021   TH1F *trdTrackOffline = (TH1F *) listcalibTask->FindObject(Form("TRDTrackOffline_%s",(const char*)calibTask->GetName()));
1022   TH1F *trdTrackStandalone = (TH1F *) listcalibTask->FindObject(Form("TRDTrackStandalone_%s",(const char*)calibTask->GetName()));
1023
1024   TH2F *tpctrdTrack = (TH2F *) listcalibTask->FindObject(Form("NbTPCTRDtrack_%s",(const char*)calibTask->GetName()));
1025
1026   TH1F *nbTimeBin = (TH1F *) listcalibTask->FindObject(Form("NbTimeBin_%s",(const char*)calibTask->GetName()));
1027   TH1F *nbTimeBinOffline = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinOffline_%s",(const char*)calibTask->GetName()));
1028   TH1F *nbTimeBinStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinStandalone_%s",(const char*)calibTask->GetName()));
1029
1030   TH1F *nbClusters = (TH1F *) listcalibTask->FindObject(Form("NbClusters_%s",(const char*)calibTask->GetName()));
1031   TH1F *nbClustersOffline = (TH1F *) listcalibTask->FindObject(Form("NbClustersOffline_%s",(const char*)calibTask->GetName()));
1032   TH1F *nbClustersStandalone = (TH1F *) listcalibTask->FindObject(Form("NbClustersStandalone_%s",(const char*)calibTask->GetName()));
1033
1034   TH1F *nbTracklets = (TH1F *) listcalibTask->FindObject(Form("NbTracklets_%s",(const char*)calibTask->GetName()));
1035   TH1F *nbTrackletsOffline = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsOffline_%s",(const char*)calibTask->GetName()));
1036   TH1F *nbTrackletsStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsStandalone_%s",(const char*)calibTask->GetName()));
1037   
1038   TH2I *ch2d = (TH2I *) listcalibTask->FindObject("CH2d");
1039   TProfile2D *ph2d = (TProfile2D *) listcalibTask->FindObject("PH2d");
1040   TProfile2D *prf2d = (TProfile2D *) listcalibTask->FindObject("PRF2d");
1041
1042   TH2I *ch2dSum = (TH2I *) listcalibTask->FindObject(Form("CH2dSum_%s",(const char*)calibTask->GetName()));
1043   TProfile2D *ph2dSum = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSum_%s",(const char*)calibTask->GetName()));
1044
1045   TH2I *ch2dSM = (TH2I *) listcalibTask->FindObject(Form("CH2dSM_%s",(const char*)calibTask->GetName()));
1046   TProfile2D *ph2dSM = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSM_%s",(const char*)calibTask->GetName()));
1047   
1048   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) listcalibTask->FindObject("AliTRDCalibraVdriftLinearFit");  
1049   AliTRDCalibraExbAltFit *exbaltfit = (AliTRDCalibraExbAltFit *) listcalibTask->FindObject("AliTRDCalibraExbAltFit");  
1050   AliTRDCalibraVector *calibraVector = (AliTRDCalibraVector *) listcalibTask->FindObject("AliTRDCalibraVector");  
1051
1052   //
1053
1054   THnSparseI *inhistoEntries = (THnSparseI *) fListHist->FindObject("NumberOfEntries");
1055
1056   TH1I *inEventsInput  = (TH1I *) fListHist->FindObject(Form("NEventsInput_%s",(const char*)fName));
1057   TH1I *inEvents  = (TH1I *) fListHist->FindObject(Form("NEvents_%s",(const char*)fName));
1058   TH2F *iabsoluteGain  = (TH2F *) fListHist->FindObject(Form("AbsoluteGain_%s",(const char*)fName));
1059
1060   TH1F *itrdTrack = (TH1F *) fListHist->FindObject(Form("TRDTrack_%s",(const char*)fName));
1061   TH1F *itrdTrackOffline = (TH1F *) fListHist->FindObject(Form("TRDTrackOffline_%s",(const char*)fName));
1062   TH1F *itrdTrackStandalone = (TH1F *) fListHist->FindObject(Form("TRDTrackStandalone_%s",(const char*)fName));
1063
1064   TH2F *itpctrdTrack = (TH2F *) fListHist->FindObject(Form("NbTPCTRDtrack_%s",(const char*)fName));
1065
1066   TH1F *inbTimeBin = (TH1F *) fListHist->FindObject(Form("NbTimeBin_%s",(const char*)fName));
1067   TH1F *inbTimeBinOffline = (TH1F *) fListHist->FindObject(Form("NbTimeBinOffline_%s",(const char*)fName));
1068   TH1F *inbTimeBinStandalone = (TH1F *) fListHist->FindObject(Form("NbTimeBinStandalone_%s",(const char*)fName));
1069
1070   TH1F *inbClusters = (TH1F *) fListHist->FindObject(Form("NbClusters_%s",(const char*)fName));
1071   TH1F *inbClustersOffline = (TH1F *) fListHist->FindObject(Form("NbClustersOffline_%s",(const char*)fName));
1072   TH1F *inbClustersStandalone = (TH1F *) fListHist->FindObject(Form("NbClustersStandalone_%s",(const char*)fName));
1073
1074   TH1F *inbTracklets = (TH1F *) fListHist->FindObject(Form("NbTracklets_%s",(const char*)fName));
1075   TH1F *inbTrackletsOffline = (TH1F *) fListHist->FindObject(Form("NbTrackletsOffline_%s",(const char*)fName));
1076   TH1F *inbTrackletsStandalone = (TH1F *) fListHist->FindObject(Form("NbTrackletsStandalone_%s",(const char*)fName));
1077   
1078   TH2I *ich2d = (TH2I *) fListHist->FindObject("CH2d");
1079   TProfile2D *iph2d = (TProfile2D *) fListHist->FindObject("PH2d");
1080   TProfile2D *iprf2d = (TProfile2D *) fListHist->FindObject("PRF2d");
1081
1082   TH2I *ich2dSum = (TH2I *) fListHist->FindObject(Form("CH2dSum_%s",(const char*)fName));
1083   TProfile2D *iph2dSum = (TProfile2D *) fListHist->FindObject(Form("PH2dSum_%s",(const char*)fName));
1084
1085   TH2I *ich2dSM = (TH2I *) fListHist->FindObject(Form("CH2dSM_%s",(const char*)fName));
1086   TProfile2D *iph2dSM = (TProfile2D *) fListHist->FindObject(Form("PH2dSM_%s",(const char*)fName));
1087   
1088   AliTRDCalibraVdriftLinearFit *ilinearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");  
1089   AliTRDCalibraExbAltFit *iexbaltfit = (AliTRDCalibraExbAltFit *) fListHist->FindObject("AliTRDCalibraExbAltFit");
1090   AliTRDCalibraVector *icalibraVector = (AliTRDCalibraVector *) fListHist->FindObject("AliTRDCalibraVector");  
1091
1092
1093   // Add
1094
1095   if(histoEntries) {
1096     if(inhistoEntries) {
1097       inhistoEntries->Add(histoEntries);
1098       //printf("Add Events\n");
1099     }
1100     else {
1101       //printf("Create new Events\n");
1102       inhistoEntries = (THnSparseI *) histoEntries->Clone();
1103       fListHist->Add(inhistoEntries);
1104     }
1105   }
1106
1107   if(nEventsInput) {
1108     if(inEventsInput) {
1109       inEventsInput->Add(nEventsInput);
1110       //printf("Add Events\n");
1111     }
1112     else {
1113       //printf("Create new Events\n");
1114       inEventsInput = new TH1I(*nEventsInput);
1115       fListHist->Add(inEventsInput);
1116     }
1117   }
1118   
1119   if(nEvents) {
1120     if(inEvents) {
1121       inEvents->Add(nEvents);
1122       //printf("Add Events\n");
1123     }
1124     else {
1125       //printf("Create new Events\n");
1126       inEvents = new TH1I(*nEvents);
1127       fListHist->Add(inEvents);
1128     }
1129   }
1130   
1131   if(absoluteGain) {
1132     if(iabsoluteGain) iabsoluteGain->Add(absoluteGain);
1133     else {
1134       iabsoluteGain = new TH2F(*absoluteGain);
1135       fListHist->Add(iabsoluteGain);
1136     }
1137   }
1138   
1139   if(trdTrack) {
1140     if(itrdTrack) itrdTrack->Add(trdTrack);
1141     else {
1142      itrdTrack = new TH1F(*trdTrack);
1143      fListHist->Add(itrdTrack);
1144     }
1145   }
1146
1147   if(trdTrackOffline) {
1148     if(itrdTrackOffline) itrdTrackOffline->Add(trdTrackOffline);
1149     else {
1150       itrdTrackOffline = new TH1F(*trdTrackOffline);
1151       fListHist->Add(itrdTrackOffline);
1152     }
1153   }
1154
1155   if(trdTrackStandalone) {
1156     if(itrdTrackStandalone) itrdTrackStandalone->Add(trdTrackStandalone);
1157     else {
1158       itrdTrackStandalone = new TH1F(*trdTrackStandalone);
1159       fListHist->Add(itrdTrackStandalone);
1160     }
1161   }
1162
1163   if(tpctrdTrack) {
1164     if(itpctrdTrack) itpctrdTrack->Add(tpctrdTrack);
1165     else {
1166       itpctrdTrack = new TH2F(*tpctrdTrack);
1167       fListHist->Add(itpctrdTrack);
1168     }
1169   }
1170
1171   if(nbTimeBin) {
1172     if(inbTimeBin) inbTimeBin->Add(nbTimeBin);
1173     else {
1174       inbTimeBin = new TH1F(*inbTimeBin);
1175       fListHist->Add(inbTimeBin);
1176     }
1177   }
1178
1179   if(nbTimeBinOffline) {
1180     if(inbTimeBinOffline) inbTimeBinOffline->Add(nbTimeBinOffline);
1181     else {
1182       inbTimeBinOffline = new TH1F(*nbTimeBinOffline);
1183       fListHist->Add(inbTimeBinOffline);
1184     }
1185   }
1186   
1187   if(nbTimeBinStandalone) {
1188     if(inbTimeBinStandalone) inbTimeBinStandalone->Add(nbTimeBinStandalone);
1189     else {
1190       inbTimeBinStandalone = new TH1F(*nbTimeBinStandalone);
1191       fListHist->Add(inbTimeBinStandalone);
1192     }
1193   }
1194
1195   if(nbClusters) {
1196     if(inbClusters) inbClusters->Add(nbClusters);
1197     else {
1198       inbClusters = new TH1F(*nbClusters);
1199       fListHist->Add(inbClusters);
1200     }
1201   }
1202   
1203   if(nbClustersOffline) {
1204     if(inbClustersOffline) inbClustersOffline->Add(nbClustersOffline);
1205     else {
1206       inbClustersOffline = new TH1F(*nbClustersOffline);
1207       fListHist->Add(inbClustersOffline);
1208     }
1209   }
1210   
1211   if(nbClustersStandalone) {
1212     if(inbClustersStandalone) inbClustersStandalone->Add(nbClustersStandalone);
1213     else {
1214       inbClustersStandalone = new TH1F(*nbClustersStandalone);
1215       fListHist->Add(inbClustersStandalone);
1216     }
1217   }
1218
1219   if(nbTracklets) {
1220     if(inbTracklets) inbTracklets->Add(nbTracklets);
1221     else {
1222       inbTracklets = new TH1F(*nbTracklets);
1223       fListHist->Add(inbTracklets);
1224     }
1225   }
1226
1227   if(nbTrackletsOffline) {
1228     if(inbTrackletsOffline) inbTrackletsOffline->Add(nbTrackletsOffline);
1229     else {
1230       inbTrackletsOffline = new TH1F(*nbTrackletsOffline);
1231       fListHist->Add(inbTrackletsOffline);
1232     }
1233   }
1234   
1235   if(nbTrackletsStandalone) {
1236     if(inbTrackletsStandalone) inbTrackletsStandalone->Add(nbTrackletsStandalone);
1237     else {
1238       inbTrackletsStandalone = new TH1F(*nbTrackletsStandalone);
1239       fListHist->Add(inbTrackletsStandalone);
1240     }
1241   }
1242   
1243   if(ch2d) {
1244     if(ich2d) ich2d->Add(ch2d);
1245     else {
1246       ich2d = new TH2I(*ch2d);
1247       fListHist->Add(ich2d);
1248     }
1249   }
1250
1251   if(ph2d) {
1252     if(iph2d) iph2d->Add(ph2d);
1253     else {
1254       iph2d = new TProfile2D(*ph2d);
1255       fListHist->Add(iph2d);
1256     }
1257   }
1258
1259   if(prf2d) {
1260     if(iprf2d) iprf2d->Add(prf2d);
1261     else {
1262       iprf2d = new TProfile2D(*prf2d);
1263       fListHist->Add(iprf2d);
1264     }
1265   }
1266
1267   if(ch2dSum) {
1268     if(ich2dSum) ich2dSum->Add(ch2dSum);
1269     else {
1270       ich2dSum = new TH2I(*ch2dSum);
1271       fListHist->Add(ich2dSum);
1272     }
1273   }
1274
1275   if(ph2dSum) {
1276     if(iph2dSum) iph2dSum->Add(ph2dSum);
1277     else {
1278       iph2dSum = new TProfile2D(*ph2dSum);
1279       fListHist->Add(iph2dSum);
1280     }
1281   }
1282
1283   if(ch2dSM) {
1284     if(ich2dSM) ich2dSM->Add(ch2dSM);
1285     else {
1286       ich2dSM = new TH2I(*ch2dSM);
1287       fListHist->Add(ich2dSM);
1288     }
1289   }
1290
1291   if(ph2dSM) {
1292     if(iph2dSM) iph2dSM->Add(ph2dSM);
1293     else {
1294       iph2dSM = new TProfile2D(*ph2dSM);
1295       fListHist->Add(iph2dSM);
1296     }
1297   }
1298   
1299   if(linearfit) {
1300     if(ilinearfit) ilinearfit->Add(linearfit);
1301     else {
1302       ilinearfit = new AliTRDCalibraVdriftLinearFit(*linearfit);
1303       fListHist->Add(ilinearfit);
1304     }
1305   } 
1306
1307   if(exbaltfit) {
1308     if(iexbaltfit) iexbaltfit->Add(exbaltfit);
1309     else {
1310       iexbaltfit = new AliTRDCalibraExbAltFit(*exbaltfit);
1311       fListHist->Add(iexbaltfit);
1312     }
1313   } 
1314
1315   if(calibraVector) {
1316     if(icalibraVector) icalibraVector->Add(calibraVector);
1317     else {
1318       icalibraVector = new AliTRDCalibraVector(*calibraVector);
1319       fListHist->Add(icalibraVector);
1320     }
1321   }
1322   
1323 }
1324 //________________________________________________________________________________
1325 Long64_t AliTRDCalibTask::Merge(TCollection *li) {
1326   
1327   //
1328   // merge component
1329   //
1330   
1331   TIterator* iter = li->MakeIterator();
1332   AliTRDCalibTask* cal = 0;
1333
1334   while ((cal = (AliTRDCalibTask*)iter->Next())) {
1335     if (!cal->InheritsFrom(AliTRDCalibTask::Class())) {
1336       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1337       return -1;
1338     }
1339
1340     // add histograms here...
1341     this->AddTask(cal);
1342     
1343   }
1344   
1345   return 0;
1346   
1347 }
1348 //_____________________________________________________
1349 Bool_t AliTRDCalibTask::SetVersionSubversion(){
1350   //
1351   // Load Chamber Gain factors into the Tender supply
1352   //
1353   
1354   //printf("SetVersionSubversion\n");
1355
1356   //find previous entry from the UserInfo
1357   TTree *tree=((TChain*)GetInputData(0))->GetTree();
1358   if (!tree) {
1359     AliError("Tree not found in ESDhandler");
1360     return kFALSE;
1361   }
1362          
1363   TList *userInfo=(TList*)tree->GetUserInfo();
1364   if (!userInfo) {
1365     AliError("No UserInfo found in tree");
1366     return kFALSE;
1367   }
1368
1369   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
1370   if (!cdbList) {
1371     AliError("No cdbList found in UserInfo");
1372     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
1373     return kFALSE;
1374   }
1375         
1376   TIter nextCDB(cdbList);
1377   TObjString *os=0x0;
1378   while ( (os=(TObjString*)nextCDB()) ){
1379     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
1380       // Get Old gain calibration
1381       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1382       fFirstRunGain = id->GetFirstRun();
1383       fVersionGainUsed = id->GetVersion();
1384       fSubVersionGainUsed = id->GetSubVersion();
1385     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
1386       // Get Old drift velocity calibration
1387       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1388       fFirstRunVdrift = id->GetFirstRun();
1389       fVersionVdriftUsed = id->GetVersion();
1390       fSubVersionVdriftUsed = id->GetSubVersion();
1391     } else if(os->GetString().Contains("TRD/Calib/LocalGainFactor")){
1392       // Get Old drift velocity calibration
1393       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1394       fFirstRunGainLocal = id->GetFirstRun();
1395       fVersionGainLocalUsed = id->GetVersion();
1396       fSubVersionGainLocalUsed = id->GetSubVersion();
1397     } else if((os->GetString().Contains("TRD/Calib/ChamberExB")) && (!os->GetString().Contains("TRD/Calib/ChamberExBAlt"))){
1398       // Get Old drift velocity calibration
1399       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1400       fFirstRunExB = id->GetFirstRun();
1401       fVersionExBUsed = id->GetVersion();
1402       fSubVersionExBUsed = id->GetSubVersion();
1403       //printf("Version %d and subversion %d\n",fVersionExBUsed,fSubVersionExBUsed);
1404     }
1405   }
1406
1407   //printf("VersionGain %d, SubversionGain %d, VersionLocalGain %d, Subversionlocalgain %d, Versionvdrift %d, Subversionvdrift %d\n",fVersionGainUsed,fSubVersionGainUsed,fVersionGainLocalUsed,fSubVersionGainLocalUsed,fVersionVdriftUsed,fSubVersionVdriftUsed);
1408
1409   // Check
1410   if((fFirstRunGain < 0)            || 
1411      (fFirstRunGainLocal < 0)       || 
1412      (fFirstRunVdrift < 0)          || 
1413      (fVersionGainUsed < 0)         || 
1414      (fVersionGainLocalUsed < 0)    || 
1415      (fSubVersionGainUsed < 0)      || 
1416      (fSubVersionGainLocalUsed < 0) || 
1417      (fVersionVdriftUsed < 0)       || 
1418      (fSubVersionVdriftUsed < 0)) {
1419     AliError("No recent calibration found");
1420     return kFALSE;
1421   }
1422   else return kTRUE;
1423
1424 }
1425 //_________________________________________________________________________________________________________________________
1426 Bool_t AliTRDCalibTask::ParticleGood(int i) const {
1427
1428   //
1429   // Definition of good tracks
1430   //
1431
1432   
1433   AliESDtrack *track = fESD->GetTrack(i);
1434   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
1435   if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
1436   if (fabs(track->Eta())>0.8) return 0;                         // fiducial pseudorapidity
1437   Float_t r,z;
1438   track->GetImpactParametersTPC(r,z);
1439   if (fabs(z)>2.0) return 0;                          // impact parameter in z
1440   if (fabs(r)>2.0) return 0;                          // impact parameter in xy
1441   if (r==0) return 0;
1442   return 1;   
1443
1444
1445 }
1446 //______________________________________________________________________________________________________________________
1447 Bool_t AliTRDCalibTask::FindP1TrackPHtrackletV1Test(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1448 {
1449   //
1450   // Drift velocity calibration:
1451   // Fit the clusters with a straight line
1452   // From the slope find the drift velocity
1453   //
1454
1455   ////////////////////////////////////////////////
1456   //Number of points: if less than 3 return kFALSE
1457   /////////////////////////////////////////////////
1458   if(nbclusters <= 2) return kFALSE;
1459
1460   ////////////
1461   //Variables
1462   ////////////
1463   // results of the linear fit
1464   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
1465   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
1466   Double_t pointError                 = 0.0;                                // error after straight line fit 
1467   // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
1468   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
1469   Int_t    rowp                       = -1;                                 // if it crosses a pad row
1470   Float_t  tnt                        = tracklet->GetTilt();                // tan tiltingangle
1471   TLinearFitter linearFitterTracklet(2,"pol1");
1472   linearFitterTracklet.StoreData(kTRUE);  
1473  
1474   
1475   ///////////////////////////////////////////
1476   // Take the parameters of the track
1477   //////////////////////////////////////////
1478   // take now the snp, tnp and tgl from the track
1479   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1480   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber 
1481   if( TMath::Abs(snp) <  1.){
1482     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1483   } 
1484   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1485   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1486   // at the entrance
1487   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1488   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1489   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1490   // at the end with correction due to linear fit
1491   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1492   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction 
1493
1494
1495   ////////////////////////////
1496   // loop over the clusters
1497   ////////////////////////////
1498   Int_t  nbli = 0;
1499   AliTRDcluster *cl                   = 0x0;
1500   //////////////////////////////
1501   // Check no shared clusters
1502   //////////////////////////////
1503   for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
1504     cl = tracklet->GetClusters(icc);
1505     if(cl)  crossrow = 1;
1506   }
1507   //////////////////////////////////
1508   // Loop clusters
1509   //////////////////////////////////
1510   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1511     if(!(cl = tracklet->GetClusters(ic))) continue;
1512     //if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
1513     
1514     Double_t ycluster                 = cl->GetY();
1515     Int_t time                        = cl->GetPadTime();
1516     Double_t timeis                   = time/10.0;
1517     //See if cross two pad rows
1518     Int_t    row                      = cl->GetPadRow();
1519     if(rowp==-1) rowp                 = row;
1520     if(row != rowp) crossrow          = 1;
1521
1522     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1523     nbli++;  
1524
1525     
1526   }
1527   
1528   ////////////////////////////////////
1529   // Do the straight line fit now
1530   ///////////////////////////////////
1531   if(nbli <= 2){ 
1532     linearFitterTracklet.ClearPoints();  
1533     return kFALSE; 
1534   }
1535   TVectorD pars;
1536   linearFitterTracklet.Eval();
1537   linearFitterTracklet.GetParameters(pars);
1538   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/(nbli-2));
1539   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
1540   dydt        = pars[1]; 
1541   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",linearFitterTracklet->GetChisquare(),nbli,pointError,linearFitterTracklet->GetParError(1),errorpar);
1542   linearFitterTracklet.ClearPoints();  
1543  
1544   /////////////////////////
1545   // Cuts quality
1546   ////////////////////////
1547   
1548   if(nbclusters < fLow) return kFALSE;
1549   if(nbclusters > fHigh) return kFALSE;
1550   if(pointError >= 0.3) return kFALSE;
1551   if(crossrow == 1) return kTRUE;
1552   
1553   ///////////////////////
1554   // Fill
1555   //////////////////////
1556
1557   if(fDebug > 0){
1558     //Add to the linear fitter of the detector
1559     if( TMath::Abs(snp) <  1.){
1560       Double_t x = tnp-dzdx*tnt; 
1561       //if(!fLinearVdriftTest) printf("Not there\n");
1562       Double_t nbentries = fLinearVdriftTest->GetEntries();
1563       if(nbentries < (5.0*32767)) fLinearVdriftTest->Fill(x,dydt);
1564     }
1565   }
1566   
1567   return kTRUE;
1568 }
1569