New reader for the pedestal run and vdrift (Julian) and some bug fixing (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibTask.cxx
1 \r
2 /**************************************************************************\r
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
4  *                                                                        *\r
5  * Author: The ALICE Off-line Project.                                    *\r
6  * Contributors are mentioned in the code where appropriate.              *\r
7  *                                                                        *\r
8  * Permission to use, copy, modify and distribute this software and its   *\r
9  * documentation strictly for non-commercial purposes is hereby granted   *\r
10  * without fee, provided that the above copyright notice appears in all   *\r
11  * copies and that both the copyright notice and this permission notice   *\r
12  * appear in the supporting documentation. The authors make no claims     *\r
13  * about the suitability of this software for any purpose. It is          *\r
14  * provided "as is" without express or implied warranty.                  *\r
15  **************************************************************************/\r
16 \r
17 /////////////////////////////////////////////////////////////////////////////////\r
18 //                                                                             \r
19 // AliTRDCalibTask                                                               \r
20 //                                                                             \r
21 // Offline TRD calibration task                                \r
22 //                        \r
23 // Author:\r
24 //   R. Bailhache (R.Bailhache@gsi.de)\r
25 //                            \r
26 //////////////////////////////////////////////////////////////////////////////////////\r
27 \r
28 \r
29 \r
30 #include "Riostream.h"\r
31 #include "TChain.h"\r
32 #include "TTree.h"\r
33 #include "TFile.h"\r
34 #include "TProfile2D.h"\r
35 #include "TH2I.h"\r
36 #include "TH1F.h"\r
37 #include "TH2F.h"\r
38 #include "TList.h"\r
39 #include "TMath.h"\r
40 #include "TObject.h"\r
41 #include "TObjArray.h"\r
42 #include "TString.h"\r
43 #include "TCanvas.h"\r
44 #include "TLegend.h"\r
45 #include "TStyle.h"\r
46 #include "TLine.h"\r
47 \r
48 #include "AliAnalysisTask.h"\r
49 #include "AliAnalysisManager.h"\r
50 \r
51 #include "AliESDVertex.h"\r
52 #include "AliESDEvent.h"\r
53 #include "AliESDfriend.h"\r
54 #include "AliESDInputHandler.h"\r
55 #include "AliESDtrack.h"\r
56 #include "AliESDfriendTrack.h"\r
57 #include "AliTRDtrackV1.h"\r
58 #include "AliTRDseedV1.h"\r
59 #include "AliTRDcluster.h"\r
60 #include "AliTRDgeometry.h"\r
61 #include "AliESDtrackCuts.h"\r
62 #include "AliESDVertex.h"\r
63 #include "AliTRDCalDet.h"\r
64 \r
65 #include "AliTRDCalibraFillHisto.h"\r
66 #include "AliTRDCalibraVdriftLinearFit.h" \r
67 \r
68 #include "AliTRDcalibDB.h"\r
69 \r
70 \r
71 #include "AliTRDCalibTask.h"\r
72 \r
73 \r
74 ClassImp(AliTRDCalibTask)\r
75 \r
76 //________________________________________________________________________\r
77   AliTRDCalibTask::AliTRDCalibTask(const char *name) \r
78     : AliAnalysisTask(name, ""), fESD(0), fESDfriend(0),\r
79       fkEsdTrack(0),\r
80       fFriendTrack(0),\r
81       fCalibObject(0),\r
82       fTrdTrack(0),\r
83       fCl(0),\r
84       fListHist(0),\r
85       fTRDCalibraFillHisto(0),\r
86       fNEvents(0),\r
87       fNbTRDTrack(0),\r
88       fNbTRDTrackOffline(0),\r
89       fNbTRDTrackStandalone(0),\r
90       fNbTPCTRDtrack(0),\r
91       fNbTimeBin(0),\r
92       fNbTimeBinOffline(0),\r
93       fNbTimeBinStandalone(0),\r
94       fNbClusters(0),\r
95       fNbClustersOffline(0),\r
96       fNbClustersStandalone(0),\r
97       fNbTracklets(0),\r
98       fNbTrackletsOffline(0),\r
99       fNbTrackletsStandalone(0),\r
100       fCH2dSum(0),\r
101       fPH2dSum(0),\r
102       fCH2dSM(0),\r
103       fPH2dSM(0),\r
104       fHisto2d(kTRUE),\r
105       fVector2d(kFALSE),\r
106       fVdriftLinear(kTRUE),\r
107       fNbTimeBins(0),\r
108       fSelectedTrigger(new TObjArray()),\r
109       fRejected(kTRUE),\r
110       fEsdTrackCuts(0),\r
111       fRequirePrimaryVertex(kFALSE),\r
112       fVtxTPC(kFALSE),\r
113       fVtxSPD(kFALSE),\r
114       fMinNbContributors(0),\r
115       fLow(0),\r
116       fHigh(30),\r
117       fFillZero(kFALSE),\r
118       fNormalizeNbOfCluster(kFALSE),\r
119       fRelativeScale(0.0),\r
120       fMaxCluster(100.0),\r
121       fNbMaxCluster(2),\r
122       fOfflineTracks(kFALSE),\r
123       fStandaloneTracks(kFALSE),\r
124       fCalDetGain(0x0),\r
125       fMaxEvent(0),\r
126       fCounter(0),\r
127       fDebug(0)\r
128 {\r
129   //\r
130   // Default constructor\r
131   //\r
132 \r
133   fNz[0] = 0;\r
134   fNz[1] = 0;\r
135   fNz[2] = 0;\r
136   \r
137   fNrphi[0] = 0;\r
138   fNrphi[1] = 0;\r
139   fNrphi[2] = 0;\r
140 \r
141   // Define input and output slots here\r
142   // Input slot #0 works with a TChain\r
143   DefineInput(0, TChain::Class());\r
144     \r
145   // Output slot #0 writes into a TList container\r
146   DefineOutput(0, TList::Class());\r
147   \r
148    \r
149 }\r
150 //____________________________________________________________________________________\r
151 AliTRDCalibTask::~AliTRDCalibTask()\r
152 {\r
153   //\r
154   // AliTRDCalibTask destructor\r
155   //\r
156 \r
157   // Pointeur\r
158   if(fNEvents) delete fNEvents;\r
159   if(fNbTRDTrack) delete fNbTRDTrack;\r
160   if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;\r
161   if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;\r
162   if(fNbTPCTRDtrack) delete fNbTPCTRDtrack;\r
163   if(fNbTimeBin) delete fNbTimeBin;\r
164   if(fNbTimeBinOffline) delete fNbTimeBinOffline;\r
165   if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;\r
166   if(fNbClusters) delete fNbClusters;\r
167   if(fNbClustersOffline) delete fNbClustersOffline;\r
168   if(fNbClustersStandalone) delete fNbClustersStandalone;\r
169   if(fNbTracklets) delete fNbTracklets;\r
170   if(fNbTrackletsOffline) delete fNbTrackletsOffline;\r
171   if(fNbTrackletsStandalone) delete fNbTrackletsStandalone;\r
172   if(fCH2dSum) delete fCH2dSum;\r
173   if(fPH2dSum) delete fPH2dSum;\r
174   if(fCH2dSM) delete fCH2dSM;\r
175   if(fPH2dSM) delete fPH2dSM;\r
176   if(fCalDetGain) delete fCalDetGain;\r
177   \r
178   if(fSelectedTrigger) {\r
179     fSelectedTrigger->Delete();\r
180     delete fSelectedTrigger;\r
181   }\r
182   if(fEsdTrackCuts) {\r
183     delete fEsdTrackCuts;\r
184   }\r
185   \r
186 }\r
187 //________________________________________________________________________\r
188 void AliTRDCalibTask::ConnectInputData(Option_t *) \r
189 {\r
190   // Connect ESD or AOD here\r
191   // Called once\r
192   \r
193   TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); //pointer wird "umgecastet" auf anderen Variablentyp\r
194   if (!tree) {\r
195     //Printf("ERROR: Could not read chain from input slot 0");\r
196   } else {\r
197     \r
198     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
199     \r
200     if (!esdH) {\r
201       //Printf("ERROR: Could not get ESDInputHandler");\r
202     } else {\r
203       fESD = esdH->GetEvent();\r
204       esdH->SetActiveBranches("ESDfriend*");\r
205       //Printf("*** CONNECTED NEW EVENT ****");\r
206     }\r
207     \r
208   }\r
209 }\r
210 //________________________________________________________________________\r
211 void AliTRDCalibTask::CreateOutputObjects() \r
212 {\r
213   //\r
214   // Create the histos\r
215   //\r
216 \r
217   // Number of time bins\r
218   if(fNbTimeBins==0) {\r
219     AliTRDcalibDB *cal = AliTRDcalibDB::Instance();\r
220     fNbTimeBins = cal->GetNumberOfTimeBinsDCS();\r
221     if(fNbTimeBins <= 0){ \r
222       AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));\r
223       fNbTimeBins = 30;\r
224     }\r
225   }\r
226   \r
227   // instance calibration \r
228   fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();\r
229   fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms\r
230   fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors\r
231   fTRDCalibraFillHisto->SetCH2dOn();  // choose to calibrate the gain\r
232   fTRDCalibraFillHisto->SetPH2dOn();  // choose to calibrate the drift velocity\r
233   fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF\r
234   fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT\r
235   fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift\r
236   for(Int_t k = 0; k < 3; k++){\r
237     if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {\r
238       fTRDCalibraFillHisto->SetNz(k,fNz[k]);                                    // Mode calibration\r
239       fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]);                             // Mode calibration\r
240     }\r
241     else {\r
242       if((fNz[k] == 100) && (fNrphi[k] == 100))  {\r
243         if(fVector2d) AliInfo("The mode all together is not supported by the vector method");\r
244         fTRDCalibraFillHisto->SetAllTogether(k);\r
245       }\r
246       if((fNz[k] == 10) && (fNrphi[k] == 10))  {\r
247         if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");\r
248         fTRDCalibraFillHisto->SetPerSuperModule(k);\r
249       }\r
250     }\r
251   }\r
252   // Variables for how to fill\r
253   fTRDCalibraFillHisto->SetFillWithZero(fFillZero);\r
254   fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fNormalizeNbOfCluster); \r
255   fTRDCalibraFillHisto->SetMaxCluster(fMaxCluster);\r
256   fTRDCalibraFillHisto->SetNbMaxCluster(fNbMaxCluster);\r
257   \r
258   // Init with 30 timebins\r
259   fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos\r
260   fTRDCalibraFillHisto->SetNumberClusters(fLow); // At least 11 clusters\r
261   fTRDCalibraFillHisto->SetNumberClustersf(fHigh); // At least 11 clusters\r
262   fRelativeScale = fTRDCalibraFillHisto->GetRelativeScale(); // Get the relative scale for the gain\r
263   \r
264   // For testing only\r
265   if(fDebug > 2) fTRDCalibraFillHisto->SetDebugLevel(1); //debug stuff\r
266   \r
267   // output list\r
268   fListHist = new TList();\r
269   if(fHisto2d) {  \r
270     fListHist->Add(fTRDCalibraFillHisto->GetCH2d());\r
271     fListHist->Add(fTRDCalibraFillHisto->GetPH2d()); \r
272     fListHist->Add(fTRDCalibraFillHisto->GetPRF2d());\r
273   } \r
274   if(fVdriftLinear) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetVdriftLinearFit());\r
275   if(fVector2d) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector  \r
276   fNEvents = new TH1I("NEvents","NEvents", 2, 0, 2);\r
277   fListHist->Add(fNEvents);\r
278   \r
279   /////////////////////////////////////////\r
280   // First debug level\r
281   ///////////////////////////////////////\r
282   if(fDebug > 0) {\r
283     \r
284     fPH2dSM = new TProfile2D("PH2dSM","Nz10Nrphi10"\r
285                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)\r
286                            ,18,0,18);\r
287     fPH2dSM->SetYTitle("Det/pad groups");\r
288     fPH2dSM->SetXTitle("time [#mus]");\r
289     fPH2dSM->SetZTitle("<PH> [a.u.]");\r
290     fPH2dSM->SetStats(0);\r
291     //\r
292     fCH2dSM = new TH2I("CH2dSM","Nz10Nrphi10",50,0,300,18,0,18);\r
293     fCH2dSM->SetYTitle("Det/pad groups");\r
294     fCH2dSM->SetXTitle("charge deposit [a.u]");\r
295     fCH2dSM->SetZTitle("counts");\r
296     fCH2dSM->SetStats(0);\r
297     fCH2dSM->Sumw2();\r
298     //\r
299     fPH2dSum = new TProfile2D("PH2dSum","Nz100Nrphi100"\r
300                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)\r
301                             ,1,0,1);\r
302     fPH2dSum->SetYTitle("Det/pad groups");\r
303     fPH2dSum->SetXTitle("time [#mus]");\r
304     fPH2dSum->SetZTitle("<PH> [a.u.]");\r
305     fPH2dSum->SetStats(0);\r
306     //\r
307     fCH2dSum = new TH2I("CH2dSum","Nz100Nrphi100",50,0,300,1,0,1);\r
308     fCH2dSum->SetYTitle("Det/pad groups");\r
309     fCH2dSum->SetXTitle("charge deposit [a.u]");\r
310     fCH2dSum->SetZTitle("counts");\r
311     fCH2dSum->SetStats(0);\r
312     fCH2dSum->Sumw2();\r
313     \r
314     // Add them\r
315     fListHist->Add(fPH2dSM);\r
316     fListHist->Add(fCH2dSM);\r
317     fListHist->Add(fPH2dSum);\r
318     fListHist->Add(fCH2dSum);\r
319   }\r
320 \r
321   /////////////////////////////////////////\r
322   // Second debug level\r
323   ///////////////////////////////////////\r
324   if(fDebug > 1) {\r
325 \r
326     fNbTRDTrack = new TH1F("TRDTrack","TRDTrack",50,0,50);\r
327     fNbTRDTrack->Sumw2();\r
328     fNbTRDTrackOffline = new TH1F("TRDTrackOffline","TRDTrackOffline",50,0,50);\r
329     fNbTRDTrackOffline->Sumw2();\r
330     fNbTRDTrackStandalone = new TH1F("TRDTrackStandalone","TRDTrackStandalone",50,0,50);\r
331     fNbTRDTrackStandalone->Sumw2();\r
332     fNbTPCTRDtrack = new TH2F("NbTPCTRDtrack","NbTPCTRDtrack",100,0,100,100,0,100);\r
333     fNbTPCTRDtrack->Sumw2();\r
334     //\r
335     fNbTimeBin = new TH1F("NbTimeBin","NbTimeBin",35,0,35);\r
336     fNbTimeBin->Sumw2();\r
337     fNbTimeBinOffline = new TH1F("NbTimeBinOffline","NbTimeBinOffline",35,0,35);\r
338     fNbTimeBinOffline->Sumw2();\r
339     fNbTimeBinStandalone = new TH1F("NbTimeBinStandalone","NbTimeBinStandalone",35,0,35);\r
340     fNbTimeBinStandalone->Sumw2();\r
341     //\r
342     fNbClusters = new TH1F("NbClusters","",35,0,35);\r
343     fNbClusters->Sumw2();\r
344     fNbClustersOffline = new TH1F("NbClustersOffline","",35,0,35);\r
345     fNbClustersOffline->Sumw2();\r
346     fNbClustersStandalone = new TH1F("NbClustersStandalone","",35,0,35);\r
347     fNbClustersStandalone->Sumw2();\r
348     //\r
349     fNbTracklets = new TH1F("NbTracklets","NbTracklets",540,0.,540.);\r
350     fNbTracklets->Sumw2();\r
351     fNbTrackletsOffline = new TH1F("NbTrackletsOffline","NbTrackletsOffline",540,0.,540.);\r
352     fNbTrackletsOffline->Sumw2();\r
353     fNbTrackletsStandalone = new TH1F("NbTrackletsStandalone","NbTrackletsStandalone",540,0.,540.);\r
354     fNbTrackletsStandalone->Sumw2();\r
355    \r
356     fListHist->Add(fNbTRDTrack);\r
357     fListHist->Add(fNbTRDTrackOffline);\r
358     fListHist->Add(fNbTRDTrackStandalone);\r
359     fListHist->Add(fNbTPCTRDtrack);\r
360     \r
361     fListHist->Add(fNbTimeBin);\r
362     fListHist->Add(fNbTimeBinOffline);\r
363     fListHist->Add(fNbTimeBinStandalone);\r
364     fListHist->Add(fNbClusters);\r
365     fListHist->Add(fNbClustersOffline);\r
366     fListHist->Add(fNbClustersStandalone);\r
367     fListHist->Add(fNbTracklets);\r
368     fListHist->Add(fNbTrackletsOffline);\r
369     fListHist->Add(fNbTrackletsStandalone);\r
370     \r
371   }\r
372   \r
373  }\r
374  //________________________________________________________________________\r
375  void AliTRDCalibTask::Exec(Option_t *) \r
376  {\r
377    //\r
378    // Filling of the histos\r
379    //\r
380 \r
381    AliLog::SetGlobalLogLevel(AliLog::kError);\r
382 \r
383    if (!fESD) {\r
384      //Printf("ERROR: fESD not available");\r
385      return;\r
386    }\r
387 \r
388    //printf("Counter %d\n",fCounter);\r
389 \r
390    fCounter++;\r
391    if((fMaxEvent != 0) && (fMaxEvent < fCounter)) return;\r
392 \r
393    //printf("Counter %d\n",fCounter);\r
394 \r
395    ///////////////////////////////\r
396    // Require a primary vertex\r
397    //////////////////////////////\r
398    if(fRequirePrimaryVertex) {\r
399      const AliESDVertex* vtxESD = 0x0;\r
400      if      (fVtxTPC) vtxESD = fESD->GetPrimaryVertexTPC() ;\r
401      else if (fVtxSPD) vtxESD = fESD->GetPrimaryVertexSPD() ;\r
402      else              vtxESD = fESD->GetPrimaryVertexTracks() ;\r
403      if(!vtxESD){\r
404        return;\r
405      }\r
406      Int_t nCtrb = vtxESD->GetNContributors();\r
407      if(nCtrb < fMinNbContributors) return;\r
408    }\r
409 \r
410    //printf("Primary vertex passed\n");\r
411 \r
412    ///////////////////\r
413    // Check trigger\r
414    ///////////////////\r
415    //printf("Class Fired %s\n",(const char*)fESD->GetFiredTriggerClasses());\r
416    Bool_t pass = kTRUE;\r
417    Int_t numberOfTriggerSelected = fSelectedTrigger->GetEntriesFast();\r
418    //printf("numberofTriggerSelected %d\n",numberOfTriggerSelected);\r
419    if(fRejected) {\r
420      pass = kTRUE;\r
421      for(Int_t k = 0; k < numberOfTriggerSelected; k++){\r
422        const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);\r
423        const TString tString=obString->GetString();\r
424        if(fESD->IsTriggerClassFired((const char*)tString)) {\r
425          pass = kFALSE;\r
426        }\r
427      }\r
428    }\r
429    else {\r
430      pass = kFALSE;\r
431      for(Int_t k = 0; k < numberOfTriggerSelected; k++){\r
432        const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);\r
433        const TString tString=obString->GetString();\r
434        if(fESD->IsTriggerClassFired((const char*)tString)) {\r
435          pass = kTRUE;\r
436        }\r
437      }\r
438    }\r
439    if(!pass) return;   \r
440 \r
441    //printf("Trigger passed\n");\r
442   \r
443    fNEvents->Fill(1);\r
444 \r
445    // In total\r
446    Int_t nbTrdTracks = 0;\r
447    // standalone\r
448    Int_t nbTrdTracksStandalone = 0;\r
449    // offline\r
450    Int_t nbTrdTracksOffline = 0;\r
451    // TPC\r
452    Int_t nbtrackTPC = 0;\r
453 \r
454    Double_t nbTracks = fESD->GetNumberOfTracks();\r
455    //printf("Number of tracks %f\n",nbTracks);  \r
456 \r
457    if (nbTracks <= 0.0) {\r
458      \r
459      if(fDebug > 1) {\r
460        fNbTRDTrack->Fill(nbTrdTracks);\r
461        fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);\r
462        fNbTRDTrackOffline->Fill(nbTrdTracksOffline);\r
463      }\r
464      PostData(0, fListHist);\r
465      return;\r
466    }\r
467 \r
468    fESDfriend = (AliESDfriend *)fESD->FindListObject("AliESDfriend");\r
469    fESD->SetESDfriend(fESDfriend);\r
470    if(!fESDfriend) return;   \r
471    \r
472    //printf("has friends\n");\r
473 \r
474    ////////////////////////////////////\r
475    // Check the number of TPC tracks\r
476    ///////////////////////////////////\r
477    //printf("Nb of tracks %f\n",nbTracks);\r
478    for(Int_t itrk = 0; itrk < nbTracks; itrk++){\r
479      // Get ESD track\r
480      fkEsdTrack = fESD->GetTrack(itrk);\r
481      ULong_t status = fkEsdTrack->GetStatus(); \r
482      if(status&(AliESDtrack::kTPCout)) nbtrackTPC++;\r
483      // Check that the calibration object is here\r
484      if(fESDfriend && (fESDfriend->GetTrack(itrk))) {\r
485        fFriendTrack = new AliESDfriendTrack(*(fESDfriend->GetTrack(itrk)));\r
486        Int_t counteer = 0;\r
487        Int_t icalib=0;\r
488        while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){\r
489          if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;\r
490          counteer++;\r
491        }\r
492        if(counteer > 0) {\r
493          nbTrdTracks++;\r
494          if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {\r
495            nbTrdTracksStandalone++;\r
496          }\r
497          if((status&(AliESDtrack::kTRDin))) {\r
498            nbTrdTracksOffline++;\r
499          }\r
500        }\r
501        if(fFriendTrack) delete fFriendTrack;\r
502      }\r
503    }\r
504    //printf("Number of TPC tracks %d, TRD %d\n",nbtrackTPC,nbTrdTracks);\r
505    if((nbtrackTPC>0) && (nbTrdTracks > (3.0*nbtrackTPC))) pass = kFALSE;\r
506    \r
507    if(fDebug > 1) {\r
508      \r
509      fNbTRDTrack->Fill(nbTrdTracks);\r
510      fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);\r
511      fNbTRDTrackOffline->Fill(nbTrdTracksOffline);\r
512      fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);\r
513    \r
514    }\r
515 \r
516    if(!pass) {\r
517      PostData(0, fListHist);\r
518      return;\r
519    }\r
520    //printf("Pass\n");  \r
521 \r
522    /////////////////////////////////////\r
523    // Loop on AliESDtrack\r
524    ////////////////////////////////////\r
525    //printf("Nb of tracks %f\n",nbTracks);      \r
526    for(int itrk=0; itrk < nbTracks; itrk++){\r
527 \r
528      // Get ESD track\r
529      fkEsdTrack = fESD->GetTrack(itrk);\r
530      if(!fkEsdTrack) continue;\r
531 \r
532      // Quality cuts on the AliESDtrack\r
533      if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {\r
534        //printf("Not a good track\n");\r
535        continue;\r
536      }\r
537      \r
538      // Other cuts\r
539      Bool_t good = kTRUE;\r
540      Bool_t standalonetrack = kFALSE;\r
541      Bool_t offlinetrack = kFALSE;\r
542      ULong_t status = fkEsdTrack->GetStatus();\r
543      \r
544      if(!(fESDfriend->GetTrack(itrk)))  continue;   \r
545      \r
546      fFriendTrack = new AliESDfriendTrack(*(fESDfriend->GetTrack(itrk)));\r
547      \r
548      //////////////////////////////////////\r
549      // Loop on calibration objects\r
550      //////////////////////////////////////\r
551      Int_t icalib=0;\r
552      while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){\r
553        if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;\r
554        //printf("Find the calibration object\n");\r
555 \r
556        if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {\r
557          standalonetrack = kTRUE;\r
558        }\r
559        if((status&(AliESDtrack::kTRDin))) {\r
560          offlinetrack = kTRUE;\r
561        }\r
562        if(fOfflineTracks){\r
563          if(!offlinetrack){\r
564            good = kFALSE;\r
565          }\r
566        }\r
567        else if(fStandaloneTracks){\r
568          if(!standalonetrack){\r
569            good = kFALSE;\r
570          }\r
571        }\r
572        \r
573        fTrdTrack = (AliTRDtrackV1 *)fCalibObject;\r
574        if(good) {\r
575          fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack);\r
576          //printf("Fill fTRDCalibraFillHisto\n");\r
577        }\r
578        \r
579        //////////////////////////////////\r
580        // Debug \r
581        ////////////////////////////////\r
582 \r
583        if(fDebug > 0) {\r
584          \r
585          //printf("Enter debug\n");\r
586          \r
587          Int_t nbtracklets = 0;\r
588          \r
589          // Check some stuff\r
590          Bool_t standalonetracklet = kFALSE;  \r
591          const AliTRDseedV1 *tracklet = 0x0;\r
592          //////////////////////////////////////\r
593          // Loop tracklets\r
594          ///////////////////////////////////// \r
595          for(Int_t itr = 0; itr < 6; itr++){\r
596            \r
597            if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;\r
598            if(!tracklet->IsOK()) continue;\r
599            nbtracklets++;\r
600            standalonetracklet = kFALSE; \r
601            if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;\r
602 \r
603            Int_t nbclusters = 0;\r
604            Double_t phtb[AliTRDseedV1::kNtb];\r
605            memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));\r
606            Double_t sum = 0.0;\r
607            Float_t normalisation = 6.67;\r
608            Int_t detector = 0;\r
609            Int_t sector = 0;\r
610            //Int_t crossrow = 0;\r
611            \r
612            // Check no shared clusters\r
613            //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){\r
614            //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;\r
615            // }\r
616            \r
617            // Loop on clusters\r
618            for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){\r
619              \r
620              if(!(fCl = tracklet->GetClusters(ic))) continue;\r
621              nbclusters++;\r
622              Int_t time = fCl->GetPadTime();\r
623              Float_t ch =  tracklet->GetdQdl(ic);\r
624              Float_t qcl = TMath::Abs(fCl->GetQ());\r
625              detector = fCl->GetDetector();       \r
626              // Add the charge if shared cluster\r
627              if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {\r
628                if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {\r
629                  qcl += TMath::Abs(fCl->GetQ());\r
630                }\r
631              }\r
632              if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;\r
633              if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation;      \r
634              else sum += ch/normalisation;\r
635              \r
636              if(fDebug > 1) {\r
637                fNbTimeBin->Fill(time);\r
638                if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);\r
639                else fNbTimeBinOffline->Fill(time);\r
640              }\r
641            }\r
642            sector = AliTRDgeometry::GetSector(detector);\r
643            \r
644            if(fDebug > 1) {\r
645              fNbTracklets->Fill(detector);\r
646              if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);\r
647              else fNbTrackletsOffline->Fill(detector);\r
648            \r
649              fNbClusters->Fill(nbclusters);\r
650              if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);\r
651              else  fNbClustersOffline->Fill(nbclusters);\r
652            }       \r
653 \r
654            if(fDebug > 0) {\r
655              if((nbclusters > fLow) && (nbclusters < fHigh)){\r
656                if(fRelativeScale > 0.0) sum = sum/fRelativeScale;              \r
657                fCH2dSM->Fill(sum,sector+0.5);\r
658                fCH2dSum->Fill(sum,0.5);\r
659                Bool_t checknoise = kTRUE;\r
660                if(fMaxCluster > 0) {\r
661                  if(phtb[0] > fMaxCluster) checknoise = kFALSE;\r
662                  if(fNbTimeBins > fNbMaxCluster) {\r
663                    for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){\r
664                      if(phtb[k] > fMaxCluster) checknoise = kFALSE;\r
665                    }\r
666                  }\r
667                }\r
668                if(checknoise) {        \r
669                  for(int ic=0; ic<fNbTimeBins; ic++){\r
670                    if(fFillZero) {\r
671                      fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);\r
672                      fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);\r
673                    }\r
674                    else {\r
675                      if(phtb[ic] > 0.0) {\r
676                        fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);\r
677                        fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);\r
678                      }\r
679                    }\r
680                  }\r
681                }\r
682              }\r
683            }\r
684          } // loop on tracklets\r
685     \r
686        } // debug\r
687        \r
688      }// while calibration objects\r
689      \r
690      delete fFriendTrack;\r
691      \r
692    } // loop ESD track\r
693    \r
694    // Post output data\r
695    PostData(0, fListHist);\r
696  }     \r
697 //________________________________________________________________________\r
698 void AliTRDCalibTask::Terminate(Option_t *) \r
699 {\r
700   //\r
701   // Terminate\r
702   //\r
703   \r
704   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();\r
705 \r
706  \r
707 }\r
708 //_______________________________________________________\r
709 Bool_t AliTRDCalibTask::Load(const Char_t *filename)\r
710 {\r
711   //\r
712   // Generic container loader\r
713   //\r
714 \r
715   if(!TFile::Open(filename)){\r
716     //AliWarning(Form("Couldn't open file %s.", filename));\r
717     return kFALSE;\r
718   }\r
719   TList *o = 0x0;\r
720   if(!(o = (TList*)gFile->Get(GetName()))){\r
721     //AliWarning("Missing histogram container.");\r
722     return kFALSE;\r
723   }\r
724   fListHist = (TList*)o->Clone(GetName());\r
725   gFile->Close();\r
726   return kTRUE;\r
727 }\r
728 //________________________________________________________________________\r
729 void AliTRDCalibTask::Plot() \r
730 {\r
731   //\r
732   // Plot the histos stored in the TList\r
733   //\r
734  \r
735   if(!fListHist) return;\r
736 \r
737   /////////////////////////////////////\r
738   // Take the debug stuff\r
739   /////////////////////////////////////\r
740 \r
741   TH1I *nEvents  = (TH1I *) fListHist->FindObject("NEvents");\r
742 \r
743   TH1F *trdTrack = (TH1F *) fListHist->FindObject("TRDTrack");\r
744   TH1F *trdTrackOffline = (TH1F *) fListHist->FindObject("TRDTrackOffline");\r
745   TH1F *trdTrackStandalone = (TH1F *) fListHist->FindObject("TRDTrackStandalone");\r
746 \r
747   TH2F *tpctrdTrack = (TH2F *) fListHist->FindObject("NbTPCTRDtrack");\r
748 \r
749   TH1F *nbTimeBin = (TH1F *) fListHist->FindObject("NbTimeBin");\r
750   TH1F *nbTimeBinOffline = (TH1F *) fListHist->FindObject("NbTimeBinOffline");\r
751   TH1F *nbTimeBinStandalone = (TH1F *) fListHist->FindObject("NbTimeBinStandalone");\r
752 \r
753   TH1F *nbClusters = (TH1F *) fListHist->FindObject("NbClusters");\r
754   TH1F *nbClustersOffline = (TH1F *) fListHist->FindObject("NbClustersOffline");\r
755   TH1F *nbClustersStandalone = (TH1F *) fListHist->FindObject("NbClustersStandalone");\r
756 \r
757   TH1F *nbTracklets = (TH1F *) fListHist->FindObject("NbTracklets");\r
758   TH1F *nbTrackletsOffline = (TH1F *) fListHist->FindObject("NbTrackletsOffline");\r
759   TH1F *nbTrackletsStandalone = (TH1F *) fListHist->FindObject("NbTrackletsStandalone");\r
760   \r
761   /////////////////////////////////////\r
762   // Take the calibration objects\r
763   /////////////////////////////////////\r
764 \r
765   TH2I *ch2d = (TH2I *) fListHist->FindObject("CH2d");\r
766   TProfile2D *ph2d = (TProfile2D *) fListHist->FindObject("PH2d");\r
767 \r
768   TH2I *ch2dSum = (TH2I *) fListHist->FindObject("CH2dSum");\r
769   TProfile2D *ph2dSum = (TProfile2D *) fListHist->FindObject("PH2dSum");\r
770 \r
771   TH2I *ch2dSM = (TH2I *) fListHist->FindObject("CH2dSM");\r
772   TProfile2D *ph2dSM = (TProfile2D *) fListHist->FindObject("PH2dSM");\r
773   \r
774   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");\r
775   \r
776   ////////////////////////////////////////////////\r
777   // Add the AliTRDCalibraVdriftLinearFit\r
778   ///////////////////////////////////////////////\r
779   \r
780   Int_t first = 0;\r
781   TH2S *histolinearfitsum = 0x0;\r
782   \r
783   if(linearfit) {\r
784     for(Int_t det = 0; det < 540; det++) {\r
785       if(linearfit->GetLinearFitterHisto(det)){\r
786         if(TMath::Abs(first)<0.0001){\r
787           histolinearfitsum = linearfit->GetLinearFitterHisto(det);\r
788           first += 1;\r
789         }\r
790         else {\r
791           histolinearfitsum ->Add(linearfit->GetLinearFitterHisto(det));\r
792         }\r
793       }\r
794     }\r
795   }\r
796 \r
797   ///////////////////////////\r
798   // Style\r
799   //////////////////////////\r
800 \r
801   gStyle->SetPalette(1);\r
802   gStyle->SetOptStat(1111);\r
803   gStyle->SetOptFit(1111);\r
804   gStyle->SetPadBorderMode(0);\r
805   gStyle->SetCanvasColor(10);\r
806   gStyle->SetPadLeftMargin(0.13);\r
807   gStyle->SetPadRightMargin(0.13);\r
808 \r
809   /////////////////////////\r
810   // Plot\r
811   ////////////////////////\r
812 \r
813  if(nEvents) {\r
814    \r
815     TCanvas *debugEvents = new TCanvas("cNEvents","cNEvents",10,10,510,510);\r
816     debugEvents->cd(1);\r
817     if(nEvents) nEvents->Draw();\r
818       \r
819   }\r
820 \r
821   if(trdTrack || tpctrdTrack) {\r
822     \r
823     TCanvas *debugtrdtpcTrack = new TCanvas("TRDtracktpctrdtrack","TRDtracktpctrdtrack",10,10,510,510);\r
824     debugtrdtpcTrack->Divide(2,1);\r
825     debugtrdtpcTrack->cd(1);\r
826     if(trdTrack) trdTrack->Draw();\r
827     if(trdTrackOffline) trdTrackOffline->Draw("same");\r
828     if(trdTrackStandalone) trdTrackStandalone->Draw("same");\r
829     TLegend *leg = new TLegend(0.4,0.6,0.89,0.89);\r
830     if(trdTrack) leg->AddEntry(trdTrack,"All","p");\r
831     if(trdTrackOffline) leg->AddEntry(trdTrackOffline,"Offline","p");\r
832     if(trdTrackStandalone) leg->AddEntry(trdTrackStandalone,"Standalone","p");\r
833     leg->Draw("same");\r
834     debugtrdtpcTrack->cd(2);\r
835     if(tpctrdTrack) tpctrdTrack->Draw();\r
836     TLine *line = new TLine(0.0,0.0,100.0,100.0);\r
837     line->Draw("same");\r
838     \r
839   }\r
840  \r
841   if(nbTimeBin || nbTracklets || nbClusters) {\r
842     \r
843     TCanvas *debugTracklets = new TCanvas("TRDtimebintrackletcluster","TRDtimebintrackletcluster",10,10,510,510);\r
844     debugTracklets->Divide(3,1);\r
845     debugTracklets->cd(1);\r
846     if(nbTimeBin) nbTimeBin->Draw();\r
847     if(nbTimeBinOffline) nbTimeBinOffline->Draw("same");\r
848     if(nbTimeBinStandalone) nbTimeBinStandalone->Draw("same");\r
849     TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);\r
850     if(nbTimeBin) lega->AddEntry(nbTimeBin,"All","p");\r
851     if(nbTimeBinOffline) lega->AddEntry(nbTimeBinOffline,"Offline","p");\r
852     if(nbTimeBinStandalone) lega->AddEntry(nbTimeBinStandalone,"Standalone","p");\r
853     lega->Draw("same");\r
854     debugTracklets->cd(2);\r
855     if(nbTracklets) nbTracklets->Draw();\r
856     if(nbTrackletsOffline) nbTrackletsOffline->Draw("same");\r
857     if(nbTrackletsStandalone) nbTrackletsStandalone->Draw("same");\r
858     TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);\r
859     if(nbTracklets) legb->AddEntry(nbTracklets,"All","p");\r
860     if(nbTrackletsOffline) legb->AddEntry(nbTrackletsOffline,"Offline","p");\r
861     if(nbTrackletsStandalone) legb->AddEntry(nbTrackletsStandalone,"Standalone","p");\r
862     legb->Draw("same");\r
863     debugTracklets->cd(3);\r
864     if(nbClusters) nbClusters->Draw();\r
865     if(nbClustersOffline) nbClustersOffline->Draw("same");\r
866     if(nbClustersStandalone) nbClustersStandalone->Draw("same");\r
867     TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);\r
868     if(nbClusters) legc->AddEntry(nbClusters,"All","p");\r
869     if(nbClustersOffline) legc->AddEntry(nbClustersOffline,"Offline","p");\r
870     if(nbClustersStandalone) legc->AddEntry(nbClustersStandalone,"Standalone","p");\r
871     legc->Draw("same");\r
872   \r
873   }\r
874 \r
875   if(ch2dSum || ph2dSum || histolinearfitsum) {\r
876     \r
877     TCanvas *debugSum = new TCanvas("SumCalibrationObjects","SumCalibrationObjects",10,10,510,510);\r
878     debugSum->Divide(3,1);\r
879     debugSum->cd(1);\r
880     if(ch2dSum) ch2dSum->Draw("lego");\r
881     debugSum->cd(2);\r
882     if(ph2dSum) ph2dSum->Draw("lego");\r
883     debugSum->cd(3);\r
884     if(histolinearfitsum) histolinearfitsum->Draw();\r
885   \r
886   }\r
887 \r
888   if(ch2dSM || ph2dSM) {\r
889     \r
890     TCanvas *debugSM = new TCanvas("SMCalibrationObjects","SMCalibrationObjects",10,10,510,510);\r
891     debugSM->Divide(2,1);\r
892     debugSM->cd(1);\r
893     if(ch2dSM) ch2dSM->Draw("lego");\r
894     debugSM->cd(2);\r
895     if(ph2dSM) ph2dSM->Draw("lego");\r
896      \r
897   }\r
898 \r
899   if(ch2d || ph2d) {\r
900     \r
901     TCanvas *debug = new TCanvas("CalibrationObjects","CalibrationObjects",10,10,510,510);\r
902     debug->Divide(2,1);\r
903     debug->cd(1);\r
904     if(ch2d) ch2d->Draw("lego");\r
905     debug->cd(2);\r
906     if(ph2d) ph2d->Draw("lego");\r
907      \r
908   }\r
909  \r
910 }\r
911 \r