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