PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaCorrectionTask.h"
   4
   5 #include <TCanvas.h>
   6 #include <TChain.h>
   7 #include <TFile.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TProfile.h>
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TDatabasePDG.h>
  15
  16 #include <AliLog.h>
  17 #include <AliESDVertex.h>
  18 #include <AliESDEvent.h>
  19 #include <AliStack.h>
  20 #include <AliHeader.h>
  21 #include <AliGenEventHeader.h>
  22 #include <AliMultiplicity.h>
  23 #include <AliAnalysisManager.h>
  24 #include <AliMCEventHandler.h>
  25 #include <AliMCEvent.h>
  26 #include <AliESDInputHandler.h>
  27
  28 #include "AliESDtrackCuts.h"
  29 #include "AliPWG0Helper.h"
  30 #include "dNdEta/dNdEtaAnalysis.h"
  31 #include "dNdEta/AlidNdEtaCorrection.h"
  32
  33 ClassImp(AlidNdEtaCorrectionTask)
  34
  35 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
  36   AliAnalysisTask(),
  37   fESD(0),
  38   fOutput(0),
  39   fOption(),
  40   fAnalysisMode(AliPWG0Helper::kTPC),
  41   fTrigger(AliPWG0Helper::kMB1),
  42   fFillPhi(kFALSE),
  43   fDeltaPhiCut(-1),
  44   fSignMode(0),
  45   fOnlyPrimaries(kFALSE),
  46   fStatError(0),
  47   fEsdTrackCuts(0),
  48   fdNdEtaCorrection(0),
  49   fdNdEtaAnalysisMC(0),
  50   fdNdEtaAnalysisESD(0),
  51   fPIDParticles(0),
  52   fPIDTracks(0),
  53   fVertexCorrelation(0),
  54   fVertexCorrelationShift(0),
  55   fVertexProfile(0),
  56   fVertexShift(0),
  57   fVertexShiftNorm(0),
  58   fEtaCorrelation(0),
  59   fEtaCorrelationShift(0),
  60   fEtaProfile(0),
  61   fEtaResolution(0),
  62   fDeltaPhiCorrelation(0),
  63   fpTResolution(0),
  64   fEsdTrackCutsPrim(0),
  65   fEsdTrackCutsSec(0),
  66   fTemp1(0),
  67   fTemp2(0),
  68   fMultAll(0),
  69   fMultTr(0),
  70   fMultVtx(0),
  71   fEventStats(0)
  72 {
  73   //
  74   // Constructor. Initialization of pointers
  75   //
  76
  77   for (Int_t i=0; i<4; i++)
  78     fdNdEtaCorrectionSpecial[i] = 0;
  79
  80   for (Int_t i=0; i<8; i++)
  81     fDeltaPhi[i] = 0;
  82 }
  83
  84 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
  85   AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
  86   fESD(0),
  87   fOutput(0),
  88   fOption(opt),
  89   fAnalysisMode(AliPWG0Helper::kTPC),
  90   fTrigger(AliPWG0Helper::kMB1),
  91   fFillPhi(kFALSE),
  92   fDeltaPhiCut(0),
  93   fSignMode(0),
  94   fOnlyPrimaries(kFALSE),
  95   fStatError(0),
  96   fEsdTrackCuts(0),
  97   fdNdEtaCorrection(0),
  98   fdNdEtaAnalysisMC(0),
  99   fdNdEtaAnalysisESD(0),
 100   fPIDParticles(0),
 101   fPIDTracks(0),
 102   fVertexCorrelation(0),
 103   fVertexCorrelationShift(0),
 104   fVertexProfile(0),
 105   fVertexShift(0),
 106   fVertexShiftNorm(0),
 107   fEtaCorrelation(0),
 108   fEtaCorrelationShift(0),
 109   fEtaProfile(0),
 110   fEtaResolution(0),
 111   fDeltaPhiCorrelation(0),
 112   fpTResolution(0),
 113   fEsdTrackCutsPrim(0),
 114   fEsdTrackCutsSec(0),
 115   fTemp1(0),
 116   fTemp2(0),
 117   fMultAll(0),
 118   fMultTr(0),
 119   fMultVtx(0),
 120   fEventStats(0)
 121 {
 122   //
 123   // Constructor. Initialization of pointers
 124   //
 125
 126   // Define input and output slots here
 127   DefineInput(0, TChain::Class());
 128   DefineOutput(0, TList::Class());
 129
 130   for (Int_t i=0; i<4; i++)
 131     fdNdEtaCorrectionSpecial[i] = 0;
 132
 133   for (Int_t i=0; i<8; i++)
 134     fDeltaPhi[i] = 0;
 135 }
 136
 137 AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
 138 {
 139   //
 140   // Destructor
 141   //
 142
 143   // histograms are in the output list and deleted when the output
 144   // list is deleted by the TSelector dtor
 145
 146   if (fOutput) {
 147     delete fOutput;
 148     fOutput = 0;
 149   }
 150 }
 151
 152 //________________________________________________________________________
 153 void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
 154 {
 155   // Connect ESD
 156   // Called once
 157
 158   Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
 159
 160   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 161
 162   if (!esdH) {
 163     Printf("ERROR: Could not get ESDInputHandler");
 164   } else {
 165     fESD = esdH->GetEvent();
 166
 167     // Enable only the needed branches
 168     esdH->SetActiveBranches("AliESDHeader Vertex");
 169
 170     if (fAnalysisMode == AliPWG0Helper::kSPD)
 171       esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
 172
 173     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
 174       esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
 175     }
 176   }
 177
 178   // disable info messages of AliMCEvent (per event)
 179   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 180 }
 181
 182 void AlidNdEtaCorrectionTask::CreateOutputObjects()
 183 {
 184   // create result objects and add to output list
 185
 186   if (fOption.Contains("only-positive"))
 187   {
 188     Printf("INFO: Processing only positive particles.");
 189     fSignMode = 1;
 190   }
 191   else if (fOption.Contains("only-negative"))
 192   {
 193     Printf("INFO: Processing only negative particles.");
 194     fSignMode = -1;
 195   }
 196
 197   if (fOption.Contains("stat-error-1"))
 198   {
 199     Printf("INFO: Evaluation statistical errors. Mode: 1.");
 200     fStatError = 1;
 201   }
 202   else if (fOption.Contains("stat-error-2"))
 203   {
 204     Printf("INFO: Evaluation statistical errors. Mode: 2.");
 205     fStatError = 2;
 206   }
 207
 208   fOutput = new TList;
 209   fOutput->SetOwner();
 210
 211   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
 212   fOutput->Add(fdNdEtaCorrection);
 213
 214   fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
 215   fOutput->Add(fPIDParticles);
 216
 217   fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
 218   fOutput->Add(fPIDTracks);
 219
 220   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
 221   fOutput->Add(fdNdEtaAnalysisMC);
 222
 223   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
 224   fOutput->Add(fdNdEtaAnalysisESD);
 225
 226   if (fEsdTrackCuts)
 227   {
 228     fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
 229     fEsdTrackCutsSec  = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
 230     fOutput->Add(fEsdTrackCutsPrim);
 231     fOutput->Add(fEsdTrackCutsSec);
 232   }
 233
 234   if (fOption.Contains("process-types")) {
 235     fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
 236     fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
 237     fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
 238
 239     fOutput->Add(fdNdEtaCorrectionSpecial[0]);
 240     fOutput->Add(fdNdEtaCorrectionSpecial[1]);
 241     fOutput->Add(fdNdEtaCorrectionSpecial[2]);
 242   }
 243
 244   if (fOption.Contains("particle-species")) {
 245     fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_pi", "dndeta_correction_pi", fAnalysisMode);
 246     fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_K", "dndeta_correction_K", fAnalysisMode);
 247     fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_p", "dndeta_correction_p", fAnalysisMode);
 248     fdNdEtaCorrectionSpecial[3] = new AlidNdEtaCorrection("dndeta_correction_other", "dndeta_correction_other", fAnalysisMode);
 249
 250     for (Int_t i=0; i<4; i++)
 251       fOutput->Add(fdNdEtaCorrectionSpecial[i]);
 252   }
 253
 254
 255   /*
 256   fTemp1 = new TH2F("fTemp1", "fTemp1", 200, -0.08, 0.08, 200, -0.08, 0.08);
 257   fOutput->Add(fTemp1);
 258   fTemp2 = new TH1F("fTemp2", "fTemp2", 2000, -5, 5);
 259   fOutput->Add(fTemp2);
 260   */
 261
 262   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
 263   fOutput->Add(fVertexCorrelation);
 264   fVertexCorrelationShift = new TH2F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx", 120, -30, 30, 100, -1, 1);
 265   fOutput->Add(fVertexCorrelationShift);
 266   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
 267   fOutput->Add(fVertexProfile);
 268   fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
 269   fOutput->Add(fVertexShift);
 270   fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
 271   fOutput->Add(fVertexShiftNorm);
 272
 273   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
 274   fOutput->Add(fEtaCorrelation);
 275   fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
 276   fOutput->Add(fEtaCorrelationShift);
 277   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
 278   fOutput->Add(fEtaProfile);
 279   fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
 280   fOutput->Add(fEtaResolution);
 281
 282   fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2);
 283   fOutput->Add(fpTResolution);
 284
 285   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
 286   fOutput->Add(fMultAll);
 287   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
 288   fOutput->Add(fMultVtx);
 289   fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
 290   fOutput->Add(fMultTr);
 291
 292   for (Int_t i=0; i<8; i++)
 293   {
 294     fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi());
 295     fOutput->Add(fDeltaPhi[i]);
 296   }
 297
 298   fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 106, -5.5, 100.5, 4, -0.5, 3.5);
 299   fOutput->Add(fEventStats);
 300   fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5
 301   fEventStats->GetXaxis()->SetBinLabel(2, "NSD");  // x = -4
 302   fEventStats->GetXaxis()->SetBinLabel(3, "ND");   // x = -3
 303   fEventStats->GetXaxis()->SetBinLabel(4, "SD");   // x = -2
 304   fEventStats->GetXaxis()->SetBinLabel(5, "DD");   // x = -1
 305
 306   for (Int_t i=0; i<100; i++)
 307     fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i));
 308
 309   fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
 310   fEventStats->GetYaxis()->SetBinLabel(2, "trg");
 311   fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
 312   fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
 313
 314   if (fEsdTrackCuts)
 315   {
 316     fEsdTrackCuts->SetName("fEsdTrackCuts");
 317     fOutput->Add(fEsdTrackCuts);
 318   }
 319 }
 320
 321 void AlidNdEtaCorrectionTask::Exec(Option_t*)
 322 {
 323   // process the event
 324
 325   // Check prerequisites
 326   if (!fESD)
 327   {
 328     AliDebug(AliLog::kError, "ESD branch not available");
 329     return;
 330   }
 331
 332   if (fOnlyPrimaries)
 333     Printf("WARNING: Processing only primaries. For systematical studies only!");
 334
 335   if (fStatError > 0)
 336     Printf("WARNING: Statistical error evaluation active!");
 337
 338   // trigger definition
 339   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
 340
 341   if (!eventTriggered)
 342     Printf("No trigger");
 343
 344   // post the data already here
 345   PostData(0, fOutput);
 346
 347   // MC info
 348   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 349   if (!eventHandler) {
 350     Printf("ERROR: Could not retrieve MC event handler");
 351     return;
 352   }
 353
 354   AliMCEvent* mcEvent = eventHandler->MCEvent();
 355   if (!mcEvent) {
 356     Printf("ERROR: Could not retrieve MC event");
 357     return;
 358   }
 359
 360   AliStack* stack = mcEvent->Stack();
 361   if (!stack)
 362   {
 363     AliDebug(AliLog::kError, "Stack not available");
 364     return;
 365   }
 366
 367   AliHeader* header = mcEvent->Header();
 368   if (!header)
 369   {
 370     AliDebug(AliLog::kError, "Header not available");
 371     return;
 372   }
 373
 374   // get process type;
 375   Int_t processType = AliPWG0Helper::GetEventProcessType(header);
 376   AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
 377
 378   if (processType == AliPWG0Helper::kInvalidProcess)
 379     AliDebug(AliLog::kError, "Unknown process type.");
 380
 381   // get the MC vertex
 382   AliGenEventHeader* genHeader = header->GenEventHeader();
 383   TArrayF vtxMC(3);
 384   genHeader->PrimaryVertex(vtxMC);
 385
 386   // get the ESD vertex
 387   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 388   Bool_t eventVertex = kFALSE;
 389   if (vtxESD)
 390   {
 391     Double_t vtx[3];
 392     vtxESD->GetXYZ(vtx);
 393
 394     Double_t diff = vtxMC[2] - vtx[2];
 395     fVertexShift->Fill(diff);
 396     if (vtxESD->GetZRes() > 0)
 397         fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
 398
 399     if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 400     {
 401         vtxESD = 0;
 402     }
 403     else
 404     {
 405       eventVertex = kTRUE;
 406
 407       if (eventTriggered)
 408       {
 409         fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
 410         fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
 411         fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
 412       }
 413     }
 414   }
 415
 416   // fill process type
 417   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
 418   // INEL
 419   fEventStats->Fill(-5, biny);
 420   // NSD
 421   if (processType != AliPWG0Helper::kSD)
 422     fEventStats->Fill(-4, biny);
 423   // SD, ND, DD
 424   if (processType == AliPWG0Helper::kND)
 425     fEventStats->Fill(-3, biny);
 426   if (processType == AliPWG0Helper::kSD)
 427     fEventStats->Fill(-2, biny);
 428   if (processType == AliPWG0Helper::kDD)
 429     fEventStats->Fill(-1, biny);
 430
 431   // create list of (label, eta, pt) tuples
 432   Int_t inputCount = 0;
 433   Int_t* labelArr = 0;
 434   Int_t* labelArr2 = 0; // only for case of SPD
 435   Float_t* etaArr = 0;
 436   Float_t* thirdDimArr = 0;
 437   Float_t* deltaPhiArr = 0;
 438   if (fAnalysisMode == AliPWG0Helper::kSPD)
 439   {
 440     // get tracklets
 441     const AliMultiplicity* mult = fESD->GetMultiplicity();
 442     if (!mult)
 443     {
 444       AliDebug(AliLog::kError, "AliMultiplicity not available");
 445       return;
 446     }
 447
 448     labelArr = new Int_t[mult->GetNumberOfTracklets()];
 449     labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
 450     etaArr = new Float_t[mult->GetNumberOfTracklets()];
 451     thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
 452     deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
 453
 454     // get multiplicity from SPD tracklets
 455     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 456     {
 457       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 458
 459       Float_t phi = mult->GetPhi(i);
 460       if (phi < 0)
 461         phi += TMath::Pi() * 2;
 462       Float_t deltaPhi = mult->GetDeltaPhi(i);
 463
 464       if (TMath::Abs(deltaPhi) > 1)
 465         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 466
 467       if (fOnlyPrimaries)
 468         if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
 469           continue;
 470
 471       if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
 472         continue;
 473
 474       etaArr[inputCount] = mult->GetEta(i);
 475       labelArr[inputCount] = mult->GetLabel(i, 0);
 476       labelArr2[inputCount] = mult->GetLabel(i, 1);
 477       thirdDimArr[inputCount] = phi;
 478       deltaPhiArr[inputCount] = deltaPhi;
 479       ++inputCount;
 480     }
 481   }
 482   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 483   {
 484     if (!fEsdTrackCuts)
 485     {
 486       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 487       return;
 488     }
 489
 490     if (vtxESD)
 491     {
 492       // get multiplicity from ESD tracks
 493       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
 494       Int_t nGoodTracks = list->GetEntries();
 495
 496       Printf("Accepted %d tracks", nGoodTracks);
 497
 498       labelArr = new Int_t[nGoodTracks];
 499       labelArr2 = new Int_t[nGoodTracks];
 500       etaArr = new Float_t[nGoodTracks];
 501       thirdDimArr = new Float_t[nGoodTracks];
 502       deltaPhiArr = new Float_t[nGoodTracks];
 503
 504       // loop over esd tracks
 505       for (Int_t i=0; i<nGoodTracks; i++)
 506       {
 507         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 508         if (!esdTrack)
 509         {
 510           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 511           continue;
 512         }
 513
 514         // TODO fOnlyPrimaries not implemented for TPC
 515
 516         etaArr[inputCount] = esdTrack->Eta();
 517         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 518         labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
 519         thirdDimArr[inputCount] = esdTrack->Pt();
 520         deltaPhiArr[inputCount] = 0; // no delta phi for tracks
 521         ++inputCount;
 522       }
 523
 524       delete list;
 525
 526       if (eventTriggered)
 527       {
 528         // collect values for primaries and secondaries
 529         for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
 530         {
 531           AliESDtrack* track = 0;
 532
 533           if (fAnalysisMode == AliPWG0Helper::kTPC)
 534             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
 535           else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
 536             track = fESD->GetTrack(iTrack);
 537
 538           if (!track)
 539             continue;
 540
 541           Int_t label = TMath::Abs(track->GetLabel());
 542           if (!stack->Particle(label) || !stack->Particle(label)->GetPDG())
 543           {
 544             Printf("WARNING: No particle for %d", label);
 545             if (stack->Particle(label))
 546               stack->Particle(label)->Print();
 547             continue;
 548           }
 549
 550           if (stack->Particle(label)->GetPDG()->Charge() == 0)
 551             continue;
 552
 553           if (TMath::Abs(track->Eta()) < 1)
 554           {
 555             if (stack->IsPhysicalPrimary(label))
 556             {
 557               // primary
 558               if (fEsdTrackCutsPrim->AcceptTrack(track))
 559               {
 560                 if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
 561                 {
 562                   Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
 563                   Float_t b[2];
 564                   Float_t r[3];
 565                   track->GetImpactParameters(b, r);
 566                   Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
 567                   track->Print("");
 568                   if (vtxESD)
 569                     vtxESD->Print();
 570                 }
 571               }
 572             }
 573             else
 574             {
 575               // secondary
 576               fEsdTrackCutsSec->AcceptTrack(track);
 577             }
 578           }
 579
 580           // TODO mem leak in the continue statements above
 581           if (fAnalysisMode == AliPWG0Helper::kTPC)
 582             delete track;
 583         }
 584       }
 585     }
 586   }
 587   else
 588     return;
 589
 590   // loop over mc particles
 591   Int_t nPrim  = stack->GetNprimary();
 592   Int_t nAccepted = 0;
 593
 594   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 595   {
 596     //Printf("Getting particle %d", iMc);
 597     TParticle* particle = stack->Particle(iMc);
 598
 599     if (!particle)
 600       continue;
 601
 602     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 603     {
 604       //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0)
 605       //  fPIDParticles->Fill(particle->GetPdgCode());
 606       continue;
 607     }
 608
 609     if (SignOK(particle->GetPDG()) == kFALSE)
 610       continue;
 611
 612     if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
 613       fPIDParticles->Fill(particle->GetPdgCode());
 614
 615     Float_t eta = particle->Eta();
 616
 617     Float_t thirdDim = -1;
 618     if (fAnalysisMode == AliPWG0Helper::kSPD)
 619     {
 620       if (fFillPhi)
 621       {
 622         thirdDim = particle->Phi();
 623       }
 624       else
 625         thirdDim = inputCount;
 626     }
 627     else
 628       thirdDim = particle->Pt();
 629
 630     // calculate y
 631     //Float_t y = 0.5 * TMath::Log((particle->Energy() + particle->Pz()) / (particle->Energy() - particle->Pz()));
 632     //fTemp1->Fill(eta);
 633     //fTemp2->Fill(y);
 634
 635     fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 636
 637     if (fOption.Contains("process-types"))
 638     {
 639       // non diffractive
 640       if (processType==AliPWG0Helper::kND)
 641         fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 642
 643       // single diffractive
 644       if (processType==AliPWG0Helper::kSD)
 645         fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 646
 647       // double diffractive
 648       if (processType==AliPWG0Helper::kDD)
 649         fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 650     }
 651
 652     if (fOption.Contains("particle-species"))
 653     {
 654       Int_t id = -1;
 655       switch (TMath::Abs(particle->GetPdgCode()))
 656       {
 657         case 211:   id = 0; break;
 658         case 321:   id = 1; break;
 659         case 2212:  id = 2; break;
 660         default:    id = 3; break;
 661       }
 662       fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 663     }
 664
 665     if (eventTriggered)
 666       if (eventVertex)
 667         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim);
 668
 669     // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
 670     if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
 671       nAccepted++;
 672   }
 673
 674   fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
 675
 676   fMultAll->Fill(nAccepted);
 677   if (eventTriggered) {
 678     fMultTr->Fill(nAccepted);
 679     if (eventVertex)
 680       fMultVtx->Fill(nAccepted);
 681   }
 682
 683   Int_t processed = 0;
 684
 685   Bool_t* primCount = 0;
 686   if (fStatError > 0)
 687   {
 688     primCount = new Bool_t[nPrim];
 689     for (Int_t i=0; i<nPrim; ++i)
 690       primCount[i] = kFALSE;
 691   }
 692
 693   for (Int_t i=0; i<inputCount; ++i)
 694   {
 695     Int_t label = labelArr[i];
 696     Int_t label2 = labelArr2[i];
 697
 698     if (label < 0)
 699     {
 700       Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
 701       continue;
 702     }
 703
 704     TParticle* particle = stack->Particle(label);
 705     if (!particle)
 706     {
 707       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 708       continue;
 709     }
 710
 711     fPIDTracks->Fill(particle->GetPdgCode());
 712
 713     // find particle that is filled in the correction map
 714     // this should be the particle that has been reconstructed
 715     // for TPC this is clear and always identified by the track's label
 716     // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken:
 717     //  L1 L2 (P = primary, S = secondary)
 718     //   P P' : --> P
 719     //   P S  : --> P
 720     //   S P  : --> P
 721     //   S S' : --> S
 722
 723     if (label != label2 && label2 >= 0)
 724     {
 725       if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
 726       {
 727         particle = stack->Particle(label2);
 728         if (!particle)
 729         {
 730           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
 731           continue;
 732         }
 733       }
 734     }
 735
 736     if (eventTriggered && eventVertex)
 737     {
 738       if (SignOK(particle->GetPDG()) == kFALSE)
 739         continue;
 740
 741       processed++;
 742
 743       // resolutions
 744       fEtaResolution->Fill(particle->Eta() - etaArr[i]);
 745
 746       if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 747         if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
 748           fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
 749
 750       Float_t eta = -999;
 751       Float_t thirdDim = -1;
 752
 753       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
 754       // in case of primary the MC values are filled, otherwise (background) the reconstructed values
 755       if (label == label2 && firstIsPrim)
 756       {
 757         eta = particle->Eta();
 758
 759         if (fAnalysisMode == AliPWG0Helper::kSPD)
 760         {
 761           if (fFillPhi)
 762           {
 763             thirdDim = particle->Phi();
 764           }
 765           else
 766             thirdDim = inputCount;
 767         }
 768         else
 769           thirdDim = particle->Pt();
 770       }
 771       else
 772       {
 773         if (fAnalysisMode == AliPWG0Helper::kSPD && !fFillPhi)
 774         {
 775           thirdDim = inputCount;
 776         }
 777         else
 778           thirdDim = thirdDimArr[i];
 779
 780         eta = etaArr[i];
 781       }
 782
 783       Bool_t fillTrack = kTRUE;
 784
 785       // statistical error evaluation active?
 786       if (fStatError > 0)
 787       {
 788         Bool_t statErrorDecision = kFALSE;
 789
 790         // primary and not yet counted
 791         if (label == label2 && firstIsPrim && !primCount[label])
 792         {
 793           statErrorDecision = kTRUE;
 794           primCount[label] = kTRUE;
 795         }
 796
 797         // in case of 1 we count only unique primaries, in case of 2 all the rest
 798         if (fStatError == 1)
 799         {
 800           fillTrack = statErrorDecision;
 801         }
 802         else if (fStatError == 2)
 803           fillTrack = !statErrorDecision;
 804       }
 805
 806       if (fillTrack)
 807         fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 808
 809       // eta comparison for tracklets with the same label (others are background)
 810       if (label == label2)
 811       {
 812         fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
 813         fEtaCorrelation->Fill(etaArr[i], particle->Eta());
 814         fEtaCorrelationShift->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
 815       }
 816
 817       fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
 818
 819       if (fOption.Contains("process-types"))
 820       {
 821         // non diffractive
 822         if (processType == AliPWG0Helper::kND)
 823           fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 824
 825         // single diffractive
 826         if (processType == AliPWG0Helper::kSD)
 827           fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 828
 829         // double diffractive
 830         if (processType == AliPWG0Helper::kDD)
 831           fdNdEtaCorrectionSpecial[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 832       }
 833
 834       if (fOption.Contains("particle-species"))
 835       {
 836         // find mother first
 837         TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label);
 838
 839         Int_t id = -1;
 840         switch (TMath::Abs(mother->GetPdgCode()))
 841         {
 842           case 211:   id = 0; break;
 843           case 321:   id = 1; break;
 844           case 2212:  id = 2; break;
 845           default:    id = 3; break;
 846         }
 847         fdNdEtaCorrectionSpecial[id]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 848       }
 849
 850       // control histograms
 851       Int_t hist = -1;
 852       if (label == label2)
 853       {
 854         if (firstIsPrim)
 855         {
 856           hist = 0;
 857         }
 858         else
 859           hist = 1;
 860       }
 861       else if (label2 >= 0)
 862       {
 863         Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
 864         if (firstIsPrim && secondIsPrim)
 865         {
 866           hist = 2;
 867         }
 868         else if (firstIsPrim && !secondIsPrim)
 869         {
 870           hist = 3;
 871
 872           // check if secondary is caused by the primary or it is a fake combination
 873           //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
 874           Int_t mother = label2;
 875           while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
 876           {
 877             //Printf("  %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
 878             if (stack->Particle(mother)->GetMother(0) == label)
 879             {
 880               /*Printf("The untraceable decay was:");
 881               Printf("   from:");
 882               particle->Print();
 883               Printf("   to (status code %d):", stack->Particle(mother)->GetStatusCode());
 884               stack->Particle(mother)->Print();*/
 885               hist = 4;
 886             }
 887             mother = stack->Particle(mother)->GetMother(0);
 888           }
 889         }
 890         else if (!firstIsPrim && secondIsPrim)
 891         {
 892           hist = 5;
 893         }
 894         else if (!firstIsPrim && !secondIsPrim)
 895         {
 896           hist = 6;
 897         }
 898
 899       }
 900       else
 901         hist = 7;
 902
 903       fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]);
 904     }
 905   }
 906
 907   if (primCount)
 908   {
 909     delete[] primCount;
 910     primCount = 0;
 911   }
 912
 913   if (processed < inputCount)
 914     Printf("Only %d out of %d track(let)s were used", processed, inputCount);
 915
 916   if (eventTriggered && eventVertex)
 917   {
 918     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
 919     fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
 920   }
 921
 922    // stuff regarding the vertex reco correction and trigger bias correction
 923   fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 924
 925   if (fOption.Contains("process-types"))
 926   {
 927     // non diffractive
 928     if (processType == AliPWG0Helper::kND )
 929       fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 930
 931     // single diffractive
 932     if (processType == AliPWG0Helper::kSD)
 933       fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 934
 935     // double diffractive
 936     if (processType == AliPWG0Helper::kDD)
 937       fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 938   }
 939
 940   if (fOption.Contains("particle-species"))
 941     for (Int_t id=0; id<4; id++)
 942       fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
 943
 944   if (etaArr)
 945     delete[] etaArr;
 946   if (labelArr)
 947     delete[] labelArr;
 948   if (labelArr2)
 949     delete[] labelArr2;
 950   if (thirdDimArr)
 951     delete[] thirdDimArr;
 952   if (deltaPhiArr)
 953     delete[] deltaPhiArr;
 954 }
 955
 956 void AlidNdEtaCorrectionTask::Terminate(Option_t *)
 957 {
 958   // The Terminate() function is the last function to be called during
 959   // a query. It always runs on the client, it can be used to present
 960   // the results graphically or save the results to file.
 961
 962   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 963   if (!fOutput) {
 964     Printf("ERROR: fOutput not available");
 965     return;
 966   }
 967
 968   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
 969   fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
 970   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
 971   if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
 972   {
 973     AliDebug(AliLog::kError, "Could not read object from output list");
 974     return;
 975   }
 976
 977   fdNdEtaCorrection->Finish();
 978
 979   TString fileName;
 980   fileName.Form("correction_map%s.root", fOption.Data());
 981   TFile* fout = new TFile(fileName, "RECREATE");
 982
 983   fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
 984   if (fEsdTrackCuts)
 985     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 986
 987   fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim"));
 988   if (fEsdTrackCutsPrim)
 989     fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries");
 990
 991   fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec"));
 992   if (fEsdTrackCutsSec)
 993     fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries");
 994
 995   fdNdEtaCorrection->SaveHistograms();
 996   fdNdEtaAnalysisMC->SaveHistograms();
 997   fdNdEtaAnalysisESD->SaveHistograms();
 998
 999   if (fOutput->FindObject("dndeta_correction_ND"))
1000   {
1001     fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND"));
1002     fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD"));
1003     fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD"));
1004   }
1005   else
1006   {
1007     fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_pi"));
1008     fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_K"));
1009     fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_p"));
1010     fdNdEtaCorrectionSpecial[3] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_other"));
1011   }
1012   for (Int_t i=0; i<4; ++i)
1013     if (fdNdEtaCorrectionSpecial[i])
1014       fdNdEtaCorrectionSpecial[i]->SaveHistograms();
1015
1016   fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1017   if (fTemp1)
1018     fTemp1->Write();
1019   fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2"));
1020   if (fTemp2)
1021     fTemp2->Write();
1022
1023   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
1024   if (fVertexCorrelation)
1025     fVertexCorrelation->Write();
1026   fVertexCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelationShift"));
1027   if (fVertexCorrelationShift)
1028     fVertexCorrelationShift->Write();
1029   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
1030   if (fVertexProfile)
1031     fVertexProfile->Write();
1032   fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
1033   if (fVertexShift)
1034     fVertexShift->Write();
1035   fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
1036   if (fVertexShiftNorm)
1037     fVertexShiftNorm->Write();
1038
1039   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
1040   if (fEtaCorrelation)
1041     fEtaCorrelation->Write();
1042   fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
1043   if (fEtaCorrelationShift)
1044     fEtaCorrelationShift->Write();
1045   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
1046   if (fEtaProfile)
1047     fEtaProfile->Write();
1048   fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
1049   if (fEtaResolution)
1050     fEtaResolution->Write();
1051   fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
1052   if (fpTResolution)
1053   {
1054     fpTResolution->FitSlicesY();
1055     fpTResolution->Write();
1056   }
1057
1058   fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
1059   if (fMultAll)
1060     fMultAll->Write();
1061
1062   fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
1063   if (fMultTr)
1064     fMultTr->Write();
1065
1066   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1067   if (fMultVtx)
1068     fMultVtx->Write();
1069
1070   for (Int_t i=0; i<8; ++i)
1071   {
1072     fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
1073     if (fDeltaPhi[i])
1074       fDeltaPhi[i]->Write();
1075   }
1076
1077   fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
1078   if (fEventStats)
1079     fEventStats->Write();
1080
1081   fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles"));
1082   if (fPIDParticles)
1083     fPIDParticles->Write();
1084
1085   fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks"));
1086   if (fPIDTracks)
1087     fPIDTracks->Write();
1088
1089  //fdNdEtaCorrection->DrawHistograms();
1090
1091   Printf("Writting result to %s", fileName.Data());
1092
1093   if (fPIDParticles && fPIDTracks)
1094   {
1095     TDatabasePDG* pdgDB = new TDatabasePDG;
1096
1097     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
1098       if (fPIDParticles->GetBinContent(i) > 0)
1099       {
1100         TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i));
1101         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), (pdgParticle) ? pdgParticle->GetName() : "not found", (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
1102       }
1103
1104     delete pdgDB;
1105     pdgDB = 0;
1106   }
1107
1108   fout->Write();
1109   fout->Close();
1110 }
1111
1112 Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
1113 {
1114   // returns if a particle with this sign should be counted
1115   // this is determined by the value of fSignMode, which should have the same sign
1116   // as the charge
1117   // if fSignMode is 0 all particles are counted
1118
1119   if (fSignMode == 0)
1120     return kTRUE;
1121
1122   if (!particle)
1123   {
1124     printf("WARNING: not counting a particle that does not have a pdg particle\n");
1125     return kFALSE;
1126   }
1127
1128   Double_t charge = particle->Charge();
1129
1130   if (fSignMode > 0)
1131     if (charge < 0)
1132       return kFALSE;
1133
1134   if (fSignMode < 0)
1135     if (charge > 0)
1136       return kFALSE;
1137
1138   return kTRUE;
1139 }
1140