using now the specification provided by the base class
[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalHistoComponent.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: Matthias Richter <Matthias.Richter@ift.uib.no>        *
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   AliHLTGlobalHistoComponent.cxx
20 /// @author Matthias Richter
21 /// @date   2010-09-16
22 /// @brief  A histogramming component for global ESD properties based
23 ///         on the AliHLTTTreeProcessor
24
25 #include "AliHLTGlobalHistoComponent.h"
26 #include "AliESDEvent.h"
27 #include "AliESDv0.h"
28 #include "AliKFParticle.h"
29 #include "AliKFVertex.h"
30
31 #include "TTree.h"
32 #include "TString.h"
33 #include <cassert>
34
35 /** ROOT macro for the implementation of ROOT specific class methods */
36 ClassImp(AliHLTGlobalHistoComponent)
37
38 AliHLTGlobalHistoComponent::AliHLTGlobalHistoComponent()
39   : AliHLTTTreeProcessor()
40   , fEvent(0)
41   , fNofTracks(0)
42   , fNofV0s(0)
43   , fNofUPCpairs(0)
44   , fNofContributors(0)
45   , fVertexX(-99)
46   , fVertexY(-99)
47   , fVertexZ(-99)
48   , fVertexStatus(kFALSE)
49   , fTrackVariables()
50   , fTrackVariablesInt()
51   , fV0Variables()
52   , fUPCVariables()
53   , fNEvents(0)
54   , fNGammas(0)
55   , fNKShorts(0)
56   , fNLambdas(0)
57   , fNPi0s(0)
58 {
59   // see header file for class documentation
60   // or
61   // refer to README to build package
62   // or
63   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
64
65 }
66
67 AliHLTGlobalHistoComponent::~AliHLTGlobalHistoComponent()
68 {
69   // see header file for class documentation
70   fTrackVariables.Reset();
71   fTrackVariablesInt.Reset();
72   fV0Variables.Reset();
73   fUPCVariables.Reset();
74 }
75
76 void AliHLTGlobalHistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list){
77   // see header file for class documentation
78   list.push_back(kAliHLTAllDataTypes);
79 }
80
81 AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOutputDataType(){
82   // see header file for class documentation
83   return kAliHLTDataTypeHistogram|kAliHLTDataOriginOut;
84 }
85
86 TTree* AliHLTGlobalHistoComponent::CreateTree(int /*argc*/, const char** /*argv*/){
87   // create the tree and branches
88   int iResult=0;
89   TTree* pTree=new TTree("ESDproperties", "HLT ESD properties");
90   if (!pTree) return NULL;
91
92   const char* trackVariableNames = {
93     // Note the black at the end of each name!
94     "Track_pt "
95     "Track_phi "
96     "Track_eta "
97     "Track_p "
98     "Track_theta "
99     "Track_Nclusters "
100     "Track_status "
101     "Track_charge "
102     "Track_DCAr "
103     "Track_DCAz "
104     "Track_dEdx "
105   };
106   
107   const char* trackIntVariableNames = {
108     // Note the black at the end of each name!
109     "Track_status "
110   };
111   
112   const char* V0VariableNames = {
113     // Note the black at the end of each name!
114     "V0_AP "
115     "V0_pt "  
116     "clust1 "
117     "clust2 "
118     "dev1 "
119     "dev2 "
120     "devPrim "
121     "length "
122     "sigmaLength "
123     "r "
124   };
125   
126   const char* UPCVariableNames = {
127     // Note the black at the end of each name!
128     "nClusters_1 "
129     "nClusters_2 "
130     "polarity_1 "
131     "polarity_2 "
132     "px_1 "
133     "py_1 "
134     "px_2 "
135     "py_2 "
136   };
137   
138   int maxTrackCount = 20000; // FIXME: make configurable
139   int maxV0Count    = 100000;
140   int maxUPCCount   = 1;
141     
142   if ((iResult=fTrackVariables.Init(maxTrackCount, trackVariableNames))<0) {
143     HLTError("failed to initialize internal structure for track properties (float)");
144   }
145   if ((iResult=fTrackVariablesInt.Init(maxTrackCount, trackIntVariableNames))<0) {
146     HLTError("failed to initialize internal structure for track properties (int)");
147   }
148   if ((iResult=fV0Variables.Init(maxV0Count, V0VariableNames))<0) {
149     HLTError("failed to initialize internal structure for V0 properties (float)");
150   }
151   if ((iResult=fUPCVariables.Init(maxUPCCount, UPCVariableNames))<0) {
152     HLTError("failed to initialize internal structure for UPC properties (float)");
153   }
154   
155   if (iResult>=0) {
156     pTree->Branch("event",        &fEvent,           "event/I");
157     pTree->Branch("trackcount",   &fNofTracks,       "trackcount/I");
158     pTree->Branch("vertexX",      &fVertexX,         "vertexX/F");
159     pTree->Branch("vertexY",      &fVertexY,         "vertexY/F");
160     pTree->Branch("vertexZ",      &fVertexZ,         "vertexZ/F");
161     pTree->Branch("V0",           &fNofV0s,          "V0/I");
162     pTree->Branch("UPC",          &fNofUPCpairs,     "UPC/I");
163     pTree->Branch("nContributors",&fNofContributors, "nContributors/I");
164     pTree->Branch("vertexStatus", &fVertexStatus,    "vertexStatus/I");
165
166     int i=0;
167     // FIXME: this is a bit ugly since type 'f' and 'i' are specified
168     // explicitely. Would be better to use a function like
169     // AliHLTGlobalHistoVariables::GetType but could not get this working
170     for (i=0; i<fTrackVariables.Variables(); i++) {
171       TString specifier=fTrackVariables.GetKey(i);
172       float* pArray=fTrackVariables.GetArray(specifier);
173       specifier+="[trackcount]/f";
174       pTree->Branch(fTrackVariables.GetKey(i), pArray, specifier.Data());
175     }
176     for (i=0; i<fTrackVariablesInt.Variables(); i++) {
177       TString specifier=fTrackVariablesInt.GetKey(i);
178       int* pArray=fTrackVariablesInt.GetArray(specifier);
179       specifier+="[trackcount]/i";
180       pTree->Branch(fTrackVariablesInt.GetKey(i), pArray, specifier.Data());
181     }    
182     for (i=0; i<fV0Variables.Variables(); i++) {
183       TString specifier=fV0Variables.GetKey(i);
184       float* pArray=fV0Variables.GetArray(specifier);
185       specifier+="[V0]/f";
186       pTree->Branch(fV0Variables.GetKey(i), pArray, specifier.Data());
187     }
188     for (i=0; i<fUPCVariables.Variables(); i++) {
189       TString specifier=fUPCVariables.GetKey(i);
190       float* pArray=fUPCVariables.GetArray(specifier);
191       specifier+="[UPC]/f";
192       pTree->Branch(fUPCVariables.GetKey(i), pArray, specifier.Data());
193     }
194   } else {
195     delete pTree;
196     pTree=NULL;
197   }
198   
199   return pTree;
200 }
201
202 void AliHLTGlobalHistoComponent::FillHistogramDefinitions(){
203   /// default histogram definitions
204 }
205
206 int AliHLTGlobalHistoComponent::FillTree(TTree* pTree, const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/ ){
207
208   /// fill the tree from the ESD
209   int iResult=0;
210   if (!IsDataEvent()) return 0;
211
212   ResetVariables();
213
214   // fetch ESD from input stream
215   const TObject* obj = GetFirstInputObject(kAliHLTAllDataTypes, "AliESDEvent");
216   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
217   esd->GetStdContent();
218
219   // fill track variables
220   fNofTracks       = esd->GetNumberOfTracks();
221   fVertexX         = esd->GetPrimaryVertexTracks()->GetX();
222   fVertexY         = esd->GetPrimaryVertexTracks()->GetY();
223   fVertexZ         = esd->GetPrimaryVertexTracks()->GetZ();
224   fNofV0s          = esd->GetNumberOfV0s();
225   fNofContributors = esd->GetPrimaryVertexTracks()->GetNContributors();
226   fVertexStatus    = esd->GetPrimaryVertexTracks()->GetStatus();
227   fNofUPCpairs     = 1;
228   
229   for (int i=0; i<fNofTracks; i++) {
230     AliESDtrack *esdTrack = esd->GetTrack(i);
231     if (!esdTrack) continue;
232     
233     Float_t DCAr, DCAz = -99;
234     esdTrack->GetImpactParametersTPC(DCAr, DCAz);
235     
236     fTrackVariables.Fill("Track_pt"        , esdTrack->Pt()                      );
237     fTrackVariables.Fill("Track_phi"       , esdTrack->Phi()*TMath::RadToDeg()   );
238     fTrackVariables.Fill("Track_eta"       , esdTrack->Theta()                   );
239     fTrackVariables.Fill("Track_p"         , esdTrack->P()                       );
240     fTrackVariables.Fill("Track_theta"     , esdTrack->Theta()*TMath::RadToDeg() );
241     fTrackVariables.Fill("Track_Nclusters" , esdTrack->GetTPCNcls()              );
242     fTrackVariables.Fill("Track_status"    , esdTrack->GetStatus()               );
243     fTrackVariables.Fill("Track_charge"    , esdTrack->Charge()                  );
244     fTrackVariables.Fill("Track_DCAr"      , DCAr                                );
245     fTrackVariables.Fill("Track_DCAz"      , DCAz                                );   
246     fTrackVariables.Fill("Track_dEdx"      , esdTrack->GetTPCsignal()            );   
247     fTrackVariablesInt.Fill("Track_status" , esdTrack->GetStatus()               );      
248   }
249   //HLTInfo("added parameters for %d tracks", fNofTracks);
250   
251   if(fNofTracks==2){
252      AliESDtrack *esdTrack1 = esd->GetTrack(0);
253      if(!esdTrack1) return 0;
254      AliESDtrack *esdTrack2 = esd->GetTrack(1); 
255      if(!esdTrack2) return 0;
256      
257      if(esdTrack1->Charge()*esdTrack2->Charge()<0){
258      
259        fUPCVariables.Fill("nClusters_1", esdTrack1->GetTPCNcls() );
260        fUPCVariables.Fill("nClusters_2", esdTrack2->GetTPCNcls() );
261        fUPCVariables.Fill("polarity_1",  esdTrack1->Charge()     );
262        fUPCVariables.Fill("polarity_2",  esdTrack2->Charge()     );
263        fUPCVariables.Fill("px_1",        esdTrack1->Px()         );
264        fUPCVariables.Fill("py_1",        esdTrack1->Py()         );
265        fUPCVariables.Fill("px_2",        esdTrack2->Px()         );
266        fUPCVariables.Fill("py_2",        esdTrack2->Py()         );
267     } 
268   }
269   
270   AliKFParticle::SetField( esd->GetMagneticField() );
271
272   //const double kKsMass = 0.49767;
273   //const double kLambdaMass = 1.11568;
274   //const double kPi0Mass = 0.13498;
275
276   //std::vector<AliKFParticle> vGammas;
277   
278   
279   for (int i=0; i<fNofV0s; i++) {
280     AliESDv0 *esdV0 = esd->GetV0(i);
281     if (!esdV0) continue;
282     
283     AliESDtrack *t1 = esd->GetTrack( esd->GetV0(i)->GetNindex());
284     AliESDtrack *t2 = esd->GetTrack( esd->GetV0(i)->GetPindex());      
285
286     AliKFParticle kf1( *t1, 11 );
287     AliKFParticle kf2( *t2, 11 );
288
289     AliKFVertex primVtx( *esd->GetPrimaryVertexTracks() );
290     double dev1 = kf1.GetDeviationFromVertex( primVtx );
291     double dev2 = kf2.GetDeviationFromVertex( primVtx );
292     
293     AliKFParticle v0( kf1, kf2 );
294     double devPrim = v0.GetDeviationFromVertex( primVtx );
295     primVtx+=v0;
296     v0.SetProductionVertex( primVtx );
297
298     Double_t length, sigmaLength;
299     if( v0.GetDecayLength( length, sigmaLength ) ) continue;
300
301     double dx = v0.GetX()-primVtx.GetX();
302     double dy = v0.GetY()-primVtx.GetY();
303     double r = sqrt(dx*dx + dy*dy);
304     
305
306     // Armenteros-Podolanski plot
307
308     double pt=0, ap=0;
309     //{
310       AliKFParticle kf01 = kf1, kf02 = kf2;
311       kf01.SetProductionVertex(v0);
312       kf02.SetProductionVertex(v0);
313       kf01.TransportToProductionVertex();
314       kf02.TransportToProductionVertex();      
315       double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
316       double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
317       double px = px1+px2, py = py1+py2, pz = pz1+pz2;
318       double p = sqrt(px*px+py*py+pz*pz);
319       double l1 = (px*px1 + py*py1 + pz*pz1)/p;
320       double l2 = (px*px2 + py*py2 + pz*pz2)/p;
321       pt = sqrt(px1*px1+py1*py1+pz1*pz1 - l1*l1);
322       ap = (l2-l1)/(l1+l2);
323     //}
324     
325 //     if( 
326 //        t1->GetTPCNcls()>=fAPCuts[0]
327 //        && t2->GetTPCNcls()>=fAPCuts[0]
328 //        && dev1>=fAPCuts[1]
329 //        && dev2>=fAPCuts[1]
330 //        && devPrim <= fAPCuts[2]
331 //        && length >= fAPCuts[3]*sigmaLength
332 //        && length >= fAPCuts[4]
333 //        && r <= fAPCuts[5]
334 //        )
335 //        //{     
336 //        //if( fAP ) fAP->Fill( ap, pt );
337 //        //} 
338 // 
339 //     // Gamma finder
340 // 
341 //     bool isGamma = 0;
342 //     
343 //     if( 
344 //        t1->GetTPCNcls()>=fGammaCuts[0]
345 //        && t2->GetTPCNcls()>=fGammaCuts[0]
346 //        && dev1>=fGammaCuts[1]
347 //        && dev2>=fGammaCuts[1]
348 //        && devPrim <= fGammaCuts[2]
349 //        && length >= fGammaCuts[3]*sigmaLength
350 //        && length >= fGammaCuts[4]
351 //        && r <= fGammaCuts[5]
352 //        ){
353 //       double mass, error;    
354 //       v0.GetMass(mass,error); 
355 //       //if( fGamma ) fGamma->Fill( mass );
356 // 
357 //       if( TMath::Abs(mass)<=fGammaCuts[6]*error || TMath::Abs(mass)<=fGammaCuts[7] ){        
358 //         AliKFParticle gamma = v0;
359 //         gamma.SetMassConstraint(0);
360 // //      if( fGammaXY
361 // //          &&  t1->GetTPCNcls()>=60
362 // //          && t2->GetTPCNcls()>=60
363 // //          ) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
364 //         isGamma = 1;
365 //         fNGammas++;
366 //         vGammas.push_back( gamma );
367 //       }         
368 //     }
369 //     
370 //     if( isGamma ) continue;
371 // 
372 // 
373 //     // KShort finder
374 //     
375 //     bool isKs = 0;
376 //     
377 //     if( 
378 //        t1->GetTPCNcls()>=fKsCuts[0]
379 //        && t2->GetTPCNcls()>=fKsCuts[0]
380 //        && dev1>=fKsCuts[1]
381 //        && dev2>=fKsCuts[1]
382 //        && devPrim <= fKsCuts[2]
383 //        && length >= fKsCuts[3]*sigmaLength
384 //        && length >= fKsCuts[4]
385 //        && r <= fKsCuts[5]
386 //        ){     
387 //     
388 //       AliKFParticle piN( *t1, 211 );  
389 //       AliKFParticle piP( *t2, 211 );  
390 // 
391 //       AliKFParticle Ks( piN, piP );
392 //       Ks.SetProductionVertex( primVtx );
393 // 
394 //       double mass, error;
395 //       Ks.GetMass( mass, error);
396 //       //if( fKShort ) fKShort->Fill( mass );  
397 //       if( TMath::Abs( mass - kKsMass )<=fKsCuts[6]*error || TMath::Abs( mass - kKsMass )<=fKsCuts[7] ){  
398 //         isKs = 1;
399 //         fNKShorts++;
400 //       }
401 //     }
402 //     
403 //     if( isKs ) continue;
404 //     
405 //     // Lambda finder 
406 //     //printf("QQQQQQQQQQQQQQQQq :%f\n",fLambdaCuts[0]);
407 //     if( 
408 //        t1->GetTPCNcls()>=fLambdaCuts[0]
409 //        && t2->GetTPCNcls()>=fLambdaCuts[0]
410 //        && dev1>=fLambdaCuts[1]
411 //        && dev2>=fLambdaCuts[1]
412 //        && devPrim <= fLambdaCuts[2]
413 //        && length >= fLambdaCuts[3]*sigmaLength
414 //        && length >= fLambdaCuts[4]
415 //        && r <= fLambdaCuts[5]
416 //        && TMath::Abs( ap )>.4
417 //        ){
418 // 
419 //       AliKFParticle kP, kpi;
420 //       if( ap<0 ){ 
421 //         kP = AliKFParticle( *t2, 2212 );
422 //         kpi = AliKFParticle( *t1, 211 );
423 //       } else {
424 //         kP = AliKFParticle( *t1, 2212 );
425 //         kpi = AliKFParticle( *t2, 211 );
426 //       }
427 // 
428 //       AliKFParticle lambda = AliKFParticle( kP, kpi );
429 //       lambda.SetProductionVertex( primVtx );  
430 //       //double mass, error;
431 //       lambda.GetMass( Lmass, Lerror);
432 //       //if( fLambda ) fLambda->Fill( mass );
433 //       if( TMath::Abs( Lmass - kLambdaMass )<=fLambdaCuts[6]*Lerror || TMath::Abs( Lmass - kLambdaMass )<=fLambdaCuts[7] ){
434 //         fNLambdas++;
435 //       }
436 //    }
437
438         
439     fV0Variables.Fill("V0_AP", ap);
440     fV0Variables.Fill("V0_pt", pt); 
441     fV0Variables.Fill("clust1", t1->GetTPCNcls()); 
442     fV0Variables.Fill("clust2", t2->GetTPCNcls()); 
443     fV0Variables.Fill("dev1", dev1); 
444     fV0Variables.Fill("dev2", dev2); 
445     fV0Variables.Fill("devPrim", devPrim); 
446     fV0Variables.Fill("length", length); 
447     fV0Variables.Fill("sigmaLength", sigmaLength); 
448     fV0Variables.Fill("r", r); 
449       
450  } // end of loop over V0s
451   
452   if (iResult<0) {
453     // fill an empty event
454     ResetVariables();
455   }
456   fEvent++;
457
458   pTree->Fill();
459   return iResult;
460 }
461
462 int AliHLTGlobalHistoComponent::ResetVariables()
463 {
464   /// reset all filling variables
465   fNofTracks=0;
466   fNofV0s=0;
467   fTrackVariables.ResetCount();
468   fTrackVariablesInt.ResetCount();
469   fV0Variables.ResetCount();
470   fUPCVariables.ResetCount();
471   return 0;
472 }
473
474 AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOriginDataType() const
475 {
476   // get the origin of the output data
477   return kAliHLTDataTypeHistogram|kAliHLTDataOriginHLT;
478 }