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