e75bc6fa0754828eb4b01cec7976b3c071e37f30
[u/mrichter/AliRoot.git] / HLT / global / physics / AliHLTV0HistoComponent.cxx
1 //**************************************************************************
2 //* This file is property of and copyright by the ALICE HLT Project        * 
3 //* ALICE Experiment at CERN, All rights reserved.                         *
4 //*                                                                        *
5 //* Primary Authors: Gaute Ovrebekk <ovrebekk@ift.uib.no>                  *
6 //*                  for The ALICE HLT Project.                            *
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 /** @file   AliHLTV0HistoComponent.cxx
18     @author Gaute Ovrebekk
19     @brief  Component for ploting charge in clusters
20 */
21
22 #if __GNUC__>= 3
23 using namespace std;
24 #endif
25
26 #include "AliHLTV0HistoComponent.h"
27 #include "AliCDBEntry.h"
28 #include "AliCDBManager.h"
29 #include <TFile.h>
30 #include <TString.h>
31 #include "TObjString.h"
32 #include "TObjArray.h"
33 #include "AliKFParticle.h"
34 #include "AliKFVertex.h"
35 #include "TH1F.h"
36 #include "TH2F.h"
37 #include "AliESDEvent.h"
38 #include "AliESDtrack.h"
39 #include "AliESDv0.h"
40 #include "AliHLTMessage.h"
41
42 //#include "AliHLTTPC.h"
43 //#include <stdlib.h>
44 //#include <cerrno>
45
46 /** ROOT macro for the implementation of ROOT specific class methods */
47 ClassImp(AliHLTV0HistoComponent)
48
49 AliHLTV0HistoComponent::AliHLTV0HistoComponent()
50 :
51   fGamma(0),
52   fKShort(0),
53   fLambda(0),
54   fAP(0),
55   fGammaXY(0),
56   fNGammas(0),
57   fNKShorts(0),
58   fNLambdas(0)
59 {
60   // see header file for class documentation
61   // or
62   // refer to README to build package
63   // or
64   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
65
66 }
67
68 AliHLTV0HistoComponent::~AliHLTV0HistoComponent()
69 {
70   // see header file for class documentation
71 }
72
73 // Public functions to implement AliHLTComponent's interface.
74 // These functions are required for the registration process
75
76 const char* AliHLTV0HistoComponent::GetComponentID()
77 {
78   // see header file for class documentation
79   
80   return "V0Histo";
81 }
82
83 void AliHLTV0HistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
84 {
85   // see header file for class documentation
86   list.clear();
87   list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut );
88 }
89
90 AliHLTComponentDataType AliHLTV0HistoComponent::GetOutputDataType()
91 {
92   // see header file for class documentation
93   return kAliHLTDataTypeHistogram;
94 }
95
96 void AliHLTV0HistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
97 {
98   // see header file for class documentation
99   // XXX TODO: Find more realistic values.
100   constBase = 80000;
101   inputMultiplier = 1;
102 }
103
104 AliHLTComponent* AliHLTV0HistoComponent::Spawn()
105 {
106   // see header file for class documentation
107   return new AliHLTV0HistoComponent;
108 }
109
110 int AliHLTV0HistoComponent::DoInit( int argc, const char** argv )
111 {
112   // init
113   
114   fGamma = new TH1F("hGamma","HLT:  #gamma inv mass",50,-.05,.2); 
115   fGamma->SetFillColor(kGreen);
116   fGamma->SetStats(0);
117   fKShort = new TH1F("hKShort","HLT:  K_{s}^{0} inv mass",80,0.4,.6); 
118   fKShort->SetFillColor(kGreen);
119   fKShort->SetStats(0);
120   fLambda = new TH1F("hLambda","HLT:  #Lambda^{0} inv mass",50,1.0,1.4); 
121   fLambda->SetFillColor(kGreen);
122   fLambda->SetStats(0);
123   fAP = new TH2F("hAP","HLT:  Armenteros-Podolanski plot",60,-1.,1.,60,0.,0.3);
124   fAP->SetMarkerStyle(8);
125   fAP->SetMarkerSize(0.4);
126   fAP->SetYTitle("p_{t}[GeV]");
127   fAP->SetXTitle("(p^{+}_{L}-p^{-}_{L})/(p^{+}_{L}+p^{-}_{L})");
128   fAP->SetStats(0);
129   fAP->SetOption("COLZ");
130
131   fGammaXY = new TH2F("hGammaXY","HLT:  #gamma conversions",100,-100.,100.,100,-100.,100.);
132   fGammaXY->SetMarkerStyle(8);
133   fGammaXY->SetMarkerSize(0.4);
134   fGammaXY->SetYTitle("X[cm]");
135   fGammaXY->SetXTitle("Y[cm]");
136   fGammaXY->SetStats(0);
137   fGammaXY->SetOption("COLZ");
138
139   fNGammas = 0;
140   fNKShorts = 0;
141   fNLambdas = 0;
142
143   int iResult=0;
144   TString configuration="";
145   TString argument="";
146   for (int i=0; i<argc && iResult>=0; i++) {
147     argument=argv[i];
148     if (!configuration.IsNull()) configuration+=" ";
149     configuration+=argument;
150   }
151   
152   if (!configuration.IsNull()) {
153     iResult=Configure(configuration.Data());
154   }  
155
156   return iResult; 
157 }
158   
159 int AliHLTV0HistoComponent::DoDeinit()
160 {
161   // see header file for class documentation
162   delete fGamma;
163   delete fKShort;
164   delete fLambda;
165   delete fAP;
166   delete fGammaXY;
167   return 0;
168 }
169
170 int AliHLTV0HistoComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/)
171 {
172   
173   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
174     return 0;
175
176   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
177
178     AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
179     event->GetStdContent();
180     Int_t nV0 = event->GetNumberOfV0s();
181     AliKFParticle::SetField( event->GetMagneticField() );
182
183     const double kKsMass = 0.49767;
184     const double kLambdaMass = 1.11568;
185
186     for (Int_t iv=0; iv<nV0; iv++) {
187         
188       AliESDtrack *t1=event->GetTrack( event->GetV0(iv)->GetNindex());
189       AliESDtrack *t2=event->GetTrack( event->GetV0(iv)->GetPindex());      
190
191       AliKFParticle kf1( *t1->GetInnerParam(), 11 );
192       AliKFParticle kf2( *t2->GetInnerParam(), 11 );
193
194       AliKFVertex primVtx( *event->GetPrimaryVertexTracks() );
195       double dev1 = kf1.GetDeviationFromVertex( primVtx );
196       double dev2 = kf2.GetDeviationFromVertex( primVtx );
197       
198       AliKFParticle v0( kf1, kf2 );
199       primVtx+=v0;
200       v0.SetProductionVertex( primVtx );
201
202       // Gamma finder
203       bool isGamma = 0;
204       {
205         double mass, error;       
206         v0.GetMass(mass,error);
207         if( fGamma ) fGamma->Fill( mass );
208
209         if( TMath::Abs(mass)<0.05 ){
210           AliKFParticle gamma = v0;
211           gamma.SetMassConstraint(0);
212           if( fGammaXY ) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
213           isGamma = 1;
214           fNGammas++;
215         }            
216       }
217       
218       if( isGamma ) continue;
219
220       double dx = v0.GetX()-primVtx.GetX();
221       double dy = v0.GetY()-primVtx.GetY();
222       double r = sqrt(dx*dx + dy*dy);
223
224       if( r>50. ) continue;
225
226       // AP plot
227
228       double ap=0;
229
230       {
231         AliKFParticle kf01 = kf1, kf02 = kf2;
232         kf01.SetProductionVertex(v0);
233         kf02.SetProductionVertex(v0);
234         kf01.TransportToProductionVertex();
235         kf02.TransportToProductionVertex();      
236         double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
237         double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
238         double px = px1+px2, py = py1+py2, pz = pz1+pz2;
239         double p = sqrt(px*px+py*py+pz*pz);
240         double l1 = (px*px1 + py*py1 + pz*pz1)/p;
241         double l2 = (px*px2 + py*py2 + pz*pz2)/p;
242         double pt = sqrt(px1*px1+py1*py1+pz1*pz1 - l1*l1);
243         ap = (l1-l2)/(l1+l2);
244         if( fAP ) fAP->Fill( ap, pt );
245       } 
246
247       double length = v0.GetDecayLength();
248
249       // KShort finder
250       
251       bool isKs = 0;
252       
253       if( t1->GetTPCNcls()>=60 && t2->GetTPCNcls()>=60 && length>=1.5 ){
254         
255         AliKFParticle piN( *t1->GetInnerParam(), 211 ); 
256         AliKFParticle piP( *t2->GetInnerParam(), 211 ); 
257         
258         AliKFParticle Ks( piN, piP );
259         Ks.SetProductionVertex( primVtx );
260         
261         double mass, error;
262         Ks.GetMass( mass, error);
263         if( fKShort ) fKShort->Fill( mass );
264         
265         if( TMath::Abs( mass - kKsMass )<3.5*error ){  
266           isKs = 1;
267           fNKShorts++;
268         }
269       }
270       
271       if( isKs ) continue;
272
273       // Lambda finder 
274      
275       if( t1->GetTPCNcls()>=60 && t2->GetTPCNcls()>=60 
276           && dev1>=3. && dev2>=3. 
277           && length>=4. 
278           ){
279         
280         AliKFParticle kP, kpi;
281         if( ap<0 ){ 
282           kP = AliKFParticle( *t2->GetInnerParam(), 2212 );
283           kpi = AliKFParticle( *t1->GetInnerParam(), 211 );
284         } else {
285           kP = AliKFParticle( *t1->GetInnerParam(), 2212 );
286           kpi = AliKFParticle( *t2->GetInnerParam(), 211 );
287         }
288
289         AliKFParticle lambda = AliKFParticle( kP, kpi );
290         lambda.SetProductionVertex( primVtx );  
291         double mass, error;
292         lambda.GetMass( mass, error);
293         if( fLambda ) fLambda->Fill( mass );
294
295         if( TMath::Abs( mass - kLambdaMass )<3.5*error ){  
296           fNLambdas++;
297         }       
298       }
299
300     }// V0's
301   
302     if( fGamma ) PushBack( (TObject*) fGamma, kAliHLTDataTypeHistogram,0);
303     
304     if( fKShort ) PushBack( (TObject*) fKShort, kAliHLTDataTypeHistogram,0);
305     
306     if( fLambda ) PushBack( (TObject*) fLambda, kAliHLTDataTypeHistogram, 0);
307     
308     if( fAP ) PushBack( (TObject*) fAP, kAliHLTDataTypeHistogram,0);    
309
310     if( fGammaXY ) PushBack( (TObject*) fGammaXY, kAliHLTDataTypeHistogram,0);
311   }  
312   
313   HLTInfo("Found %d Gammas, %d KShorts and %d Lambdas total", fNGammas, fNKShorts, fNLambdas );
314   
315   return 0;
316 }
317
318 int AliHLTV0HistoComponent::Configure(const char* arguments)
319 {
320   
321   int iResult=0;
322   if (!arguments) return iResult;
323   
324   TString allArgs=arguments;
325   TString argument;
326   
327   TObjArray* pTokens=allArgs.Tokenize(" ");
328   
329   if (pTokens) {
330     for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
331       argument=((TObjString*)pTokens->At(i))->GetString();
332       if (argument.IsNull()) continue;
333       
334       if (argument.CompareTo("-plot-all")==0) {
335         HLTInfo("Ploting charge of all clusters");
336         //fPlotAll = kTRUE;
337         continue;
338       }
339       
340       else if (argument.CompareTo("-plot-trackclusters")==0) {
341         HLTInfo("Ploting charge of clusters used on a track");
342         //fPlotAll = kFALSE;
343         continue;
344       }
345       else {
346         HLTError("unknown argument %s", argument.Data());
347         iResult=-EINVAL;
348         break;
349       }
350     }
351     delete pTokens;
352   }
353   
354   return iResult;
355 }
356
357 int AliHLTV0HistoComponent::Reconfigure(const char* cdbEntry, const char* chainId)
358 {
359   // see header file for class documentation
360   int iResult=0;
361   const char* path="HLT/ConfigTPC/KryptonHistoComponent";
362   const char* defaultNotify="";
363   if (cdbEntry) {
364     path=cdbEntry;
365     defaultNotify=" (default)";
366   }
367   if (path) {
368     HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
369     AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
370     if (pEntry) {
371       TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
372       if (pString) {
373         HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
374         iResult=Configure(pString->GetString().Data());
375       } else {
376         HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
377       }
378     } else {
379       HLTError("can not fetch object \"%s\" from CDB", path);
380     }
381   }
382
383   return iResult;
384 }