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