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