HLT/trigger/AliHLTD0Trigger.cxx

   1 // $Id: AliHLTD0Trigger.cxx
   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: Gaute Ovrebekk                                        *
   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   AliHLTD0Trigger.cxx
  19 /// @author Gaute Ovrebekk
  20 /// @date   2009-10-28
  21 /// @brief  HLT trigger component for D0->Kpi
  22
  23 // see header file for class documentation
  24 // or
  25 // refer to README to build package
  26 // or
  27 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  28
  29 #include "AliHLTD0Trigger.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDv0.h"
  32 #include "AliHLTTriggerDecision.h"
  33 #include "AliHLTDomainEntry.h"
  34 #include "AliHLTGlobalBarrelTrack.h"
  35 #include "TObjArray.h"
  36 #include "TObjString.h"
  37 #include "TDatabasePDG.h"
  38 #include "AliESDVertex.h"
  39 #include "TH1F.h"
  40 #include "AliHLTD0toKpi.h"
  41 #include "AliAODVertex.h"
  42 #include "AliESDVertex.h"
  43
  44 /** ROOT macro for the implementation of ROOT specific class methods */
  45 ClassImp(AliHLTD0Trigger)
  46
  47 AliHLTD0Trigger::AliHLTD0Trigger()
  48   : AliHLTTrigger()
  49   , fPtMin(0.0)
  50   , fdca(0.0)
  51   , finvMass(0.0)
  52   , fcosThetaStar(0.0)
  53   , fd0(0.0)
  54   , fd0d0(0.0)
  55   , fcosPoint(0.0)
  56   , fplothisto(false)
  57   , fUseV0(false)
  58   , mD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  59   , fD0mass(NULL)
  60   , fPos()
  61   , fNeg()
  62   , fd0calc(NULL)
  63   , ftwoTrackArray(NULL)
  64   , fTotalD0(0)
  65   , fVertex(NULL)
  66   , fField(0)
  67 {
  68
  69   // see header file for class documentation
  70   // or
  71   // refer to README to build package
  72   // or
  73   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  74
  75 }
  76
  77 const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
  78
  79 AliHLTD0Trigger::~AliHLTD0Trigger()
  80 {
  81   //if(fd0calc){delete fd0calc;}
  82   //if(fD0mass){delete fD0mass;}
  83   //if(ftwoTrackArray){delete ftwoTrackArray;}
  84   // see header file for class documentation
  85 }
  86
  87 const char* AliHLTD0Trigger::GetTriggerName() const
  88 {
  89   // see header file for class documentation
  90   return "D0Trigger";
  91 }
  92
  93 AliHLTComponent* AliHLTD0Trigger::Spawn()
  94 {
  95   // see header file for class documentation
  96   return new AliHLTD0Trigger;
  97 }
  98
  99 int AliHLTD0Trigger::DoTrigger()
 100 {
 101   // -- Iterator over Data Blocks --
 102   //const AliHLTComponentBlockData* iter = NULL;
 103
 104   if (!IsDataEvent()) return 0;
 105
 106   Int_t nD0=0;
 107   TString description;
 108
 109   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
 110     if(fUseV0){
 111       nD0=RecV0(iter);
 112     }
 113     else{
 114        AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 115        event->GetStdContent();
 116        fField = event->GetMagneticField();
 117        const AliESDVertex* pv = event->GetPrimaryVertexTracks();
 118        fVertex =  new AliESDVertex(*pv);
 119
 120        for(Int_t it=0;it<event->GetNumberOfTracks();it++){
 121          SingleTrackSelect(event->GetTrack(it));
 122        }
 123
 124        nD0=RecD0();
 125     }
 126   }
 127
 128   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS);
 129         iter != NULL; iter = GetNextInputObject() ) {
 130     fVertex = dynamic_cast<AliESDVertex*>(const_cast<TObject*>( iter ));
 131     if(!fVertex){
 132       HLTError("ITS SPD vertex object is corrupted");
 133       //iResult = -EINVAL;
 134     }
 135   }
 136
 137   for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
 138         iter != NULL; iter = GetNextInputBlock() ) {
 139     vector<AliHLTGlobalBarrelTrack> tracksVector;
 140     AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
 141
 142     fField = GetBz();
 143     for(UInt_t i=0;i<tracksVector.size();i++){
 144       SingleTrackSelect(&tracksVector[i]);
 145     }
 146     nD0=RecD0();
 147   }
 148
 149   fTotalD0+=nD0;
 150
 151   ftwoTrackArray->Clear();
 152   fPos.clear();
 153   fNeg.clear();
 154
 155   HLTInfo("Number of D0 found: %d",nD0);
 156   HLTInfo("Total Number of D0 found: %d",fTotalD0);
 157
 158   if(fplothisto){PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);}
 159
 160   //if (iResult>=0) {
 161   if (1) {
 162
 163     if (nD0>=1) {
 164       description.Form("Event contains %d D0(s)", nD0);
 165       SetDescription(description.Data());
 166       // Enable the central detectors for readout.
 167       GetReadoutList().Enable(
 168                               AliHLTReadoutList::kITSSPD |
 169                               AliHLTReadoutList::kITSSDD |
 170                               AliHLTReadoutList::kITSSSD |
 171                               AliHLTReadoutList::kTPC |
 172                               AliHLTReadoutList::kTRD |
 173                               AliHLTReadoutList::kTOF |
 174                               AliHLTReadoutList::kHMPID |
 175                               AliHLTReadoutList::kPHOS
 176                               );
 177       // Add the available HLT information for readout too.
 178       GetTriggerDomain().Add("CLUSTERS", "TPC ");
 179       TriggerEvent(true);
 180       return 0;
 181     }
 182     description.Form("No D0");
 183   } else {
 184     description.Form("No input blocks found");
 185   }
 186   SetDescription(description.Data());
 187   TriggerEvent(false);
 188   //return iResult;
 189   return 0;
 190 }
 191
 192 int AliHLTD0Trigger::DoInit(int argc, const char** argv)
 193 {
 194
 195   fd0calc = new AliHLTD0toKpi();
 196   ftwoTrackArray = new TObjArray(2);
 197
 198   fplothisto=false;
 199   // see header file for class documentation
 200   fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 201   // first configure the default
 202   int iResult=0;
 203   if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
 204
 205   // configure from the command line parameters if specified
 206   if (iResult>=0 && argc>0)
 207     iResult=ConfigureFromArgumentString(argc, argv);
 208   return iResult;
 209 }
 210
 211 int AliHLTD0Trigger::DoDeinit()
 212 {
 213   // see header file for class documentation
 214   if(fd0calc){delete fd0calc;}
 215   if(fD0mass){delete fD0mass;}
 216   if(ftwoTrackArray){delete ftwoTrackArray;}
 217   if(fVertex){delete fVertex;}
 218   return 0;
 219 }
 220
 221 int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
 222 {
 223   // see header file for class documentation
 224
 225   // configure from the specified antry or the default one
 226   const char* entry=cdbEntry;
 227   if (!entry || entry[0]==0) {
 228     entry=fgkOCDBEntry;
 229   }
 230
 231   return ConfigureFromCDBTObjString(entry);
 232 }
 233
 234 int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
 235 {
 236   // see header file for class documentation
 237   if (argc<=0) return 0;
 238   int i=0;
 239   TString argument=argv[i];
 240   // -minpt for decay
 241   if (argument.CompareTo("-pt")==0) {
 242     if (++i>=argc) return -EPROTO;
 243     argument=argv[i];
 244     fPtMin=argument.Atof();
 245     return 2;
 246   }
 247   // minimum dca for decay tracks
 248   if (argument.CompareTo("-dca")==0) {
 249     if (++i>=argc) return -EPROTO;
 250     argument=argv[i];
 251     fdca=argument.Atof();
 252     return 2;
 253   }
 254   // inv. mass half width.
 255   if (argument.CompareTo("-invmass")==0) {
 256     if (++i>=argc) return -EPROTO;
 257     argument=argv[i];
 258     finvMass=argument.Atof();
 259     return 2;
 260   }
 261
 262 // cos theta for decay angle
 263   if (argument.CompareTo("-costhetastar")==0) {
 264     if (++i>=argc) return -EPROTO;
 265     argument=argv[i];
 266     fcosThetaStar=argument.Atof();
 267     return 2;
 268   }
 269
 270   // impact parameter for decay
 271   if (argument.CompareTo("-d0")==0) {
 272     if (++i>=argc) return -EPROTO;
 273     argument=argv[i];
 274     fd0=argument.Atof();
 275     return 2;
 276   }
 277   // product of impact parameter
 278   if (argument.CompareTo("-d0d0")==0) {
 279     if (++i>=argc) return -EPROTO;
 280     argument=argv[i];
 281     fd0d0=argument.Atof();
 282     return 2;
 283   }
 284   // product of impact parameter
 285   if (argument.CompareTo("-cospoint")==0) {
 286     if (++i>=argc) return -EPROTO;
 287     argument=argv[i];
 288     fcosPoint=argument.Atof();
 289     return 2;
 290   }
 291   if (argument.CompareTo("-plothistogram")==0) {
 292     fplothisto=true;
 293     return 1;
 294   }
 295   if (argument.CompareTo("-usev0")==0) {
 296     fUseV0=true;
 297     return 1;
 298   }
 299
 300   // unknown argument
 301   return -EINVAL;
 302 }
 303
 304 void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t){
 305   // Offline har || på disse kuttene på de to henfallsproduktene
 306   Double_t pv[3];
 307   fVertex->GetXYZ(pv);
 308
 309   if(t->Pt()<fPtMin){return;}
 310   if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > fd0){return;}
 311
 312   if(t->Charge()>0){
 313     fPos.push_back(t);
 314   }
 315   else{
 316     fNeg.push_back(t);
 317   }
 318 }
 319
 320 Int_t AliHLTD0Trigger::RecD0(){
 321
 322   int nD0=0;
 323   Double_t D0,D0bar,xdummy,ydummy;
 324   Double_t d0[2];
 325   Double_t svpos[3];
 326   Double_t pvpos[3];
 327
 328   if(!fVertex){
 329     HLTError("No Vertex is set");
 330     return 0;
 331   }
 332   fVertex->GetXYZ(pvpos);
 333
 334   for(UInt_t ip=0;ip<fPos.size();ip++){
 335     AliExternalTrackParam *tP=fPos[ip];
 336     for(UInt_t in=0;in<fNeg.size();in++){
 337       AliExternalTrackParam *tN=fNeg[in];
 338
 339       tP->PropagateToDCA(fVertex,fField,kVeryBig);  //do I need this??????
 340       tN->PropagateToDCA(fVertex,fField,kVeryBig);
 341
 342       Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
 343       if(dcatPtN>fdca) { continue; }
 344
 345       ftwoTrackArray->AddAt(tP,0);
 346       ftwoTrackArray->AddAt(tN,1);
 347       AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,fField,fVertex);
 348       if(!vertexp1n1) {
 349         ftwoTrackArray->Clear();
 350         continue;
 351       }
 352
 353       vertexp1n1->GetXYZ(svpos);
 354
 355       tP->PropagateToDCA(vertexp1n1,fField,kVeryBig);
 356       tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
 357
 358       if((TMath::Abs(fd0calc->InvMass(tN,tP)-mD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-mD0PDG) > finvMass){continue;}
 359       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
 360       if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){continue;}
 361       d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
 362       d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
 363       if((d0[0]*d0[1]) > fd0d0){continue;}
 364       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
 365
 366       if(fplothisto){
 367         if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
 368           fD0mass->Fill(fd0calc->InvMass(tN,tP));
 369         }
 370         else{
 371           fD0mass->Fill(fd0calc->InvMass(tP,tN));
 372         }
 373       }
 374       nD0++;
 375       delete vertexp1n1;
 376     }
 377   }
 378   return nD0;
 379 }
 380
 381 Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
 382   int nD0=0;
 383   Double_t d0[2];
 384   Double_t D0,D0bar;
 385   Double_t svpos[3];
 386   Double_t pvpos[3];
 387
 388   AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 389   event->GetStdContent();
 390   Int_t nV0 = event->GetNumberOfV0s();
 391   Double_t field = event->GetMagneticField();
 392   const AliESDVertex* pv = event->GetPrimaryVertexTracks();
 393   pv->GetXYZ(pvpos);
 394
 395   for (Int_t iv=0; iv<nV0; iv++) {
 396
 397       AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
 398       AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());
 399
 400       if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
 401
 402       d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
 403       d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
 404
 405       if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
 406
 407       event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
 408
 409       if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
 410         HLTInfo("Tracks could not be propagated to secondary vertex");
 411         continue;
 412       }
 413
 414       Double_t tmp1, tmp2;
 415       if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
 416
 417       if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - mD0PDG) > finvMass){continue;}
 418       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
 419       if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
 420       if((d0[0]*d0[1]) > fd0d0){continue;}
 421       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
 422
 423       nD0++;
 424       if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
 425         fD0mass->Fill(fd0calc->InvMass(tN,tP));
 426       }
 427       else{
 428         fD0mass->Fill(fd0calc->InvMass(tP,tN));
 429       }
 430   }
 431   return nD0;
 432
 433 }