PWG3/hfe/AliAnalysisTaskCheckV0tender.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Task fir checking the performance of the V0 tender
  17 //
  18 //
  19 // Authors
  20 //   Matus Kalisky <matus.kalisky@cern.ch>
  21 //
  22
  23 #include <sstream>
  24
  25 #include <TH1F.h>
  26 #include <TList.h>
  27
  28 #include "AliAnalysisManager.h"
  29 #include "AliMCEventHandler.h"
  30 #include "AliESDInputHandler.h"
  31 #include "AliESDv0.h"
  32 #include "AliESDtrack.h"
  33 #include "AliMCParticle.h"
  34 #include "AliMCEvent.h"
  35
  36 #include "AliHFEtools.h"
  37 #include "AliHFEcollection.h"
  38
  39 #include "AliAnalysisTaskCheckV0tender.h"
  40
  41 ClassImp(AliAnalysisTaskCheckV0tender)
  42
  43 //__________________________________________________________
  44 AliAnalysisTaskCheckV0tender::AliAnalysisTaskCheckV0tender():
  45   AliAnalysisTaskSE("CheckV0tenderTask")
  46   , fOutput(0x0)
  47   , fColl(0x0)
  48   , fCollMC(0x0)
  49   , fEvents(0x0)
  50 {
  51   //
  52   // Default Constructor
  53   //
  54 }
  55 //__________________________________________________________
  56 AliAnalysisTaskCheckV0tender::AliAnalysisTaskCheckV0tender(const Char_t *name):
  57   AliAnalysisTaskSE(name)
  58   , fOutput(0x0)
  59   , fColl(0x0)
  60   , fCollMC(0x0)
  61   , fEvents(0x0)
  62 {
  63   //
  64   // Default Constructor
  65   //
  66   DefineOutput(1, TH1F::Class());
  67   DefineOutput(2, TList::Class());
  68
  69 }
  70 //__________________________________________________________
  71 AliAnalysisTaskCheckV0tender::~AliAnalysisTaskCheckV0tender(){
  72   //
  73   // Destructor
  74   //
  75
  76   if (fOutput) delete fOutput;
  77   //if (fColl) delete fColl;
  78   //if (fCollMC) delete fCollMC;
  79   if (fEvents) delete fEvents;
  80 }
  81 //__________________________________________________________
  82 void AliAnalysisTaskCheckV0tender::UserCreateOutputObjects(){
  83   //
  84   // prepare output objects
  85   //
  86
  87   fOutput = new TList();
  88   fOutput->SetOwner();
  89   // Counter for number of events
  90   fEvents = new TH1I("nEvents", "NumberOfEvents", 1, 1, 2);
  91
  92   fColl = new AliHFEcollection("V0_QA", "tender V0s for data");
  93   fCollMC = new AliHFEcollection("V0_MC_QA", "tender V0s for MC");
  94
  95   //
  96   // Data histos
  97   //
  98
  99   // total number of tagged V0s
 100   fColl->CreateTH1F("h_NumberOf_V0s", "Number of tagged V0s; type; counts", 4, -0.5, 3.5);
 101   // pT spectra of the tagged V0s
 102   fColl->CreateTH1F("h_Gamma_pt", "p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 103   fColl->CreateTH1F("h_K0_pt", "p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 104   fColl->CreateTH1F("h_Lambda_pt", "p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 105   fColl->CreateTH1F("h_ALambda_pt", "p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 106   // invariant mass of the V0s
 107   fColl->CreateTH1F("h_Gamma_mass", "Inv. Mass of the ", 100, 0., 0.1);
 108   fColl->CreateTH1F("h_K0_mass", "Inv. Mass of the ", 100, 0.45, 0.55);
 109   fColl->CreateTH1F("h_Lambda_mass", "Inv. Mass of the ", 100, 1.08, 1.14);
 110   fColl->CreateTH1F("h_ALambda_mass", "Inv. Mass of the ", 100, 1.08, 1.14);
 111
 112   // total number of tagged daughter particles (should correlate with number of V0s !
 113   fColl->CreateTH1F("h_NumberOfDaughters", "Number of tagged daughters; type; counts", 3, -0.5, 2.5);
 114   // pT spectra of tagged daughter particles
 115   fColl->CreateTH1F("h_Electron_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 116   fColl->CreateTH1F("h_Pion_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 117   fColl->CreateTH1F("h_Proton_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 118
 119   //
 120   // MC histos
 121   //
 122
 123   // pT spectra of the tagged V0s
 124   fCollMC->CreateTH1F("h_Gamma_pt_S", "MC-S: p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 125   fCollMC->CreateTH1F("h_K0_pt_S", "MC-S: p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 126   fCollMC->CreateTH1F("h_Lambda_pt_S", "MC-S: p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 127   fCollMC->CreateTH1F("h_ALambda_pt_S", "MC-S: p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 128   fCollMC->CreateTH1F("h_Gamma_pt_B", "MC-B: p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 129   fCollMC->CreateTH1F("h_K0_pt_B", "MC-B: p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 130   fCollMC->CreateTH1F("h_Lambda_pt_B", "MC-B: p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 131   fCollMC->CreateTH1F("h_ALambda_pt_B", "MC-B: p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 132   // invariant mass of the V0s
 133   fCollMC->CreateTH1F("h_Gamma_mass_S", "MC-S: Inv. Mass of the gamma; m (GeV/c^{2}); counts", 100, 0., 0.1);
 134   fCollMC->CreateTH1F("h_K0_mass_S", "MC-S: Inv. Mass of the K0; m (GeV/c^{2}); counts", 100, 0.45, 0.55);
 135   fCollMC->CreateTH1F("h_Lambda_mass_S", "MC-S: Inv. Mass of the Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
 136   fCollMC->CreateTH1F("h_ALambda_mass_S", "MC-S: Inv. Mass of the A-Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
 137   fCollMC->CreateTH1F("h_Gamma_mass_B", "MC-B: Inv. Mass of the gamma; m (GeV/c^{2}); counts", 100, 0., 0.1);
 138   fCollMC->CreateTH1F("h_K0_mass_B", "MC-B: Inv. Mass of the K0; m (GeV/c^{2}); counts", 100, 0.45, 0.55);
 139   fCollMC->CreateTH1F("h_Lambda_mass_B", "MC-B: Inv. Mass of the Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
 140   fCollMC->CreateTH1F("h_ALambda_mass_B", "MC-B: Inv. Mass of the A-Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
 141   // pT spectra of tagged daughter particles
 142   fCollMC->CreateTH1F("h_Electron_pt_S", "MC-S: p_{T} spectrum of tagged electrons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 143   fCollMC->CreateTH1F("h_Pion_pt_S", "MC-S: p_{T} spectrum of tagged pions; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 144   fCollMC->CreateTH1F("h_Proton_pt_S", "MC-S: p_{T} spectrum of tagged protons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 145   fCollMC->CreateTH1F("h_Electron_pt_B", "MC-B: p_{T} spectrum of tagged electrons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 146   fCollMC->CreateTH1F("h_Pion_pt_B", "MC-B: p_{T} spectrum of tagged pions; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 147   fCollMC->CreateTH1F("h_Proton_pt_B", "MC-B: p_{T} spectrum of tagged protons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
 148
 149
 150   TList *tmp = fColl->GetList();
 151   tmp->SetName(fColl->GetName());
 152   fOutput->Add(tmp);
 153   tmp = 0x0;
 154   tmp = fCollMC->GetList();
 155   tmp->SetName(fCollMC->GetName());
 156   fOutput->Add(tmp);
 157
 158
 159 }
 160 //__________________________________________________________
 161 void AliAnalysisTaskCheckV0tender::UserExec(Option_t *){
 162   //
 163   // Event Loop
 164   //
 165
 166   AliMCEventHandler* mcHandler = (dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()));
 167   //AliESDInputHandler *inh = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 168   //AliESDpid *workingPID = NULL;
 169   //if(inh && (workingPID = inh->GetESDpid())) workingPID = inh->GetESDpid();
 170   //else workingPID = AliHFEtools::GetDefaultPID(mcHandler ? kTRUE : kFALSE);
 171
 172   // check the MC data
 173   if(fMCEvent && !mcHandler ) return;
 174   if(fMCEvent &&  !mcHandler->InitOk() ) return;
 175   if(fMCEvent &&  !mcHandler->TreeK() ) return;
 176   if(fMCEvent &&  !mcHandler->TreeTR() ) return;
 177
 178   ProcessV0s();
 179   ProcessDaughters();
 180
 181   if(fMCEvent){
 182     ProcessV0sMC();
 183     ProcessDaughtersMC();
 184   }
 185
 186   fEvents->Fill(1.1);
 187   PostData(1, fEvents);
 188   PostData(2, fOutput);
 189 }
 190 //__________________________________________________________
 191 void AliAnalysisTaskCheckV0tender::Terminate(Option_t *){
 192   //
 193   // Do Post Processing
 194   //
 195 }
 196 //__________________________________________________________
 197 void AliAnalysisTaskCheckV0tender::ProcessV0s(){
 198   //
 199   // loop over the V0s and extract the information about
 200   // the V0s tagged by the V0 tender
 201   //
 202
 203   const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
 204   TString name;
 205
 206   Int_t nV0s = fInputEvent->GetNumberOfV0s();
 207   for(Int_t i=0; i<nV0s; ++i){
 208     AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(i);
 209     if(!esdV0) continue;
 210     if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 211     Int_t pid = GetTenderPidV0(esdV0);
 212     //if(pid < 0) continue;
 213     pid = 2;
 214     fColl->Fill("h_NumberOf_V0s", pid);
 215     Float_t pT = esdV0->Pt();
 216     name = "h_" + type[pid] + "_pt";
 217     fColl->Fill(name, pT);
 218     Float_t mass = MassV0(esdV0, pid);
 219     name = "h_" + type[pid] + "_mass";
 220     //printf(" -D: name: %s \n", name.Data());
 221     fColl->Fill(name, mass);
 222   }
 223
 224 }
 225 //__________________________________________________________
 226 void AliAnalysisTaskCheckV0tender::ProcessDaughters(){
 227   //
 228   // produce some check plots for V0 tender tagged single tracks
 229   //
 230
 231   const TString type[3] = {"Electron", "Pion", "Proton"};
 232   TString name;
 233
 234   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 235   for(Int_t i=0; i<nTracks; ++i){
 236     AliESDtrack *track = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(i));
 237     if(!track) continue;
 238     Int_t  pid = GetTenderPidDaughter(track);
 239     if(pid < 0) continue;
 240     fColl->Fill("h_NumberOfDaughters", pid*1.0);
 241     Float_t pT = track->Pt();
 242     name = "h_" + type[pid] + "_pt";
 243     fColl->Fill(name, pT);
 244
 245   }
 246 }
 247 //__________________________________________________________
 248 void AliAnalysisTaskCheckV0tender::ProcessV0sMC(){
 249   //
 250   // check all V0tender selected V0 on their true identity
 251   //
 252
 253   const Int_t pid2pdg[4] = {22, 310, 3122, -3122};
 254   const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
 255   TString name;
 256   Int_t nV0s = fInputEvent->GetNumberOfV0s();
 257
 258   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 259
 260   // V0 loop
 261   for(Int_t i=0; i<nV0s; ++i){
 262     Bool_t id = kFALSE;
 263     AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(i);
 264     if(!esdV0) continue;
 265     if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 266     Int_t pid = GetTenderPidV0(esdV0);
 267     if(pid < 0) continue;
 268
 269     // both ESD daughtr tracks
 270     Int_t iN, iP;
 271     iN = iP = -1;
 272     iP = esdV0->GetPindex();
 273     iN = esdV0->GetNindex();
 274     if(iN < 0 || iP < 0) continue;
 275     if(iN >= nTracks || iP >= nTracks) continue;
 276     AliESDtrack *dP, *dN;
 277     dP = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iP));
 278     dN = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iN));
 279     if(!dN || !dP) continue;
 280
 281     // MC labels of the daughter tracks
 282     Int_t lN, lP;
 283     lN = dN->GetLabel();
 284     lP = dP->GetLabel();
 285     if(lN < 0 || lP < 0) continue;
 286
 287     // MC daughter particles
 288     AliMCParticle *mcP, *mcN;
 289     mcP = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lP));
 290     mcN = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lN));
 291     if(!mcP || !mcN) continue;
 292
 293     // labels of the mother particles
 294     Int_t lPm, lNm;
 295     lPm = mcP->GetMother();
 296     lNm = mcN->GetMother();
 297     if(lPm < 0) continue;
 298     AliMCParticle *m = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lPm));
 299     if(!m) continue;
 300     Int_t pdg = m->PdgCode();
 301     if((lPm == lNm) && (pdg == pid2pdg[pid])) id = kTRUE;
 302
 303     if(id) name = "h_" + type[pid] + "_pt_S";
 304     else name = "h_" + type[pid] + "_pt_B";
 305     Float_t pT = esdV0->Pt();
 306     fCollMC->Fill(name, pT);
 307
 308     if(id) name = "h_" + type[pid] + "_mass_S";
 309     else name = "h_" + type[pid] + "_mass_B";
 310     Float_t mass = MassV0(esdV0, pid);
 311     fCollMC->Fill(name, mass);
 312
 313    }
 314 }
 315 //__________________________________________________________
 316 void AliAnalysisTaskCheckV0tender::ProcessDaughtersMC(){
 317   //
 318   // check the identity of the V0tender selected V0 daughters
 319   // !!! for positive check only the true identity plays a role here,
 320   // not the true V0 mother identity (e.g. selected electron could come
 321   // from primary vertex or pion dalitz deca or true gamma conversion) !!!
 322   //
 323
 324   const Int_t pid2pdg [3] = {11, 211, 2212};
 325   const TString type[3] = {"Electron", "Pion", "Proton"};
 326   TString name;
 327
 328
 329   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 330   for(Int_t i=0; i<nTracks; ++i){
 331     Bool_t id = kFALSE;
 332     AliESDtrack *track = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(i));
 333     if(!track) continue;
 334     Int_t  pid = GetTenderPidDaughter(track);
 335     if(pid < 0) continue;
 336     Float_t pT = track->Pt();
 337     Int_t label = track->GetLabel();
 338     if(label < 0) continue;
 339     AliMCParticle *mcp = 0x0;
 340     mcp = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
 341     if(!mcp) continue;
 342     Int_t pdg = TMath::Abs(mcp->PdgCode());
 343     if(pdg == pid2pdg[pid]) id = kTRUE;
 344     if(id) name = "h_" + type[pid] + "_pt_S";
 345     else name = "h_" + type[pid] + "_pt_B";
 346     fCollMC->Fill(name, pT);
 347   }
 348 }
 349 //__________________________________________________________
 350 Int_t AliAnalysisTaskCheckV0tender::GetTenderPidV0(AliESDv0 * const v0){
 351   //
 352   // retrieve the PID nformation stored in the status flags by the train tender
 353   //
 354   Int_t pid = -1;
 355   if(!v0){
 356     return pid;
 357   }
 358   Int_t nTimes = 0;
 359   if(v0->TestBit(BIT(14))){
 360     pid = 0;
 361     nTimes++;
 362   }
 363   if(v0->TestBit(BIT(15))){
 364     pid = 1;
 365     nTimes++;
 366   }
 367   if(v0->TestBit(BIT(16))){
 368     pid = 2;
 369     nTimes++;
 370   }
 371   if(v0->TestBit(BIT(17))){
 372     pid = 3;
 373     nTimes++;
 374   }
 375   if(nTimes > 1){
 376     AliWarning("V0 track labeled multiple times by the V0 tender");
 377     pid = -1;
 378   }
 379
 380   //printf(" -D: pid: %i \n", pid);
 381
 382   return pid;
 383 }
 384 //__________________________________________________________
 385 Int_t AliAnalysisTaskCheckV0tender::GetTenderPidDaughter(AliESDtrack * const track){
 386   //
 387   // retrieve the PID nformation stored in the status flags by the train tender
 388   //
 389
 390   Int_t pid = -1;
 391   if(!track){
 392     return pid;
 393   }
 394   Int_t nTimes = 0;
 395   if(track->TestBit(BIT(14))){
 396     pid = 0;
 397     nTimes++;
 398   }
 399   if(track->TestBit(BIT(15))){
 400     pid = 1;
 401     nTimes++;
 402   }
 403   if(track->TestBit(BIT(16))){
 404     pid = 2;
 405     nTimes++;
 406   }
 407   if(nTimes > 1){
 408     AliWarning("V0 track labeled multiple times by the V0 tender");
 409     pid = -1;
 410   }
 411   return pid;
 412 }
 413 //__________________________________________________________
 414 Float_t AliAnalysisTaskCheckV0tender::MassV0(AliESDv0 * const v0, Int_t id){
 415   //
 416   // Get the V0 effective mass
 417   //
 418
 419   Float_t mass = -0.1;
 420   Bool_t sign = CheckSigns(v0);
 421   if(0 == id){
 422     mass = v0->GetEffMass(0, 0);
 423   }
 424   else if(1 == id){
 425     mass = v0->GetEffMass(2, 2);
 426   }
 427   else if(2 == id){
 428     mass = (sign) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
 429   }
 430   else if(3 == id){
 431     mass = (sign) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
 432   }
 433   else{
 434     AliWarning(Form("Unrecognized V0 id: %i", id));
 435   }
 436
 437   return mass;
 438
 439 }
 440 //__________________________________________________________
 441 Bool_t AliAnalysisTaskCheckV0tender::CheckSigns(AliESDv0 * const v0){
 442   //
 443   // check wheter the sign was correctly applied to
 444   // V0 daughter tracks
 445   // This function should become obsolete once the V0 finder will be updated (fixed)
 446   //
 447
 448   Bool_t correct = kFALSE;
 449
 450   Int_t pIndex = 0, nIndex = 0;
 451   pIndex = v0->GetPindex();
 452   nIndex = v0->GetNindex();
 453
 454   AliESDtrack* d[2];
 455   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
 456   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
 457
 458   Int_t sign[2];
 459   sign[0] = (int)d[0]->GetSign();
 460   sign[1] = (int)d[1]->GetSign();
 461
 462   if(-1 == sign[0] && 1 == sign[1]){
 463     correct = kFALSE;
 464     //v0->SetIndex(0, pIndex);  // set the index of the negative v0 track
 465     //v0->SetIndex(1, nIndex);  // set the index of the positive v0 track
 466   }
 467   else{
 468     correct = kTRUE;
 469   }
 470
 471   //pIndex = v0->GetPindex();
 472   //nIndex = v0->GetNindex();
 473   //printf("-D2: P: %i, N: %i\n", pIndex, nIndex);
 474
 475   return correct;
 476 }