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 }
 479
 480 // AliHLTUInt32_t AliHLTGlobalHistoComponent::GetDataSpec() const;
 481 // {
 482 //   // get specifications of the output data
 483 //   return 0;
 484 // }