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