[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalHistoComponent.cxx

// $Id$

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
//*                  for The ALICE HLT Project.                            *
//*                                                                        *
//* Permission to use, copy, modify and distribute this software and its   *
//* documentation strictly for non-commercial purposes is hereby granted   *
//* without fee, provided that the above copyright notice appears in all   *
//* copies and that both the copyright notice and this permission notice   *
//* appear in the supporting documentation. The authors make no claims     *
//* about the suitability of this software for any purpose. It is          *
//* provided "as is" without express or implied warranty.                  *
//**************************************************************************

/// @file   AliHLTGlobalHistoComponent.cxx
/// @author Matthias Richter
/// @date   2010-09-16
/// @brief  A histogramming component for global ESD properties based
///         on the AliHLTTTreeProcessor

#include "AliHLTGlobalHistoComponent.h"
#include "AliESDEvent.h"
#include "AliESDv0.h"
#include "AliKFParticle.h"
#include "AliKFVertex.h"

#include "TTree.h"
#include "TString.h"
#include <cassert>

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTGlobalHistoComponent)

AliHLTGlobalHistoComponent::AliHLTGlobalHistoComponent()
  : AliHLTTTreeProcessor()
  , fEvent(0)
  , fNofTracks(0)
  , fNofV0s(0)
  , fNofUPCpairs(0)
  , fNofContributors(0)
  , fVertexX(-99)
  , fVertexY(-99)
  , fVertexZ(-99)
  , fVertexStatus(kFALSE)
  , fTrackVariables()
  , fTrackVariablesInt()
  , fV0Variables()
  , fUPCVariables()
  , fNEvents(0)
  , fNGammas(0)
  , fNKShorts(0)
  , fNLambdas(0)
  , fNPi0s(0)
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt

}

AliHLTGlobalHistoComponent::~AliHLTGlobalHistoComponent()
{
  // see header file for class documentation
  fTrackVariables.Reset();
  fTrackVariablesInt.Reset();
  fV0Variables.Reset();
  fUPCVariables.Reset();
}

void AliHLTGlobalHistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list){
  // see header file for class documentation
  list.push_back(kAliHLTAllDataTypes);
}

AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOutputDataType(){
  // see header file for class documentation
  return kAliHLTDataTypeHistogram|kAliHLTDataOriginOut;
}

TTree* AliHLTGlobalHistoComponent::CreateTree(int /*argc*/, const char** /*argv*/){
  // create the tree and branches
  int iResult=0;
  TTree* pTree=new TTree("ESDproperties", "HLT ESD properties");
  if (!pTree) return NULL;

  const char* trackVariableNames = {
    // Note the black at the end of each name!
    "Track_pt "
    "Track_phi "
    "Track_eta "
    "Track_p "
    "Track_theta "
    "Track_Nclusters "
    "Track_status "
    "Track_charge "
    "Track_DCAr "
    "Track_DCAz "
    "Track_dEdx "
  };
  
  const char* trackIntVariableNames = {
    // Note the black at the end of each name!
    "Track_status "
  };
  
  const char* V0VariableNames = {
    // Note the black at the end of each name!
    "V0_AP "
    "V0_pt "  
    "clust1 "
    "clust2 "
    "dev1 "
    "dev2 "
    "devPrim "
    "length "
    "sigmaLength "
    "r "
  };
  
  const char* UPCVariableNames = {
    // Note the black at the end of each name!
    "nClusters_1 "
    "nClusters_2 "
    "polarity_1 "
    "polarity_2 "
    "px_1 "
    "py_1 "
    "px_2 "
    "py_2 "
  };
  
  int maxTrackCount = 20000; // FIXME: make configurable
  int maxV0Count    = 100000;
  int maxUPCCount   = 1;
    
  if ((iResult=fTrackVariables.Init(maxTrackCount, trackVariableNames))<0) {
    HLTError("failed to initialize internal structure for track properties (float)");
  }
  if ((iResult=fTrackVariablesInt.Init(maxTrackCount, trackIntVariableNames))<0) {
    HLTError("failed to initialize internal structure for track properties (int)");
  }
  if ((iResult=fV0Variables.Init(maxV0Count, V0VariableNames))<0) {
    HLTError("failed to initialize internal structure for V0 properties (float)");
  }
  if ((iResult=fUPCVariables.Init(maxUPCCount, UPCVariableNames))<0) {
    HLTError("failed to initialize internal structure for UPC properties (float)");
  }
  
  if (iResult>=0) {
    pTree->Branch("event",        &fEvent,           "event/I");
    pTree->Branch("trackcount",   &fNofTracks,       "trackcount/I");
    pTree->Branch("vertexX",      &fVertexX,         "vertexX/F");
    pTree->Branch("vertexY",      &fVertexY,         "vertexY/F");
    pTree->Branch("vertexZ",      &fVertexZ,         "vertexZ/F");
    pTree->Branch("V0",           &fNofV0s,          "V0/I");
    pTree->Branch("UPC",          &fNofUPCpairs,     "UPC/I");
    pTree->Branch("nContributors",&fNofContributors, "nContributors/I");
    pTree->Branch("vertexStatus", &fVertexStatus,    "vertexStatus/I");

    int i=0;
    // FIXME: this is a bit ugly since type 'f' and 'i' are specified
    // explicitely. Would be better to use a function like
    // AliHLTGlobalHistoVariables::GetType but could not get this working
    for (i=0; i<fTrackVariables.Variables(); i++) {
      TString specifier=fTrackVariables.GetKey(i);
      float* pArray=fTrackVariables.GetArray(specifier);
      specifier+="[trackcount]/f";
      pTree->Branch(fTrackVariables.GetKey(i), pArray, specifier.Data());
    }
    for (i=0; i<fTrackVariablesInt.Variables(); i++) {
      TString specifier=fTrackVariablesInt.GetKey(i);
      int* pArray=fTrackVariablesInt.GetArray(specifier);
      specifier+="[trackcount]/i";
      pTree->Branch(fTrackVariablesInt.GetKey(i), pArray, specifier.Data());
    }    
    for (i=0; i<fV0Variables.Variables(); i++) {
      TString specifier=fV0Variables.GetKey(i);
      float* pArray=fV0Variables.GetArray(specifier);
      specifier+="[V0]/f";
      pTree->Branch(fV0Variables.GetKey(i), pArray, specifier.Data());
    }
    for (i=0; i<fUPCVariables.Variables(); i++) {
      TString specifier=fUPCVariables.GetKey(i);
      float* pArray=fUPCVariables.GetArray(specifier);
      specifier+="[UPC]/f";
      pTree->Branch(fUPCVariables.GetKey(i), pArray, specifier.Data());
    }
  } else {
    delete pTree;
    pTree=NULL;
  }
  
  return pTree;
}

void AliHLTGlobalHistoComponent::FillHistogramDefinitions(){
  /// default histogram definitions
}

int AliHLTGlobalHistoComponent::FillTree(TTree* pTree, const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/ ){

  /// fill the tree from the ESD
  int iResult=0;
  if (!IsDataEvent()) return 0;

  ResetVariables();

  // fetch ESD from input stream
  const TObject* obj = GetFirstInputObject(kAliHLTAllDataTypes, "AliESDEvent");
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
  esd->GetStdContent();

  // fill track variables
  fNofTracks       = esd->GetNumberOfTracks();
  fVertexX         = esd->GetPrimaryVertexTracks()->GetX();
  fVertexY         = esd->GetPrimaryVertexTracks()->GetY();
  fVertexZ         = esd->GetPrimaryVertexTracks()->GetZ();
  fNofV0s          = esd->GetNumberOfV0s();
  fNofContributors = esd->GetPrimaryVertexTracks()->GetNContributors();
  fVertexStatus    = esd->GetPrimaryVertexTracks()->GetStatus();
  fNofUPCpairs     = 1;
  
  for (int i=0; i<fNofTracks; i++) {
    AliESDtrack *esdTrack = esd->GetTrack(i);
    if (!esdTrack) continue;
    
    Float_t DCAr, DCAz = -99;
    esdTrack->GetImpactParametersTPC(DCAr, DCAz);
    
    fTrackVariables.Fill("Track_pt"        , esdTrack->Pt()                      );
    fTrackVariables.Fill("Track_phi"       , esdTrack->Phi()*TMath::RadToDeg()   );
    fTrackVariables.Fill("Track_eta"       , esdTrack->Theta()                   );
    fTrackVariables.Fill("Track_p"         , esdTrack->P()                       );
    fTrackVariables.Fill("Track_theta"     , esdTrack->Theta()*TMath::RadToDeg() );
    fTrackVariables.Fill("Track_Nclusters" , esdTrack->GetTPCNcls()              );
    fTrackVariables.Fill("Track_status"    , esdTrack->GetStatus()               );
    fTrackVariables.Fill("Track_charge"    , esdTrack->Charge()                  );
    fTrackVariables.Fill("Track_DCAr"      , DCAr                        	 );
    fTrackVariables.Fill("Track_DCAz"      , DCAz                                );   
    fTrackVariables.Fill("Track_dEdx"      , esdTrack->GetTPCsignal()            );   
    fTrackVariablesInt.Fill("Track_status" , esdTrack->GetStatus()               );      
  }
  //HLTInfo("added parameters for %d tracks", fNofTracks);
  
  if(fNofTracks==2){
     AliESDtrack *esdTrack1 = esd->GetTrack(0);
     if(!esdTrack1) return 0;
     AliESDtrack *esdTrack2 = esd->GetTrack(1); 
     if(!esdTrack2) return 0;
     
     if(esdTrack1->Charge()*esdTrack2->Charge()<0){
     
       fUPCVariables.Fill("nClusters_1", esdTrack1->GetTPCNcls() );
       fUPCVariables.Fill("nClusters_2", esdTrack2->GetTPCNcls() );
       fUPCVariables.Fill("polarity_1",  esdTrack1->Charge()     );
       fUPCVariables.Fill("polarity_2",  esdTrack2->Charge()     );
       fUPCVariables.Fill("px_1",        esdTrack1->Px()         );
       fUPCVariables.Fill("py_1",        esdTrack1->Py()         );
       fUPCVariables.Fill("px_2",        esdTrack2->Px()         );
       fUPCVariables.Fill("py_2",        esdTrack2->Py()         );
    } 
  }
  
  AliKFParticle::SetField( esd->GetMagneticField() );

  //const double kKsMass = 0.49767;
  //const double kLambdaMass = 1.11568;
  //const double kPi0Mass = 0.13498;

  //std::vector<AliKFParticle> vGammas;
  
  
  for (int i=0; i<fNofV0s; i++) {
    AliESDv0 *esdV0 = esd->GetV0(i);
    if (!esdV0) continue;
    
    AliESDtrack *t1 = esd->GetTrack( esd->GetV0(i)->GetNindex());
    AliESDtrack *t2 = esd->GetTrack( esd->GetV0(i)->GetPindex());      

    AliKFParticle kf1( *t1, 11 );
    AliKFParticle kf2( *t2, 11 );

    AliKFVertex primVtx( *esd->GetPrimaryVertexTracks() );
    double dev1 = kf1.GetDeviationFromVertex( primVtx );
    double dev2 = kf2.GetDeviationFromVertex( primVtx );
    
    AliKFParticle v0( kf1, kf2 );
    double devPrim = v0.GetDeviationFromVertex( primVtx );
    primVtx+=v0;
    v0.SetProductionVertex( primVtx );

    Double_t length, sigmaLength;
    if( v0.GetDecayLength( length, sigmaLength ) ) continue;

    double dx = v0.GetX()-primVtx.GetX();
    double dy = v0.GetY()-primVtx.GetY();
    double r = sqrt(dx*dx + dy*dy);
    

    // Armenteros-Podolanski plot

    double pt=0, ap=0;
    //{
      AliKFParticle kf01 = kf1, kf02 = kf2;
      kf01.SetProductionVertex(v0);
      kf02.SetProductionVertex(v0);
      kf01.TransportToProductionVertex();
      kf02.TransportToProductionVertex();      
      double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
      double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
      double px = px1+px2, py = py1+py2, pz = pz1+pz2;
      double p = sqrt(px*px+py*py+pz*pz);
      double l1 = (px*px1 + py*py1 + pz*pz1)/p;
      double l2 = (px*px2 + py*py2 + pz*pz2)/p;
      pt = sqrt(px1*px1+py1*py1+pz1*pz1 - l1*l1);
      ap = (l2-l1)/(l1+l2);
    //}
    
//     if( 
//        t1->GetTPCNcls()>=fAPCuts[0]
//        && t2->GetTPCNcls()>=fAPCuts[0]
//        && dev1>=fAPCuts[1]
//        && dev2>=fAPCuts[1]
//        && devPrim <= fAPCuts[2]
//        && length >= fAPCuts[3]*sigmaLength
//        && length >= fAPCuts[4]
//        && r <= fAPCuts[5]
//        )
//        //{     
//        //if( fAP ) fAP->Fill( ap, pt );
//        //} 
// 
//     // Gamma finder
// 
//     bool isGamma = 0;
//     
//     if( 
//        t1->GetTPCNcls()>=fGammaCuts[0]
//        && t2->GetTPCNcls()>=fGammaCuts[0]
//        && dev1>=fGammaCuts[1]
//        && dev2>=fGammaCuts[1]
//        && devPrim <= fGammaCuts[2]
//        && length >= fGammaCuts[3]*sigmaLength
//        && length >= fGammaCuts[4]
//        && r <= fGammaCuts[5]
//        ){
//       double mass, error;	
//       v0.GetMass(mass,error); 
//       //if( fGamma ) fGamma->Fill( mass );
// 
//       if( TMath::Abs(mass)<=fGammaCuts[6]*error || TMath::Abs(mass)<=fGammaCuts[7] ){	
//         AliKFParticle gamma = v0;
//         gamma.SetMassConstraint(0);
// //      if( fGammaXY
// //          &&  t1->GetTPCNcls()>=60
// //          && t2->GetTPCNcls()>=60
// //          ) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
//         isGamma = 1;
//         fNGammas++;
//         vGammas.push_back( gamma );
//       } 	   
//     }
//     
//     if( isGamma ) continue;
// 
// 
//     // KShort finder
//     
//     bool isKs = 0;
//     
//     if( 
//        t1->GetTPCNcls()>=fKsCuts[0]
//        && t2->GetTPCNcls()>=fKsCuts[0]
//        && dev1>=fKsCuts[1]
//        && dev2>=fKsCuts[1]
//        && devPrim <= fKsCuts[2]
//        && length >= fKsCuts[3]*sigmaLength
//        && length >= fKsCuts[4]
//        && r <= fKsCuts[5]
//        ){     
//     
//       AliKFParticle piN( *t1, 211 );  
//       AliKFParticle piP( *t2, 211 );  
// 
//       AliKFParticle Ks( piN, piP );
//       Ks.SetProductionVertex( primVtx );
// 
//       double mass, error;
//       Ks.GetMass( mass, error);
//       //if( fKShort ) fKShort->Fill( mass );  
//       if( TMath::Abs( mass - kKsMass )<=fKsCuts[6]*error || TMath::Abs( mass - kKsMass )<=fKsCuts[7] ){  
//         isKs = 1;
//         fNKShorts++;
//       }
//     }
//     
//     if( isKs ) continue;
//     
//     // Lambda finder 
//     //printf("QQQQQQQQQQQQQQQQq :%f\n",fLambdaCuts[0]);
//     if( 
//        t1->GetTPCNcls()>=fLambdaCuts[0]
//        && t2->GetTPCNcls()>=fLambdaCuts[0]
//        && dev1>=fLambdaCuts[1]
//        && dev2>=fLambdaCuts[1]
//        && devPrim <= fLambdaCuts[2]
//        && length >= fLambdaCuts[3]*sigmaLength
//        && length >= fLambdaCuts[4]
//        && r <= fLambdaCuts[5]
//        && TMath::Abs( ap )>.4
//        ){
// 
//       AliKFParticle kP, kpi;
//       if( ap<0 ){ 
//         kP = AliKFParticle( *t2, 2212 );
//         kpi = AliKFParticle( *t1, 211 );
//       } else {
//         kP = AliKFParticle( *t1, 2212 );
//         kpi = AliKFParticle( *t2, 211 );
//       }
// 
//       AliKFParticle lambda = AliKFParticle( kP, kpi );
//       lambda.SetProductionVertex( primVtx );  
//       //double mass, error;
//       lambda.GetMass( Lmass, Lerror);
//       //if( fLambda ) fLambda->Fill( mass );
//       if( TMath::Abs( Lmass - kLambdaMass )<=fLambdaCuts[6]*Lerror || TMath::Abs( Lmass - kLambdaMass )<=fLambdaCuts[7] ){
//         fNLambdas++;
//       }
//    }

        
    fV0Variables.Fill("V0_AP", ap);
    fV0Variables.Fill("V0_pt", pt); 
    fV0Variables.Fill("clust1", t1->GetTPCNcls()); 
    fV0Variables.Fill("clust2", t2->GetTPCNcls()); 
    fV0Variables.Fill("dev1", dev1); 
    fV0Variables.Fill("dev2", dev2); 
    fV0Variables.Fill("devPrim", devPrim); 
    fV0Variables.Fill("length", length); 
    fV0Variables.Fill("sigmaLength", sigmaLength); 
    fV0Variables.Fill("r", r); 
      
 } // end of loop over V0s
  
  if (iResult<0) {
    // fill an empty event
    ResetVariables();
  }
  fEvent++;

  pTree->Fill();
  return iResult;
}

int AliHLTGlobalHistoComponent::ResetVariables()
{
  /// reset all filling variables
  fNofTracks=0;
  fNofV0s=0;
  fTrackVariables.ResetCount();
  fTrackVariablesInt.ResetCount();
  fV0Variables.ResetCount();
  fUPCVariables.ResetCount();
  return 0;
}

AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOriginDataType() const
{
  // get the origin of the output data
  return kAliHLTDataTypeHistogram|kAliHLTDataOriginHLT;
}
Commit	Line	Data
a4c1f5dd	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"
7a37a034	27	#include "AliESDv0.h"
	28	#include "AliKFParticle.h"
	29	#include "AliKFVertex.h"
	30
a4c1f5dd	31	#include "TTree.h"
a4c1f5dd	32	#include "TString.h"
a4c1f5dd	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)
1b89b93a	42	, fNofV0s(0)
1c1af02a	43	, fNofUPCpairs(0)
1b89b93a	44	, fNofContributors(0)
54da0d34	45	, fVertexX(-99)
	46	, fVertexY(-99)
	47	, fVertexZ(-99)
1b89b93a	48	, fVertexStatus(kFALSE)
a4c1f5dd	49	, fTrackVariables()
e7f3baeb	50	, fTrackVariablesInt()
7a37a034	51	, fV0Variables()
1c1af02a	52	, fUPCVariables()
7a37a034	53	, fNEvents(0)
	54	, fNGammas(0)
	55	, fNKShorts(0)
	56	, fNLambdas(0)
	57	, fNPi0s(0)
a4c1f5dd	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();
e7f3baeb	71	fTrackVariablesInt.Reset();
7a37a034	72	fV0Variables.Reset();
1c1af02a	73	fUPCVariables.Reset();
a4c1f5dd	74	}
a4c1f5dd	75
7a37a034	76	void AliHLTGlobalHistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list){
a4c1f5dd	77	// see header file for class documentation
	78	list.push_back(kAliHLTAllDataTypes);
	79	}
	80
7a37a034	81	AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOutputDataType(){
af747e4b	82	// see header file for class documentation
	83	return kAliHLTDataTypeHistogram\|kAliHLTDataOriginOut;
	84	}
	85
7a37a034	86	TTree* AliHLTGlobalHistoComponent::CreateTree(int /argc/, const char** /argv/){
a4c1f5dd	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 = {
e7f3baeb	93	// Note the black at the end of each name!
a4c1f5dd	94	"Track_pt "
	95	"Track_phi "
	96	"Track_eta "
54da0d34	97	"Track_p "
	98	"Track_theta "
	99	"Track_Nclusters "
97491309	100	"Track_status "
54da0d34	101	"Track_charge "
97491309	102	"Track_DCAr "
	103	"Track_DCAz "
	104	"Track_dEdx "
a4c1f5dd	105	};
54da0d34	106
e7f3baeb	107	const char* trackIntVariableNames = {
	108	// Note the black at the end of each name!
	109	"Track_status "
	110	};
	111
7a37a034	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 "
1c1af02a	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 "
49da4a19	132	"px_1 "
	133	"py_1 "
	134	"px_2 "
	135	"py_2 "
1c1af02a	136	};
7a37a034	137
af747e4b	138	int maxTrackCount = 20000; // FIXME: make configurable
7a37a034	139	int maxV0Count = 100000;
1c1af02a	140	int maxUPCCount = 1;
af747e4b	141
a4c1f5dd	142	if ((iResult=fTrackVariables.Init(maxTrackCount, trackVariableNames))<0) {
e7f3baeb	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)");
a4c1f5dd	147	}
7a37a034	148	if ((iResult=fV0Variables.Init(maxV0Count, V0VariableNames))<0) {
	149	HLTError("failed to initialize internal structure for V0 properties (float)");
	150	}
1c1af02a	151	if ((iResult=fUPCVariables.Init(maxUPCCount, UPCVariableNames))<0) {
	152	HLTError("failed to initialize internal structure for UPC properties (float)");
	153	}
a4c1f5dd	154
a4c1f5dd	155	if (iResult>=0) {
1b89b93a	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");
7a37a034	161	pTree->Branch("V0", &fNofV0s, "V0/I");
1c1af02a	162	pTree->Branch("UPC", &fNofUPCpairs, "UPC/I");
1b89b93a	163	pTree->Branch("nContributors",&fNofContributors, "nContributors/I");
1b89b93a	164	pTree->Branch("vertexStatus", &fVertexStatus, "vertexStatus/I");
e7f3baeb	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++) {
a4c1f5dd	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	}
e7f3baeb	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());
7a37a034	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());
e7f3baeb	187	}
1c1af02a	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	}
a4c1f5dd	194	} else {
	195	delete pTree;
	196	pTree=NULL;
	197	}
7a37a034	198
a4c1f5dd	199	return pTree;
	200	}
	201
7a37a034	202	void AliHLTGlobalHistoComponent::FillHistogramDefinitions(){
a4c1f5dd	203	/// default histogram definitions
	204	}
	205
7a37a034	206	int AliHLTGlobalHistoComponent::FillTree(TTree* pTree, const AliHLTComponentEventData& /evtData/, AliHLTComponentTriggerData& /trigData/ ){
7a37a034	207
a4c1f5dd	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
1b89b93a	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();
1c1af02a	227	fNofUPCpairs = 1;
54da0d34	228
a4c1f5dd	229	for (int i=0; i<fNofTracks; i++) {
	230	AliESDtrack *esdTrack = esd->GetTrack(i);
	231	if (!esdTrack) continue;
	232
97491309	233	Float_t DCAr, DCAz = -99;
	234	esdTrack->GetImpactParametersTPC(DCAr, DCAz);
	235
54da0d34	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() );
97491309	242	fTrackVariables.Fill("Track_status" , esdTrack->GetStatus() );
54da0d34	243	fTrackVariables.Fill("Track_charge" , esdTrack->Charge() );
97491309	244	fTrackVariables.Fill("Track_DCAr" , DCAr );
	245	fTrackVariables.Fill("Track_DCAz" , DCAz );
	246	fTrackVariables.Fill("Track_dEdx" , esdTrack->GetTPCsignal() );
1c1af02a	247	fTrackVariablesInt.Fill("Track_status" , esdTrack->GetStatus() );
7a37a034	248	}
1b89b93a	249	//HLTInfo("added parameters for %d tracks", fNofTracks);
1c1af02a	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() );
49da4a19	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() );
1c1af02a	267	}
	268	}
	269
7a37a034	270	AliKFParticle::SetField( esd->GetMagneticField() );
7a37a034	271
8dabf5a9	272	//const double kKsMass = 0.49767;
8dabf5a9	273	//const double kLambdaMass = 1.11568;
7a37a034	274	//const double kPi0Mass = 0.13498;
7a37a034	275
1c1af02a	276	//std::vector<AliKFParticle> vGammas;
7a37a034	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(dxdx + dydy);
	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(pxpx+pypy+pz*pz);
	319	double l1 = (pxpx1 + pypy1 + pz*pz1)/p;
	320	double l2 = (pxpx2 + pypy2 + pz*pz2)/p;
	321	pt = sqrt(px1px1+py1py1+pz1pz1 - l1l1);
	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
a4c1f5dd	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;
4868ab3f	466	fNofV0s=0;
a4c1f5dd	467	fTrackVariables.ResetCount();
e7f3baeb	468	fTrackVariablesInt.ResetCount();
c3176951	469	fV0Variables.ResetCount();
1c1af02a	470	fUPCVariables.ResetCount();
a4c1f5dd	471	return 0;
	472	}
	473
	474	AliHLTComponentDataType AliHLTGlobalHistoComponent::GetOriginDataType() const
	475	{
	476	// get the origin of the output data
4bb46850	477	return kAliHLTDataTypeHistogram\|kAliHLTDataOriginHLT;
a4c1f5dd	478	}