HLT/global/physics/AliHLTV0HistoComponent.cxx

   1 // $Id$
   2 //**************************************************************************
   3 //* This file is property of and copyright by the ALICE HLT Project        *
   4 //* ALICE Experiment at CERN, All rights reserved.                         *
   5 //*                                                                        *
   6 //* Primary Authors: S.Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>    *
   7 //*                  for The ALICE HLT Project.                            *
   8 //*                                                                        *
   9 //* Permission to use, copy, modify and distribute this software and its   *
  10 //* documentation strictly for non-commercial purposes is hereby granted   *
  11 //* without fee, provided that the above copyright notice appears in all   *
  12 //* copies and that both the copyright notice and this permission notice   *
  13 //* appear in the supporting documentation. The authors make no claims     *
  14 //* about the suitability of this software for any purpose. It is          *
  15 //* provided "as is" without express or implied warranty.                  *
  16 //**************************************************************************
  17
  18 /** @file   AliHLTV0HistoComponent.cxx
  19     @author Sergey Gorbunov
  20     @brief  Component for ploting charge in clusters
  21 */
  22
  23 #if __GNUC__>= 3
  24 using namespace std;
  25 #endif
  26
  27 #include "AliHLTV0HistoComponent.h"
  28 #include "AliCDBEntry.h"
  29 #include "AliCDBManager.h"
  30 #include <TFile.h>
  31 #include <TString.h>
  32 #include "TObjString.h"
  33 #include "TObjArray.h"
  34 #include "AliKFParticle.h"
  35 #include "AliKFVertex.h"
  36 #include "TH1F.h"
  37 #include "TH2F.h"
  38 #include "TSystem.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDv0.h"
  42 #include "AliHLTMessage.h"
  43 #include "TTimeStamp.h"
  44
  45 //#include "AliHLTTPC.h"
  46 //#include <stdlib.h>
  47 //#include <cerrno>
  48
  49 /** ROOT macro for the implementation of ROOT specific class methods */
  50 ClassImp(AliHLTV0HistoComponent)
  51
  52 AliHLTV0HistoComponent::AliHLTV0HistoComponent() :
  53   fUID(0),
  54   fGamma(0),
  55   fKShort(0),
  56   fLambda(0),
  57   fPi0(0),
  58   fAP(0),
  59   fGammaXY(0),
  60   fNEvents(0),
  61   fNGammas(0),
  62   fNKShorts(0),
  63   fNLambdas(0),
  64   fNPi0s(0)
  65 {
  66   // see header file for class documentation
  67   // or
  68   // refer to README to build package
  69   // or
  70   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  71   for( int i=0; i<8; i++){
  72     fGammaCuts[i] = 0;
  73     fKsCuts[i] = 0;
  74     fLambdaCuts[i] = 0;
  75     fAPCuts[i] = 0;
  76   }
  77 }
  78
  79 AliHLTV0HistoComponent::~AliHLTV0HistoComponent()
  80 {
  81   // see header file for class documentation
  82 }
  83
  84 // Public functions to implement AliHLTComponent's interface.
  85 // These functions are required for the registration process
  86
  87 const char* AliHLTV0HistoComponent::GetComponentID()
  88 {
  89   // see header file for class documentation
  90
  91   return "V0Histo";
  92 }
  93
  94 void AliHLTV0HistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  95 {
  96   // see header file for class documentation
  97   list.clear();
  98   list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut );
  99 }
 100
 101 AliHLTComponentDataType AliHLTV0HistoComponent::GetOutputDataType()
 102 {
 103   // see header file for class documentation
 104   return kAliHLTDataTypeHistogram  | kAliHLTDataOriginOut;
 105 }
 106
 107 void AliHLTV0HistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 108 {
 109   // see header file for class documentation
 110   // XXX TODO: Find more realistic values.
 111   constBase = 80000;
 112   inputMultiplier = 0;
 113 }
 114
 115 AliHLTComponent* AliHLTV0HistoComponent::Spawn()
 116 {
 117   // see header file for class documentation
 118   return new AliHLTV0HistoComponent;
 119 }
 120
 121 int AliHLTV0HistoComponent::DoInit( int argc, const char** argv ) {
 122   // init
 123
 124   fUID = 0;
 125
 126   fGamma = new TH1F("hGamma","HLT:  #gamma inv mass",50,-.06,.2);
 127   fGamma->SetFillColor(kGreen);
 128   fGamma->SetStats(0);
 129   fKShort = new TH1F("hKShort","HLT:  K_{s}^{0} inv mass",80,0.4,.6);
 130   fKShort->SetFillColor(kGreen);
 131   fKShort->SetStats(0);
 132   fLambda = new TH1F("hLambda","HLT:  #Lambda^{0} inv mass",50,1.0,1.36);
 133   fLambda->SetFillColor(kGreen);
 134   fLambda->SetStats(0);
 135
 136   fPi0 = new TH1F("hPi0","HLT:  #Pi^{0} inv mass",50,0.0,0.38);
 137   fPi0->SetFillColor(kGreen);
 138   fPi0->SetStats(0);
 139
 140   fAP = new TH2F("hAP","HLT:  Armenteros-Podolanski plot",60,-1.,1.,60,-0.02,0.3);
 141   fAP->SetMarkerStyle(8);
 142   fAP->SetMarkerSize(0.4);
 143   fAP->SetYTitle("p_{t}[GeV]");
 144   fAP->SetXTitle("(p^{+}_{L}-p^{-}_{L})/(p^{+}_{L}+p^{-}_{L})");
 145   fAP->SetStats(0);
 146   fAP->SetOption("CONT4 Z");
 147
 148   fGammaXY = new TH2F("hGammaXY","HLT:  #gamma conversions",100,-100.,100.,100,-100.,100.);
 149   fGammaXY->SetMarkerStyle(8);
 150   fGammaXY->SetMarkerSize(0.4);
 151   fGammaXY->SetYTitle("X[cm]");
 152   fGammaXY->SetXTitle("Y[cm]");
 153   fGammaXY->SetStats(0);
 154   fGammaXY->SetOption("COLZ");
 155
 156   fNEvents =0;
 157   fNGammas = 0;
 158   fNKShorts = 0;
 159   fNLambdas = 0;
 160   fNPi0s = 0;
 161
 162   // cuts:
 163   // [0] == 0   --- N clusters on each daughter track
 164   // [1] == 2.5 --- (daughter-primVtx)/sigma >= cut
 165   // [2] == 3.5 --- (v0 - primVtx)/sigma <= cut
 166   // [3] == 3.0 --- (decay_length)/sigma >= cut
 167   // [4] == 0.0 --- (decay_length)[cm]   >= cut
 168   // [5] == 300.0 --- (v0 radius)[cm]    <= cut
 169   // [6] == 3.5  --- (v0 mass - true value)/sigma <= cut (for identification)
 170   // [7] == 0.05 --- (v0 mass - true value)       <= cut (for identification)
 171
 172   fGammaCuts[0] = 0;
 173   fGammaCuts[1] = 2.5;
 174   fGammaCuts[2] = 3.5;
 175   fGammaCuts[3] = 3.0;
 176   fGammaCuts[4] = 0.0;
 177   fGammaCuts[5] = 300.0;
 178   fGammaCuts[6] = 3.5;
 179   fGammaCuts[7] = 0.05;
 180
 181   fAPCuts[0] = 60;
 182   fAPCuts[1] = 2.5;
 183   fAPCuts[2] = 3.5;
 184   fAPCuts[3] = 3.0;
 185   fAPCuts[4] = 0.0;
 186   fAPCuts[5] = 50.0;
 187   fAPCuts[6] = 4.0;
 188   fAPCuts[7] = 0.05;
 189
 190   fKsCuts[0] = 60;
 191   fKsCuts[1] = 2.5;
 192   fKsCuts[2] = 3.5;
 193   fKsCuts[3] = 3.0;
 194   fKsCuts[4] = 1.5;
 195   fKsCuts[5] = 50.0;
 196   fKsCuts[6] = 4.0;
 197   fKsCuts[7] = 0.015;
 198
 199   fLambdaCuts[0] = 60;   // [0] == 60  --- N clusters on each daughter track
 200   fLambdaCuts[1] = 3.0;  // [1] == 3.0 --- (daughter-primVtx)/sigma >= cut
 201   fLambdaCuts[2] = 3.0;  // [2] == 3.0 --- (v0 - primVtx)/sigma <= cut
 202   fLambdaCuts[3] = 3.5;  // [3] == 3.5 --- (decay_length)/sigma >= cut
 203   fLambdaCuts[4] = 4.0;  // [4] == 0.0 --- (decay_length)[cm]   >= cut
 204   fLambdaCuts[5] = 50.0; // [5] == 300.0 --- (v0 radius)[cm]    <= cut
 205   fLambdaCuts[6] = 4.0;  // [6] == 3.5  --- (v0 mass - true value)/sigma <= cut (for identification)
 206   fLambdaCuts[7] = 0.03; // [7] == 0.05 --- (v0 mass - true value)       <= cut (for identification)
 207
 208   int iResult=0;
 209   TString configuration="";
 210   TString argument="";
 211   for (int i=0; i<argc && iResult>=0; i++) {
 212     argument=argv[i];
 213     if (!configuration.IsNull()) configuration+=" ";
 214     configuration+=argument;
 215   }
 216
 217   if (!configuration.IsNull()) {
 218     iResult=Configure(configuration.Data());
 219   }
 220
 221   return iResult;
 222 }
 223
 224 int AliHLTV0HistoComponent::DoDeinit()
 225 {
 226   // see header file for class documentation
 227   delete fGamma;
 228   delete fKShort;
 229   delete fLambda;
 230   delete fAP;
 231   delete fGammaXY;
 232   fUID = 0;
 233   return 0;
 234 }
 235
 236 int AliHLTV0HistoComponent::DoEvent(const AliHLTComponentEventData& evtData, AliHLTComponentTriggerData& /*trigData*/)
 237 {
 238
 239   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
 240     return 0;
 241
 242   if( fUID == 0 ){
 243     TTimeStamp t;
 244     fUID = ( gSystem->GetPid() + t.GetNanoSec())*10 + evtData.fEventID;
 245   }
 246
 247   fNEvents++;
 248
 249   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
 250
 251     AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 252     event->GetStdContent();
 253     Int_t nV0 = event->GetNumberOfV0s();
 254     AliKFParticle::SetField( event->GetMagneticField() );
 255
 256     const double kKsMass = 0.49767;
 257     const double kLambdaMass = 1.11568;
 258     const double kPi0Mass = 0.13498;
 259
 260     std::vector<AliKFParticle> vGammas;
 261
 262     for (Int_t iv=0; iv<nV0; iv++) {
 263
 264       AliESDtrack *t1=event->GetTrack( event->GetV0(iv)->GetNindex());
 265       AliESDtrack *t2=event->GetTrack( event->GetV0(iv)->GetPindex());
 266
 267       AliKFParticle kf1( *t1, 11 );
 268       AliKFParticle kf2( *t2, 11 );
 269
 270       AliKFVertex primVtx( *event->GetPrimaryVertexTracks() );
 271       double dev1 = kf1.GetDeviationFromVertex( primVtx );
 272       double dev2 = kf2.GetDeviationFromVertex( primVtx );
 273
 274       AliKFParticle v0( kf1, kf2 );
 275       double devPrim = v0.GetDeviationFromVertex( primVtx );
 276       primVtx+=v0;
 277       v0.SetProductionVertex( primVtx );
 278
 279       Double_t length, sigmaLength;
 280       if( v0.GetDecayLength( length, sigmaLength ) ) continue;
 281
 282       double dx = v0.GetX()-primVtx.GetX();
 283       double dy = v0.GetY()-primVtx.GetY();
 284       double r = sqrt(dx*dx + dy*dy);
 285
 286       // AP plot
 287
 288       double pt=0, ap=0;
 289       {
 290         AliKFParticle kf01 = kf1, kf02 = kf2;
 291         kf01.SetProductionVertex(v0);
 292         kf02.SetProductionVertex(v0);
 293         kf01.TransportToProductionVertex();
 294         kf02.TransportToProductionVertex();
 295         double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
 296         double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
 297         double px = px1+px2, py = py1+py2, pz = pz1+pz2;
 298         double p = sqrt(px*px+py*py+pz*pz);
 299         double l1 = (px*px1 + py*py1 + pz*pz1)/p;
 300         double l2 = (px*px2 + py*py2 + pz*pz2)/p;
 301         pt = sqrt(px1*px1+py1*py1+pz1*pz1 - l1*l1);
 302         ap = (l2-l1)/(l1+l2);
 303       }
 304
 305       if(
 306          t1->GetTPCNcls()>=fAPCuts[0]
 307          && t2->GetTPCNcls()>=fAPCuts[0]
 308          && dev1>=fAPCuts[1]
 309          && dev2>=fAPCuts[1]
 310          && devPrim <= fAPCuts[2]
 311          && length >= fAPCuts[3]*sigmaLength
 312          && length >= fAPCuts[4]
 313          && r <= fAPCuts[5]
 314          ){
 315         if( fAP ) fAP->Fill( ap, pt );
 316       }
 317
 318       // Gamma finder
 319
 320       bool isGamma = 0;
 321
 322       if(
 323          t1->GetTPCNcls()>=fGammaCuts[0]
 324          && t2->GetTPCNcls()>=fGammaCuts[0]
 325          && dev1>=fGammaCuts[1]
 326          && dev2>=fGammaCuts[1]
 327          && devPrim <= fGammaCuts[2]
 328          && length >= fGammaCuts[3]*sigmaLength
 329          && length >= fGammaCuts[4]
 330          && r <= fGammaCuts[5]
 331          ){
 332         double mass, error;
 333         v0.GetMass(mass,error);
 334         if( fGamma ) fGamma->Fill( mass );
 335
 336         if( TMath::Abs(mass)<=fGammaCuts[6]*error || TMath::Abs(mass)<=fGammaCuts[7] ){
 337           AliKFParticle gamma = v0;
 338           gamma.SetMassConstraint(0);
 339           if( fGammaXY
 340               &&  t1->GetTPCNcls()>=60
 341               && t2->GetTPCNcls()>=60
 342               ) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
 343           isGamma = 1;
 344           fNGammas++;
 345           vGammas.push_back( gamma );
 346         }
 347       }
 348
 349       if( isGamma ) continue;
 350
 351
 352       // KShort finder
 353
 354       bool isKs = 0;
 355
 356       if(
 357          t1->GetTPCNcls()>=fKsCuts[0]
 358          && t2->GetTPCNcls()>=fKsCuts[0]
 359          && dev1>=fKsCuts[1]
 360          && dev2>=fKsCuts[1]
 361          && devPrim <= fKsCuts[2]
 362          && length >= fKsCuts[3]*sigmaLength
 363          && length >= fKsCuts[4]
 364          && r <= fKsCuts[5]
 365          ){
 366
 367         AliKFParticle piN( *t1, 211 );
 368         AliKFParticle piP( *t2, 211 );
 369
 370         AliKFParticle Ks( piN, piP );
 371         Ks.SetProductionVertex( primVtx );
 372
 373         double mass, error;
 374         Ks.GetMass( mass, error);
 375         if( fKShort ) fKShort->Fill( mass );
 376         if( TMath::Abs( mass - kKsMass )<=fKsCuts[6]*error || TMath::Abs( mass - kKsMass )<=fKsCuts[7] ){
 377           isKs = 1;
 378           fNKShorts++;
 379         }
 380       }
 381
 382       if( isKs ) continue;
 383
 384       // Lambda finder
 385
 386       if(
 387          t1->GetTPCNcls()>=fLambdaCuts[0]
 388          && t2->GetTPCNcls()>=fLambdaCuts[0]
 389          && dev1>=fLambdaCuts[1]
 390          && dev2>=fLambdaCuts[1]
 391          && devPrim <= fLambdaCuts[2]
 392          && length >= fLambdaCuts[3]*sigmaLength
 393          && length >= fLambdaCuts[4]
 394          && r <= fLambdaCuts[5]
 395          && TMath::Abs( ap )>.4
 396          ){
 397
 398         AliKFParticle kP, kpi;
 399         if( ap<0 ){
 400           kP = AliKFParticle( *t2, 2212 );
 401           kpi = AliKFParticle( *t1, 211 );
 402         } else {
 403           kP = AliKFParticle( *t1, 2212 );
 404           kpi = AliKFParticle( *t2, 211 );
 405         }
 406
 407         AliKFParticle lambda = AliKFParticle( kP, kpi );
 408         lambda.SetProductionVertex( primVtx );
 409         double mass, error;
 410         lambda.GetMass( mass, error);
 411         if( fLambda ) fLambda->Fill( mass );
 412         if( TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[6]*error || TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[7] ){
 413           fNLambdas++;
 414         }
 415       }
 416
 417     }// V0's
 418
 419
 420     // Pi0 finder
 421
 422     for(UInt_t g1=0;g1<vGammas.size();g1++){
 423       for(UInt_t g2=g1+1;g2<vGammas.size();g2++){
 424         AliKFParticle pi0(vGammas.at(g1),vGammas.at(g2));
 425         double mass, error;
 426         pi0.GetMass(mass,error);
 427         fPi0->Fill(mass);
 428         if( TMath::Abs( mass - kPi0Mass )<=0.03 ){
 429           fNPi0s++;
 430         }
 431       }
 432     }
 433
 434
 435     if( fGamma ) PushBack( (TObject*) fGamma, kAliHLTDataTypeHistogram,fUID);
 436
 437     if( fKShort ) PushBack( (TObject*) fKShort, kAliHLTDataTypeHistogram,fUID);
 438
 439     if( fLambda ) PushBack( (TObject*) fLambda, kAliHLTDataTypeHistogram, fUID);
 440
 441     if( fPi0 ) PushBack( (TObject*) fPi0, kAliHLTDataTypeHistogram, fUID);
 442
 443     if( fAP ) PushBack( (TObject*) fAP, kAliHLTDataTypeHistogram, fUID);
 444
 445     if( fGammaXY ) PushBack( (TObject*) fGammaXY, kAliHLTDataTypeHistogram, fUID);
 446   }
 447   if( fNPi0s>0 ){
 448     HLTInfo("Found %d Gammas, %d KShorts, %d Lambdas and %d Pi0's in %d events", fNGammas, fNKShorts, fNLambdas, fNPi0s, fNEvents );
 449   }
 450   else HLTInfo("Found %d Gammas, %d KShorts and %d Lambdas in %d events", fNGammas, fNKShorts, fNLambdas, fNEvents );
 451
 452   return 0;
 453 }
 454
 455 int AliHLTV0HistoComponent::Configure(const char* arguments)
 456 {
 457
 458   int iResult=0;
 459   if (!arguments) return iResult;
 460
 461   TString allArgs=arguments;
 462   TString argument;
 463
 464   TObjArray* pTokens=allArgs.Tokenize(" ");
 465   int bMissingParam=0;
 466
 467   if (pTokens) {
 468     for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 469       argument=((TObjString*)pTokens->At(i))->GetString();
 470       if (argument.IsNull()) continue;
 471
 472       if (argument.CompareTo("-cutsGamma")==0) {
 473         TString spar = "";
 474         for( int j=0; j<8; j++ ){
 475           if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 476           spar+=" ";
 477           spar+=((TObjString*)pTokens->At(i))->GetString();
 478           fGammaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
 479         }
 480         if( !bMissingParam ){
 481           HLTInfo("Gamma cuts are set to: %s", spar.Data());
 482           continue;
 483         }
 484       } else if (argument.CompareTo("-cutsAP")==0) {
 485         TString spar = "";
 486         for( int j=0; j<8; j++ ){
 487           if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 488           spar+=" ";
 489           spar+=((TObjString*)pTokens->At(i))->GetString();
 490           fAPCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
 491         }
 492         if( !bMissingParam ){
 493           HLTInfo("AP cuts are set to: %s", spar.Data());
 494           continue;
 495         }
 496       }
 497       else if (argument.CompareTo("-cutsKs")==0) {
 498         TString spar = "";
 499         for( int j=0; j<8; j++ ){
 500           if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 501           spar+=" ";
 502           spar+=((TObjString*)pTokens->At(i))->GetString();
 503           fKsCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
 504         }
 505         if( !bMissingParam ){
 506           HLTInfo("KShort cuts are set to: %s", spar.Data());
 507           continue;
 508         }
 509       } else if (argument.CompareTo("-cutsLambda")==0) {
 510         TString spar = "";
 511         for( int j=0; j<8; j++ ){
 512           if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 513           spar+=" ";
 514           spar+=((TObjString*)pTokens->At(i))->GetString();
 515           fLambdaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
 516         }
 517         if( !bMissingParam ){
 518           HLTInfo("Lampda cuts are set to: %s", spar.Data());
 519           continue;
 520         }
 521       }else {
 522         HLTError("unknown argument %s", argument.Data());
 523         iResult=-EINVAL;
 524         break;
 525       }
 526     }
 527     delete pTokens;
 528   }
 529
 530   return iResult;
 531 }
 532
 533 int AliHLTV0HistoComponent::Reconfigure(const char* cdbEntry, const char* chainId)
 534 {
 535   // see header file for class documentation
 536   int iResult=0;
 537   const char* path="HLT/ConfigHLT/V0Histo";
 538   const char* defaultNotify="";
 539   if (cdbEntry) {
 540     path=cdbEntry;
 541     defaultNotify=" (default)";
 542   }
 543   if (path) {
 544     HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
 545     AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
 546     if (pEntry) {
 547       TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
 548       if (pString) {
 549         HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
 550         iResult=Configure(pString->GetString().Data());
 551       } else {
 552         HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
 553       }
 554     } else {
 555       HLTError("can not fetch object \"%s\" from CDB", path);
 556     }
 557   }
 558
 559   return iResult;
 560 }