exclude unnecessary cuts; remove cuts on ITS clusters as those remove all data points...
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCdEdxMonitoringComponent.cxx
1 // $Id$
2
3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project        * 
5 //* ALICE Experiment at CERN, All rights reserved.                         *
6 //*                                                                        *
7 //* Primary Authors: Per-Ivar Lønne <ploenne@cern.ch>                      *
8 //*                  for The ALICE HLT Project.                            *
9 //*                                                                        *
10 //* Permission to use, copy, modify and distribute this software and its   *
11 //* documentation strictly for non-commercial purposes is hereby granted   *
12 //* without fee, provided that the above copyright notice appears in all   *
13 //* copies and that both the copyright notice and this permission notice   *
14 //* appear in the supporting documentation. The authors make no claims     *
15 //* about the suitability of this software for any purpose. It is          *
16 //* provided "as is" without express or implied warranty.                  *
17 //**************************************************************************/
18
19 /// @file   AliHLTTPCdEdxMonitoringComponent.cxx
20 /// @author Per-Ivar Lønne, Jochen Thaeder, Matthias Richter, Alexander Kalweit
21 /// @date   21.08.2011
22 /// @brief  Component for reading ESD from chain and produce a dEdx monitoring plot
23 ///
24
25 #if __GNUC__>= 3
26 using namespace std;
27 #endif
28
29
30 #include "TSystem.h"
31 #include "AliESDtrackCuts.h"
32 #include <AliHLTDAQ.h>
33 #include "AliESDEvent.h"
34 #include "AliESDtrack.h"
35 #include "AliLog.h"
36 #include "TMap.h"
37 #include "TObjString.h"
38 #include "AliHLTDataTypes.h"
39 #include "AliHLTComponentBenchmark.h"
40 #include "TROOT.h"
41 #include "TStyle.h"
42 #include "TColor.h"
43 #include "TH2F.h"
44
45 #include "AliHLTTPCdEdxMonitoringComponent.h"
46
47 /** ROOT macro for the implementation of ROOT specific class methods */
48 ClassImp( AliHLTTPCdEdxMonitoringComponent )
49   
50 AliHLTTPCdEdxMonitoringComponent::AliHLTTPCdEdxMonitoringComponent() 
51   :
52   AliHLTProcessor(),
53   fESDTrackCuts(NULL),
54   fHist(),
55   fxbins(),
56   fxmin(),
57   fxmax(),
58   fybins(),
59   fymin(),
60   fymax()
61 {
62   //Constructor
63 }
64
65 AliHLTTPCdEdxMonitoringComponent::~AliHLTTPCdEdxMonitoringComponent()
66 {
67   //Destructor
68 }
69
70 // ######################################################################### //
71 const char* AliHLTTPCdEdxMonitoringComponent::GetComponentID()
72 {
73   // get component id
74   return "TPCdEdxMonitoring";
75 }
76
77 // ######################################################################### //
78 void AliHLTTPCdEdxMonitoringComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list)
79 {
80   // get list of input data types
81   list.push_back(kAliHLTDataTypeESDObject|kAliHLTDataOriginAny);
82 }
83
84 // ######################################################################### //
85 AliHLTComponentDataType AliHLTTPCdEdxMonitoringComponent::GetOutputDataType()
86 {
87   // get the output data size
88   return kAliHLTDataTypeHistogram|kAliHLTDataOriginHLT;
89 }
90
91 // ######################################################################### //
92 void AliHLTTPCdEdxMonitoringComponent::GetOutputDataSize( unsigned long& constBase, Double_t& inputMultiplier )
93 {
94   // get output size estimator
95   constBase = 4096;
96   inputMultiplier = 0;
97 }
98
99 // ######################################################################### //
100 void AliHLTTPCdEdxMonitoringComponent::GetOCDBObjectDescription( TMap* const targetMap)
101 {
102   // Get a list of OCDB object description.
103   // The list of objects is provided in a TMap
104   // - key: complete OCDB path, e.g. GRP/GRP/Data
105   // - value: short description why the object is needed
106   // Key and value objects created inside this class go into ownership of
107   // target TMap.
108   
109     if (!targetMap) return;
110     targetMap->Add(new TObjString("HLT/ConfigTPC/TPCdEdxMonitoring"),
111     new TObjString("configuration object"));
112     targetMap->Add(new TObjString("GRP/GRP/Data"),
113     new TObjString("GRP object"));
114   
115 }
116
117 // ######################################################################### //
118 AliHLTComponent* AliHLTTPCdEdxMonitoringComponent::Spawn()
119 {
120   // Spawn function, return new class instance
121   return new AliHLTTPCdEdxMonitoringComponent;
122 }
123
124 // ######################################################################### //
125 Int_t AliHLTTPCdEdxMonitoringComponent::DoInit( Int_t argc, const char** argv )
126 {
127   // init the component
128   Int_t iResult=0;
129   fESDTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","HLT");
130
131   fxbins=Int_t(300);
132   fxmin=Double_t(-2);
133   fxmax=Double_t(2);
134   fybins=Int_t(500);
135   fymin=Double_t(0);
136   fymax=Double_t(500);
137   
138 if (!fESDTrackCuts)
139     {
140       iResult=-ENOMEM;
141     }
142
143   if (iResult<0)
144     {
145       if (fESDTrackCuts)
146         delete fESDTrackCuts;
147       fESDTrackCuts = NULL;
148     }      
149        
150   if (iResult>=0) 
151     {
152       SetDefaultConfiguration();
153       TString cdbPath="HLT/ConfigTPC/"; 
154       cdbPath+=GetComponentID();
155       iResult=ConfigureFromCDBTObjString(cdbPath);
156     
157   
158       if (iResult>=0) 
159         {
160           iResult=ConfigureFromArgumentString(argc, argv);      
161         }
162     }
163
164   if (iResult>=0) {
165     HLTInfo("ESD track cuts : %s",fESDTrackCuts->GetTitle() );
166   }
167   fHist = new TH2F("hHLT", "HLT", fxbins, fxmin, fxmax, fybins, fymin, fymax);
168   fHist->GetXaxis()->SetTitle("momentum/charge #frac{p}{z}  (GeV/c)");
169   fHist->GetYaxis()->SetTitle("dE/dx in TPC (a.u.)");
170   Plotstyle();
171   return iResult;
172 }
173   
174
175 // ######################################################################### //
176 Int_t AliHLTTPCdEdxMonitoringComponent::DoDeinit()
177 {
178   // component cleanup, delete all instances of helper classes here
179   if (fESDTrackCuts)
180     delete fESDTrackCuts;
181   fESDTrackCuts = NULL;
182   if (fHist)
183     delete fHist;
184   fHist=NULL;
185   return 0;
186 }
187
188
189 // ######################################################################### //
190 Int_t AliHLTTPCdEdxMonitoringComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/,
191                                               AliHLTComponentTriggerData& /*trigData*/)
192 {
193   // event processing function
194   float sig;
195   Double_t ptot;
196
197   // check if this is a data event, there are a couple of special events
198   // which should be ignored for normal processing
199   if (!IsDataEvent()) return 0;
200
201   const TObject* obj = GetFirstInputObject(kAliHLTAllDataTypes, "AliESDEvent");
202
203   // input objects are not supposed to be changed by the component, so they
204   // are defined const. However, the implementation of AliESDEvent does not
205   // support this and we need the const_cast
206   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
207   if (esd != NULL) {
208     esd->GetStdContent();
209     for (Int_t i = 0; i < esd->GetNumberOfTracks(); i++) {
210       AliESDtrack* track = esd->GetTrack(i);
211       sig=track->GetTPCsignal();
212       if(!fESDTrackCuts->AcceptTrack(track)) continue;
213       if (!track->GetInnerParam()) continue;
214       ptot = track->GetInnerParam()->GetP()*track->GetSign();
215       fHist->Fill(ptot, sig);
216     }
217   }
218   // publish the histogram
219   PushBack(fHist, kAliHLTDataTypeHistogram | kAliHLTDataOriginHLT);
220
221   return 0;
222 }
223
224 // ######################################################################### //
225 Int_t AliHLTTPCdEdxMonitoringComponent::ScanConfigurationArgument(Int_t argc, const char** argv)
226 {
227   // Scan configuration arguments
228   // Return the number of processed arguments
229   //        -EPROTO if argument format error (e.g. number expected but not found)
230   //
231   // The AliHLTComponent base class implements a parsing loop for argument strings and
232   // arrays of strings which is invoked by ConfigureFromArgumentString/ConfigureFromCDBTObjString
233   // The component needs to implement ScanConfigurationArgument in order to decode the arguments.
234   
235   if (argc<=0) return 0;
236   Int_t ii =0;
237   TString argument=argv[ii];
238   
239   if (argument.IsNull()) return 0;
240   
241   if( !fESDTrackCuts){
242     HLTError("No ESD track cuts availible");
243     return -ENOMEM;
244   }
245   
246
247   //**********************************//
248   //        Histogram Binning         //
249   //**********************************//
250
251
252   // -xbins
253   if (argument.CompareTo("-xbins")==0) 
254     {
255       if (++ii>=argc) return -EPROTO;
256       argument=argv[ii];
257
258       fxbins = argument.Atoi();
259
260       return 2;
261     }
262
263   // -xmin
264   if (argument.CompareTo("-xmin")==0) 
265     {
266       if (++ii>=argc) return -EPROTO;
267       argument=argv[ii];
268
269       fxmin = argument.Atoi();
270
271       return 2;
272     }
273
274   if (argument.CompareTo("-xmax")==0) 
275     {
276       if (++ii>=argc) return -EPROTO;
277       argument=argv[ii];
278
279       fxmax = argument.Atoi();
280
281       return 2;
282     }
283
284   // -xbins
285   if (argument.CompareTo("-ybins")==0) 
286     {
287       if (++ii>=argc) return -EPROTO;
288       argument=argv[ii];
289
290       fybins = argument.Atoi();
291
292       return 2;
293     }
294
295   // -xmin
296   if (argument.CompareTo("-ymin")==0) 
297     {
298       if (++ii>=argc) return -EPROTO;
299       argument=argv[ii];
300
301       fymin = argument.Atoi();
302
303       return 2;
304     }
305
306   if (argument.CompareTo("-ymax")==0) 
307     {
308       if (++ii>=argc) return -EPROTO;
309       argument=argv[ii];
310
311       fymax = argument.Atoi();
312
313       return 2;
314     }
315
316
317   //**********************************//
318   //           Track Cuts             //
319   //**********************************//
320
321  // -maxpt
322   if (argument.CompareTo("-maxpt")==0) {
323     if (++ii>=argc) return -EPROTO;
324     argument=argv[ii];
325
326     Float_t minPt, maxPt;
327     fESDTrackCuts->GetPtRange(minPt,maxPt);
328     maxPt = argument.Atof(); 
329     fESDTrackCuts->SetPtRange(minPt,maxPt);
330
331     TString title = fESDTrackCuts->GetTitle();
332     if (!title.CompareTo("No track cuts")) title = "";
333     else title += " && ";
334     title += Form("p_t < %f", maxPt);
335     fESDTrackCuts->SetTitle(title);
336     return 2;
337   }    
338
339   // -minpt
340   if (argument.CompareTo("-minpt")==0) {
341     if (++ii>=argc) return -EPROTO;
342     argument=argv[ii];
343
344     Float_t minPt, maxPt;
345     fESDTrackCuts->GetPtRange(minPt,maxPt);
346     minPt = argument.Atof(); 
347     fESDTrackCuts->SetPtRange(minPt,maxPt);
348
349     TString title = fESDTrackCuts->GetTitle();
350     if (!title.CompareTo("No track cuts")) title = "";
351     else title += " && ";
352     title += Form("p_t > %f", minPt);
353     fESDTrackCuts->SetTitle(title);
354     return 2;
355   }    
356
357   // -min-ldca
358   // minimum longitudinal dca to vertex
359   if (argument.CompareTo("-min-ldca")==0) {
360     if (++ii>=argc) return -EPROTO;
361     argument=argv[ii];
362
363     fESDTrackCuts->SetMinDCAToVertexZ(argument.Atof());
364     TString title = fESDTrackCuts->GetTitle();
365     if (!title.CompareTo("No track cuts")) title = "";
366     else title += " && ";
367     title += Form("DCAz > %f", argument.Atof());
368     fESDTrackCuts->SetTitle(title);
369     return 2;
370   }
371   
372   // -max-ldca
373   // maximum longitudinal dca to vertex
374   if (argument.CompareTo("-max-ldca")==0) {
375     if (++ii>=argc) return -EPROTO;
376     argument=argv[ii];
377
378     fESDTrackCuts->SetMaxDCAToVertexZ(argument.Atof());
379     TString title = fESDTrackCuts->GetTitle();
380     if (!title.CompareTo("No track cuts")) title = "";
381     else title += " && ";
382     title += Form("DCAz < %f", argument.Atof());
383     fESDTrackCuts->SetTitle(title);
384     return 2;
385   }
386
387   // -min-tdca
388   // minimum transverse dca to vertex
389   if (argument.CompareTo("-min-tdca")==0) {
390     if (++ii>=argc) return -EPROTO;
391     argument=argv[ii];
392
393     fESDTrackCuts->SetMinDCAToVertexXY(argument.Atof());
394     TString title = fESDTrackCuts->GetTitle();
395     if (!title.CompareTo("No track cuts")) title = "";
396     else title += " && ";
397     title += Form("DCAr > %f", argument.Atof());
398     fESDTrackCuts->SetTitle(title);
399     return 2;
400   }
401   
402   // -max-tdca
403   // maximum transverse dca to vertex
404   if (argument.CompareTo("-max-tdca")==0) {
405     if (++ii>=argc) return -EPROTO;
406     argument=argv[ii];
407
408     fESDTrackCuts->SetMaxDCAToVertexXY(argument.Atof());
409     TString title = fESDTrackCuts->GetTitle();
410     if (!title.CompareTo("No track cuts")) title = "";
411     else title += " && ";
412     title += Form("DCAr < %f", argument.Atof());
413     fESDTrackCuts->SetTitle(title);
414     return 2;
415   }
416
417   // -etarange
418   // +/- eta 
419   if (argument.CompareTo("-etarange")==0) {
420     if (++ii>=argc) return -EPROTO;
421     argument=argv[ii];
422     Float_t eta = argument.Atof();
423
424     fESDTrackCuts->SetEtaRange(-eta,eta);     
425     TString title = fESDTrackCuts->GetTitle();
426     if (!title.CompareTo("No track cuts")) title = "";
427     else title += " && ";
428     title += Form("Eta[%f,%f]", argument.Atof(),argument.Atof());
429     fESDTrackCuts->SetTitle(title);
430     return 2;
431   }
432
433   // -minNClsTPC
434   // minimum clusters in TPC
435   if (argument.CompareTo("-minNClsTPC")==0) {
436     if (++ii>=argc) return -EPROTO;
437     argument=argv[ii];
438
439     Int_t ncls;
440     ncls = Int_t(argument.Atof()); 
441     fESDTrackCuts->SetMinNClustersTPC(ncls);  
442
443     TString title = fESDTrackCuts->GetTitle();
444     if (!title.CompareTo("No track cuts")) title = "";
445     else title += " && ";
446     title += Form("minNClsTPC < %i", ncls);
447     fESDTrackCuts->SetTitle(title);
448     return 2;
449   }    
450
451
452
453   // unknown argument
454   return -EINVAL;
455   
456   return 0;
457   
458 }
459
460 // ######################################################################### //
461 Int_t AliHLTTPCdEdxMonitoringComponent::Reconfigure(const char* cdbEntry, const char* chainId)
462 {
463   // reconfigure the component from the specified CDB entry, or default CDB entry
464   HLTInfo("reconfigure '%s' from entry %s", chainId, cdbEntry);
465
466   Int_t iResult=0;
467   TString cdbPath;
468   if (cdbEntry)
469     {
470       cdbPath=cdbEntry;
471     } else {
472       cdbPath="HLT/ConfigTPC/";
473       cdbPath+=GetComponentID();
474     }
475
476   iResult=ConfigureFromCDBTObjString(cdbPath); //// Or use return 0, and skip this line?
477
478   return iResult;
479 }
480
481 // ######################################################################### //
482 Int_t AliHLTTPCdEdxMonitoringComponent::ReadPreprocessorValues(const char* modules)
483 {
484   // read the preprocessor values for the detectors in the modules list
485   Int_t iResult=0;
486   TString detectors(modules!=NULL?modules:"");
487   HLTInfo("read preprocessor values for detector(s): %s", detectors.IsNull()?"none":detectors.Data());
488   return iResult; //Done differently in AliHLTMultiplicityCorrelationsComponent...
489 }
490
491
492
493 // ######################################################################### //
494 void AliHLTTPCdEdxMonitoringComponent::SetDefaultConfiguration() 
495 {
496   if (fESDTrackCuts)
497     {
498       //fESDTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
499       fESDTrackCuts->SetEtaRange(-0.8,+0.8);
500       fESDTrackCuts->SetPtRange(0.15,1e10);
501       // 2011-10-28 investigation by Per Ivar and Alexander
502       // the following cuts are not needed
503       /*
504       fESDTrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2);  // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
505       fESDTrackCuts->SetMaxNsigmaToVertex(3);
506       fESDTrackCuts->SetRequireSigmaToVertex(kTRUE);
507       */
508       fESDTrackCuts->SetAcceptKinkDaughters(kFALSE);
509       fESDTrackCuts->SetMinNClustersTPC(70);
510       fESDTrackCuts->SetMaxChi2PerClusterTPC(4);
511       fESDTrackCuts->SetMaxDCAToVertexXY(3);
512       fESDTrackCuts->SetMaxDCAToVertexZ(3);
513       fESDTrackCuts->SetRequireTPCRefit(kTRUE);
514       //fESDTrackCuts->SetRequireITSRefit(kTRUE); //Kills HLT simulated reconstructions?
515       // 2011-10-28 investigation by Per Ivar and Alexander
516       // those cuts remove all data points because the filling is different in the
517       // HLT, needs further investigation
518       /*
519       fESDTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
520       fESDTrackCuts->SetMinNClustersITS(3);
521       */
522       
523     }
524  fxbins=300;
525  fxmin=-2;
526  fxmax=2;
527  fybins=500;  
528  fymin=0;
529  fymax=500;
530   return;
531 }
532
533
534
535 // ######################################################################### //
536 void AliHLTTPCdEdxMonitoringComponent::Plotstyle()
537 {
538   gROOT->SetStyle("Plain");
539   
540   gStyle->SetCanvasBorderMode(0);
541   gStyle->SetCanvasColor(10);
542   gStyle->SetCanvasDefH(550);
543   gStyle->SetCanvasDefW(575);
544   gStyle->SetPadBorderMode(0);
545   gStyle->SetPadColor(10);
546   gStyle->SetPadTickX(1);
547   gStyle->SetPadTickY(1);
548   gStyle->SetStatColor(10);
549   gStyle->SetFrameFillColor(10);
550   gStyle->SetPalette(1,0);
551   
552   //
553   
554   //gStyle->SetStatX(0.7);
555   //gStyle->SetStatW(0.2);
556   //gStyle->SetLabelOffset(1.2);
557   //gStyle->SetLabelFont(72);
558   //gStyle->SetLabelSize(0.6);
559   //gStyle->SetTitleOffset(1.2);
560   gStyle->SetTitleFontSize(0.04);
561   
562   
563   gStyle->SetOptStat(10);
564   gStyle->SetLineWidth(2);
565   gStyle->SetMarkerSize(1.0);
566   gStyle->SetTextSize(0.04);
567   gStyle->SetTitleSize(0.04,"xyz");
568   gStyle->SetLabelSize(0.04,"xyz");
569   gStyle->SetLabelOffset(0.02,"xyz");
570   gStyle->SetLabelFont(42,"xyz");
571   gStyle->SetTitleOffset(1.3,"x");
572   gStyle->SetTitleOffset(1.6,"y");
573   gStyle->SetTitleOffset(1.6,"z");
574   gStyle->SetTitleFont(42,"xyz");
575   gStyle->SetTitleColor(1,"xyz");
576   //gStyle->SetPadTopMargin(0.1);
577   //gStyle->SetPadRightMargin(0.1);
578   gStyle->SetPadBottomMargin(0.2);
579   gStyle->SetPadLeftMargin(0.2);
580   
581   
582   const Int_t NCont=255;
583   const Int_t NRGBs = 5;
584   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
585   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
586   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
587   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
588   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
589   gStyle->SetNumberContours(NCont); 
590 }