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