PWGCF/Correlations/Base/AliAnalyseLeadingTrackUE.cxx

   1 /*************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: A.Abrahantes, E.Lopez, S.Vallero                               *
   5  * Version 1.0                                                            *
   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 //#include <TBranch.h>
  16 //#include <TCanvas.h>
  17 //#include <TChain.h>
  18 //#include <TFile.h>
  19 //#include <TH1F.h>
  20 //#include <TH1I.h>
  21 //#include <TH2F.h>
  22 #include <TList.h>
  23 //#include <TLorentzVector.h>
  24 #include <TMath.h>
  25 #include <TObjArray.h>
  26 #include <TObject.h>
  27 //#include <TProfile.h>
  28 //#include <TRandom.h>
  29 //#include <TSystem.h>
  30 //#include <TTree.h>
  31 #include <TVector3.h>
  32
  33 #include "AliAnalyseLeadingTrackUE.h"
  34 //#include "AliAnalysisTask.h"
  35
  36 //#include "AliAnalysisHelperJetTasks.h"
  37 //#include "AliAnalysisManager.h"
  38 #include "AliAODEvent.h"
  39 //#include "AliAODHandler.h"
  40 //#include "AliAODJet.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliAODTrack.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDtrackCuts.h"
  46 //#include "AliGenPythiaEventHeader.h"
  47 #include "AliInputEventHandler.h"
  48 //#include "AliKFVertex.h"
  49 //#include "AliLog.h"
  50 #include "AliMCEvent.h"
  51 #include "AliVParticle.h"
  52 #include "AliAODMCHeader.h"
  53
  54 #include "AliAnalysisManager.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliStack.h"
  57 #include "AliPIDResponse.h"
  58
  59
  60 ////////////////////////////////////////////////
  61 //---------------------------------------------
  62 // Class for transverse regions analysis
  63 //---------------------------------------------
  64 ////////////////////////////////////////////////
  65
  66
  67 using namespace std;
  68
  69 ClassImp(AliAnalyseLeadingTrackUE)
  70
  71 //-------------------------------------------------------------------
  72 AliAnalyseLeadingTrackUE::AliAnalyseLeadingTrackUE() :
  73   TObject(),
  74   fDebug(0),
  75   fFilterBit(16),
  76   fOnlyHadrons(kFALSE),
  77   fTrackEtaCut(0.8),
  78   fTrackPtMin(0),
  79   fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
  80   fEsdTrackCuts(0x0),
  81   fEsdTrackCutsExtra1(0x0),
  82   fEsdTrackCutsExtra2(0x0)
  83 {
  84   // constructor
  85 }
  86
  87
  88 //-------------------------------------------------------------------
  89 AliAnalyseLeadingTrackUE & AliAnalyseLeadingTrackUE::operator = (const AliAnalyseLeadingTrackUE & /*source*/)
  90 {
  91   // assignment operator
  92   return *this;
  93 }
  94
  95
  96 //-------------------------------------------------------------------
  97 AliAnalyseLeadingTrackUE::~AliAnalyseLeadingTrackUE()
  98 {
  99
 100   //clear memory
 101
 102 }
 103
 104
 105 //____________________________________________________________________
 106 Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
 107 {
 108
 109   // select track according to set of cuts
 110   if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
 111   if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2 && !fEsdTrackCutsExtra1->IsSelected(track) && !fEsdTrackCutsExtra2->IsSelected(track)) return kFALSE;
 112
 113   return kTRUE;
 114 }
 115
 116
 117 //____________________________________________________________________
 118 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t filterbit) {
 119
 120   // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
 121   // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
 122
 123   if (filterbit == -1)
 124     filterbit = fFilterBit;
 125
 126   if (filterbit == 128)
 127   {
 128     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 129     fEsdTrackCuts->SetMinNClustersTPC(70);
 130   }
 131   else if (filterbit == 256)
 132   {
 133     // syst study
 134     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 135     fEsdTrackCuts->SetMinNClustersTPC(80);
 136     fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
 137     fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
 138     fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
 139   }
 140   else if (filterbit == 512)
 141   {
 142     // syst study
 143     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 144     fEsdTrackCuts->SetMinNClustersTPC(60);
 145     fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
 146     fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
 147     fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
 148   }
 149   else if (filterbit == 1024)
 150   {
 151     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 152     fEsdTrackCuts->SetMinNClustersTPC(-1);
 153     fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
 154     fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 155   }
 156   else if (filterbit == 2048) // mimic hybrid tracks
 157   {
 158     // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
 159     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 160     fEsdTrackCuts->SetName("Global Hybrid tracks, loose DCA");
 161     fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
 162     fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
 163     fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 164     fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 165     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 166
 167     // Add SPD requirement
 168     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 169     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 170     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
 171
 172     fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
 173     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
 174     // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
 175   }
 176   else if (filterbit == 4096) // require TOF matching
 177   {
 178     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
 179     // FIXME: TOF REQUIREMENTS ARE IN GetParticleSpecies FOR THE MOMENT
 180   }
 181   else
 182   {
 183     fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 184     fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 185
 186     // Add SPD requirement
 187     fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
 188     fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 189
 190     // Add SDD requirement
 191     fEsdTrackCutsExtra2 = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
 192     fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
 193   }
 194 }
 195
 196 //____________________________________________________________________
 197 TObjArray*  AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
 198 {
 199
 200   // Returns an array of charged particles (or jets) ordered according to their pT.
 201
 202   Int_t nTracks = NParticles(obj);
 203
 204
 205   if( !nTracks ) return 0;
 206
 207   // Define array of AliVParticle objects
 208   TObjArray* tracks = new TObjArray(nTracks);
 209
 210   // Loop over tracks or jets
 211   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 212         AliVParticle* part = ParticleWithCuts( obj, ipart );
 213         if (!part) continue;
 214         // Accept leading-tracks in a limited pseudo-rapidity range
 215         if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
 216         tracks->AddLast( part );
 217         }
 218   // Order tracks by pT
 219   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
 220
 221   nTracks = tracks->GetEntriesFast();
 222   if( !nTracks ) return 0;
 223
 224   return tracks;
 225   }
 226
 227
 228 void AliAnalyseLeadingTrackUE::RemoveInjectedSignals(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
 229 {
 230   // remove injected signals (primaries above <maxLabel>)
 231   // <tracks> can be the following cases:
 232   // a. tracks: in this case the label is taken and then case b.
 233   // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
 234   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
 235
 236   TClonesArray* arrayMC = 0;
 237   AliMCEvent* mcEvent = 0;
 238   if (mcObj->InheritsFrom("AliMCEvent"))
 239     mcEvent = static_cast<AliMCEvent*>(mcObj);
 240   else if (mcObj->InheritsFrom("TClonesArray"))
 241     arrayMC = static_cast<TClonesArray*>(mcObj);
 242   else
 243   {
 244     arrayMC->Dump();
 245     AliFatal("Invalid object passed");
 246   }
 247
 248   Int_t before = tracks->GetEntriesFast();
 249
 250   for (Int_t i=0; i<before; ++i)
 251   {
 252     AliVParticle* part = (AliVParticle*) tracks->At(i);
 253
 254     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
 255       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
 256
 257     AliVParticle* mother = part;
 258     if (mcEvent)
 259     {
 260       while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 261       {
 262         if (((AliMCParticle*)mother)->GetMother() < 0)
 263         {
 264           mother = 0;
 265           break;
 266         }
 267
 268         mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
 269         if (!mother)
 270           break;
 271       }
 272     }
 273     else
 274     {
 275       // find the primary mother
 276       while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 277       {
 278         if (((AliAODMCParticle*)mother)->GetMother() < 0)
 279         {
 280           mother = 0;
 281           break;
 282         }
 283
 284         mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 285         if (!mother)
 286           break;
 287       }
 288     }
 289
 290     if (!mother)
 291     {
 292       Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 293       continue;
 294     }
 295
 296     if (mother->GetLabel() > maxLabel)
 297     {
 298 //       Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
 299       TObject* object = tracks->RemoveAt(i);
 300       if (tracks->IsOwner())
 301         delete object;
 302     }
 303   }
 304
 305   tracks->Compress();
 306
 307   AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 308 }
 309
 310 //-------------------------------------------------------------------
 311 void AliAnalyseLeadingTrackUE::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
 312 {
 313   // remove particles from weak decays
 314   // <tracks> can be the following cases:
 315   // a. tracks: in this case the label is taken and then case b.
 316   // b. particles: it is checked if IsSecondaryFromWeakDecay is true
 317   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
 318
 319   TClonesArray* arrayMC = 0;
 320   AliMCEvent* mcEvent = 0;
 321   if (mcObj->InheritsFrom("AliMCEvent"))
 322     mcEvent = static_cast<AliMCEvent*>(mcObj);
 323   else if (mcObj->InheritsFrom("TClonesArray"))
 324     arrayMC = static_cast<TClonesArray*>(mcObj);
 325   else
 326   {
 327     arrayMC->Dump();
 328     AliFatal("Invalid object passed");
 329   }
 330
 331   Int_t before = tracks->GetEntriesFast();
 332
 333   for (Int_t i=0; i<before; ++i)
 334   {
 335     AliVParticle* part = (AliVParticle*) tracks->At(i);
 336
 337     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
 338       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
 339
 340     if (part->InheritsFrom("AliAODMCParticle"))
 341     {
 342       if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
 343         continue;
 344     }
 345     else if (part->InheritsFrom("AliMCParticle") && mcEvent)
 346     {
 347       if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
 348         continue;
 349     }
 350     else
 351     {
 352       part->Dump();
 353       AliFatal("Unknown particle");
 354     }
 355
 356 //     Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
 357     TObject* object = tracks->RemoveAt(i);
 358     if (tracks->IsOwner())
 359       delete object;
 360   }
 361
 362   tracks->Compress();
 363
 364   if (before > tracks->GetEntriesFast())
 365     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
 366 }
 367
 368 //-------------------------------------------------------------------
 369 TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
 370 {
 371   // Returns an array of particles that pass the cuts, if arrayMC is given each reconstructed particle is replaced by its corresponding MC particles, depending on the parameter onlyprimaries only for primaries
 372   // particleSpecies: -1 all particles are returned
 373   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 374
 375   Int_t nTracks = NParticles(obj);
 376   TObjArray* tracks = new TObjArray;
 377
 378   // for TPC only tracks
 379   Bool_t hasOwnership = kFALSE;
 380   if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
 381     hasOwnership = kTRUE;
 382
 383   if (!arrayMC)
 384     tracks->SetOwner(hasOwnership);
 385
 386   // Loop over tracks or jets
 387   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 388     AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
 389     if (!part) continue;
 390
 391     if (useEtaPtCuts)
 392       if (TMath::Abs(part->Eta()) > fTrackEtaCut || part->Pt() < fTrackPtMin)
 393       {
 394         if (hasOwnership)
 395           delete part;
 396         continue;
 397       }
 398
 399     if (arrayMC) {
 400       Int_t label = part->GetLabel();
 401       if (hasOwnership)
 402         delete part;
 403       // re-define part as the matched MC particle
 404       part = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
 405       if (!part)continue;
 406     }
 407
 408     tracks->Add(part);
 409   }
 410
 411   return tracks;
 412 }
 413
 414 //-------------------------------------------------------------------
 415 TObjArray* AliAnalyseLeadingTrackUE::GetFakeParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
 416 {
 417   // particleSpecies: -1 all particles are returned
 418   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 419
 420   Int_t nTracks = NParticles(obj);
 421   TObjArray* tracksReconstructed = new TObjArray;
 422   TObjArray* tracksOriginal = new TObjArray;
 423   TObjArray* tracksFake = new TObjArray;
 424
 425   // for TPC only tracks
 426   Bool_t hasOwnership = kFALSE;
 427   if ((fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024 || fFilterBit == 2048) && obj->InheritsFrom("AliESDEvent"))
 428     hasOwnership = kTRUE;
 429
 430   tracksReconstructed->SetOwner(hasOwnership);
 431   tracksFake->SetOwner(hasOwnership);
 432
 433   // Loop over tracks or jets
 434   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 435     AliVParticle* partReconstructed = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
 436     if (!partReconstructed) continue;
 437
 438     if (useEtaPtCuts)
 439       if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || partReconstructed->Pt() < fTrackPtMin)
 440       {
 441         if (hasOwnership)
 442           delete partReconstructed;
 443         continue;
 444       }
 445
 446     Int_t label = partReconstructed->GetLabel();
 447     if (label == 0)
 448     {
 449       /*
 450       Printf(">>> TPC only track:");
 451       partReconstructed->Print();
 452       partReconstructed->Dump();
 453       Printf(">>> Global track:");
 454       ((AliESDEvent*) obj)->GetTrack(ipart)->Dump();
 455       Printf("Fake (TPC only): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", partReconstructed->Eta(), partReconstructed->Phi(), partReconstructed->Pt(), ((AliESDtrack*) partReconstructed)->GetTPCclusters(0), ((AliESDtrack*) partReconstructed)->GetTPCsignal());
 456       Printf("Fake (global  ): eta = %f, phi = %f, pT = %f, ncl = %d, dedx = %f", ((AliESDEvent*) obj)->GetTrack(ipart)->Eta(), ((AliESDEvent*) obj)->GetTrack(ipart)->Phi(), ((AliESDEvent*) obj)->GetTrack(ipart)->Pt(), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCclusters(0), ((AliESDEvent*) obj)->GetTrack(ipart)->GetTPCsignal());
 457       */
 458       tracksFake->AddLast(partReconstructed);
 459       continue;
 460     }
 461
 462     AliVParticle* partOriginal = ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries, particleSpecies);
 463     if (!partOriginal)
 464     {
 465       if (hasOwnership)
 466         delete partReconstructed;
 467       continue;
 468     }
 469
 470     tracksReconstructed->AddLast(partReconstructed);
 471     tracksOriginal->AddLast(partOriginal);
 472   }
 473   TObjArray* pairs = new TObjArray;
 474   pairs->SetOwner(kTRUE);
 475   pairs->Add(tracksReconstructed);
 476   pairs->Add(tracksOriginal);
 477   pairs->Add(tracksFake);
 478   return pairs;
 479 }
 480
 481 //-------------------------------------------------------------------
 482 TObjArray* AliAnalyseLeadingTrackUE::GetMinMaxRegion(TList *transv1, TList *transv2)
 483 {
 484
 485  // Returns two lists of particles, one for MIN and one for MAX region
 486   Double_t sumpT1 = 0.;
 487   Double_t sumpT2 = 0.;
 488
 489   Int_t particles1 = transv1->GetEntries();
 490   Int_t particles2 = transv2->GetEntries();
 491
 492   // Loop on transverse region 1
 493   for (Int_t i=0; i<particles1; i++){
 494         AliVParticle *part = (AliVParticle*)transv1->At(i);
 495         sumpT1 +=  part->Pt();
 496         }
 497
 498   // Loop on transverse region 2
 499   for (Int_t i=0; i<particles2; i++){
 500         AliVParticle *part = (AliVParticle*)transv2->At(i);
 501         sumpT2 +=  part->Pt();
 502         }
 503
 504   TObjArray *regionParticles = new TObjArray;
 505   if ( sumpT2 >= sumpT1 ){
 506         regionParticles->AddLast(transv1); // MIN
 507         regionParticles->AddLast(transv2); // MAX
 508   }else {
 509         regionParticles->AddLast(transv2); // MIN
 510         regionParticles->AddLast(transv1); // MAX
 511         }
 512
 513  return regionParticles;
 514 }
 515
 516 //-------------------------------------------------------------------
 517 Int_t  AliAnalyseLeadingTrackUE::NParticles(TObject* obj)
 518 {
 519
 520   //Returns the number of particles in AliAODMCParticle array  or AliAODTracks or AliESDTracks
 521
 522   Int_t nTracks;
 523
 524   if (obj->InheritsFrom("TClonesArray")){ // MC particles
 525         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 526         nTracks = arrayMC->GetEntriesFast();
 527   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 528         TObjArray *array = static_cast<TObjArray*>(obj);
 529         nTracks = array->GetEntriesFast();
 530   }else if (obj->InheritsFrom("AliAODEvent")){  // RECO AOD tracks
 531         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 532         nTracks = aodEvent->GetNTracks();
 533   }else if (obj->InheritsFrom("AliESDEvent")){  // RECO ESD tracks
 534         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 535         nTracks = esdEvent->GetNumberOfTracks();
 536   }else if (obj->InheritsFrom("AliMCEvent")){  // RECO ESD tracks
 537         AliMCEvent *mcEvent = static_cast<AliMCEvent*>(obj);
 538         nTracks = mcEvent->GetNumberOfTracks();
 539   }else {
 540         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 541         return 0;
 542         }
 543
 544   return nTracks;
 545 }
 546
 547 //-------------------------------------------------------------------
 548 AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ipart, Bool_t onlyprimaries, Int_t particleSpecies)
 549 {
 550   // Returns track or MC particle at position "ipart" if passes selection criteria
 551   // particleSpecies: -1 all particles are returned
 552   //                  0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
 553   AliVParticle *part=0;
 554
 555   if (obj->InheritsFrom("TClonesArray")){ // AOD-MC PARTICLE
 556         TClonesArray *arrayMC = static_cast<TClonesArray*>(obj);
 557         part = (AliVParticle*)arrayMC->At( ipart );
 558         if (!part)return 0;
 559         // eventually only primaries
 560         if (onlyprimaries && !( ((AliAODMCParticle*)part)->IsPhysicalPrimary()) )return 0;
 561         // eventually only hadrons
 562         if (fOnlyHadrons){
 563                 Int_t pdgCode = ((AliAODMCParticle*)part)->GetPdgCode();
 564                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 565                                   TMath::Abs(pdgCode)==2212 || // Proton
 566                                   TMath::Abs(pdgCode)==321;    // Kaon
 567                 if (!isHadron) return 0;
 568                 }
 569         if (particleSpecies != -1) {
 570                 // find the primary mother
 571                 AliVParticle* mother = part;
 572                 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 573                 {
 574                   if (((AliAODMCParticle*)mother)->GetMother() < 0)
 575                   {
 576                     mother = 0;
 577                     break;
 578                   }
 579
 580                   mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 581                   if (!mother)
 582                     break;
 583                 }
 584
 585                 if (mother)
 586                 {
 587                   Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
 588                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 589                           return 0;
 590                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 591                           return 0;
 592                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 593                           return 0;
 594                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 595                           return 0;
 596                 }
 597                 else
 598                 {
 599                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 600                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 601                   part->Print();
 602
 603                   /*
 604                   // this code prints the details of the mother that is missing in the AOD
 605                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 606
 607                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 608
 609                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 610                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 611                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 612                   */
 613
 614                   if (particleSpecies != 3)
 615                     return 0;
 616                 }
 617         }
 618
 619   }else if (obj->InheritsFrom("TObjArray")){ // list of AliVParticle
 620         TObjArray *array = static_cast<TObjArray*>(obj);
 621         part = (AliVParticle*)array->At( ipart );
 622         if (!part)return 0;
 623   }else if (obj->InheritsFrom("AliMCEvent")){ // MC PARTICLE
 624         AliMCEvent* mcEvent =  static_cast<AliMCEvent*>(obj);
 625         part = mcEvent->GetTrack( ipart );
 626         if (!part) return 0;
 627         // eventually only primaries
 628         if (onlyprimaries && !( mcEvent->IsPhysicalPrimary(ipart)) )return 0;
 629         // eventually only hadrons
 630         //TO-DO
 631         /*if (fOnlyHadrons){
 632                 Int_t pdgCode =  part->GetPdgCode();
 633                 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||  // Pion
 634                                   TMath::Abs(pdgCode)==2212 || // Proton
 635                                   TMath::Abs(pdgCode)==321;    // Kaon
 636                 if (!isHadron) return 0;
 637                 }
 638         */
 639        if (particleSpecies != -1) {
 640                 // find the primary mother
 641                 AliMCParticle* mother = (AliMCParticle*) part;
 642 //              Printf("");
 643                 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
 644                 {
 645 //                Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
 646                   if (mother->GetMother() < 0)
 647                   {
 648                     mother = 0;
 649                     break;
 650                   }
 651
 652                   mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
 653                   if (!mother)
 654                     break;
 655                 }
 656
 657                 if (mother)
 658                 {
 659                   Int_t pdgCode = mother->PdgCode();
 660                   if (particleSpecies == 0 && TMath::Abs(pdgCode)!=211)
 661                           return 0;
 662                   if (particleSpecies == 1 && TMath::Abs(pdgCode)!=321)
 663                           return 0;
 664                   if (particleSpecies == 2 && TMath::Abs(pdgCode)!=2212)
 665                           return 0;
 666                   if (particleSpecies == 3 && (TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212))
 667                           return 0;
 668                 }
 669                 else
 670                 {
 671                   // if mother not found, accept particle only in case of particleSpecies == 3. To include it in all or no sample is no solution
 672                   Printf("WARNING: No mother found for particle %d:", part->GetLabel());
 673                   //part->Dump();
 674                   //part->Print();
 675
 676                   /*
 677                   // this code prints the details of the mother that is missing in the AOD
 678                   AliMCEventHandler* fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 679
 680                   AliMCEvent* fMcEvent = fMcHandler->MCEvent();
 681
 682                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->Print();
 683                   fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Print();
 684                   Printf("eta = %f", fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(fMcEvent->Stack()->Particle(part->GetLabel())->GetMother(0))->GetMother(0))->Eta());
 685                   */
 686
 687                   if (particleSpecies != 3)
 688                     return 0;
 689                 }
 690         }
 691   }else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
 692         AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
 693         part = aodEvent->GetTrack(ipart);
 694         // track selection cuts
 695         if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
 696         //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
 697         // only primary candidates
 698         //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
 699         // eventually only hadrons
 700         if (fOnlyHadrons){
 701                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 702                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 703                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 704                 if (!isHadron) return 0;
 705                 }
 706
 707   }else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
 708         AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
 709         part = esdEvent->GetTrack(ipart);
 710         if (!part)return 0;
 711         // track selection cuts
 712
 713         if (!( ApplyCuts(part)) )
 714          return 0;
 715
 716         if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
 717         {
 718           // create TPC only tracks constrained to the SPD vertex
 719
 720           const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
 721
 722           AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
 723           if(!track) return 0;
 724
 725 //        Printf(">%f %f %f", track->Eta(), track->Phi(), track->Pt());
 726
 727           if(track->Pt()>0.){
 728             // only constrain tracks above threshold
 729             AliExternalTrackParam exParam;
 730             // take the B-feild from the ESD, no 3D fieldMap available at this point
 731             Bool_t relate = kFALSE;
 732             relate = track->RelateToVertexTPC(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
 733             if(!relate)
 734             {
 735 //                 Printf("relating failed");
 736               delete track;
 737               return 0;
 738             }
 739             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 740           }
 741
 742 //        Printf(">%f %f %f\n", track->Eta(), track->Phi(), track->Pt());
 743
 744           part = track;
 745         }
 746         else if (fFilterBit == 2048)
 747         {
 748           // hybrid tracks
 749
 750           // clone
 751           AliESDtrack* esdTrack = new AliESDtrack(*((AliESDtrack*) part));
 752 //        Printf("%d %d %d %d %d", fEsdTrackCuts->IsSelected(esdTrack), fEsdTrackCutsExtra1->IsSelected(esdTrack), fEsdTrackCutsExtra2->IsSelected(esdTrack), esdTrack->HasPointOnITSLayer(0), esdTrack->HasPointOnITSLayer(1));
 753
 754           if (fEsdTrackCutsExtra2->IsSelected(esdTrack))
 755           {
 756 //          Float_t ptBefore = esdTrack->Pt();
 757             // set constrained pT as default pT
 758             if (!esdTrack->GetConstrainedParam())
 759               return 0;
 760             esdTrack->CopyFromVTrack(esdTrack->GetConstrainedParam());
 761 //          Printf("%f %f", ptBefore, esdTrack->Pt());
 762           }
 763           part = esdTrack;
 764         }
 765
 766         // eventually only hadrons
 767         //TO-DO
 768         /*if (fOnlyHadrons){
 769                 Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
 770                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kKaon ||
 771                                   ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kProton;
 772                 if (!isHadron) return 0;
 773                 }
 774         */
 775         if (particleSpecies != -1 && GetParticleSpecies((AliVTrack*) part) != particleSpecies) return 0; // If it is -1 you take all the particles
 776
 777   }else {
 778         if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
 779         return 0;
 780         }
 781
 782   // only charged
 783   if (!part->Charge())return 0;
 784
 785   return part;
 786 }
 787
 788
 789 //-------------------------------------------------------------------
 790 void  AliAnalyseLeadingTrackUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 791 {
 792   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 793
 794   static TObject *tmp;
 795   static int i;           // "static" to save stack space
 796   int j;
 797
 798   while (last - first > 1) {
 799     i = first;
 800     j = last;
 801     for (;;) {
 802       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 803         ;
 804       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 805         ;
 806       if (i >= j)
 807         break;
 808
 809       tmp  = a[i];
 810       a[i] = a[j];
 811       a[j] = tmp;
 812     }
 813     if (j == first) {
 814       ++first;
 815       continue;
 816     }
 817     tmp = a[first];
 818     a[first] = a[j];
 819     a[j] = tmp;
 820     if (j - first < last - (j + 1)) {
 821       QSortTracks(a, first, j);
 822       first = j + 1;   // QSortTracks(j + 1, last);
 823     } else {
 824       QSortTracks(a, j + 1, last);
 825       last = j;        // QSortTracks(first, j);
 826     }
 827   }
 828 }
 829
 830 //____________________________________________________________________
 831 TObjArray*  AliAnalyseLeadingTrackUE::SortRegions(const AliVParticle* leading, TObject* obj, TObject* arrayMC, Bool_t onlyprimaries)
 832 {
 833
 834   // Assign particles to towards, away or transverse regions.
 835   // Returns a lists of particles for each region.
 836
 837   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 838   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 839
 840   // Define output lists of particles
 841   TList *toward = new TList();
 842   TList *away = new TList();
 843   // Two transverse regions, for the moment those are not yet MIN and MAX!!!
 844   // MIN and MAX can be sorted in GetMinMaxRegion function
 845   TList *transverse1 = new TList();
 846   TList *transverse2 = new TList();
 847
 848   TObjArray *regionParticles = new TObjArray;
 849   regionParticles->SetOwner(kTRUE);
 850
 851   regionParticles->AddLast(toward);
 852   regionParticles->AddLast(away);
 853   regionParticles->AddLast(transverse1);
 854   regionParticles->AddLast(transverse2);
 855
 856   if (!leading)
 857     return regionParticles;
 858
 859   // Switch to vector for leading particle
 860   TVector3 leadVect(leading->Px(),leading->Py(),leading->Pz());
 861
 862   Int_t nTracks = NParticles(obj);
 863   if( !nTracks ) return 0;
 864   // Loop over tracks
 865   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 866         AliVParticle* part = ParticleWithCuts(obj, ipart);
 867         if (!part)continue;
 868         //Switch to vectors for particles
 869         TVector3 partVect(part->Px(), part->Py(), part->Pz());
 870
 871         Int_t region = 0;
 872         if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
 873         // transverse regions
 874         if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
 875         if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1;   //right
 876
 877         if (TMath::Abs(leadVect.DeltaPhi(partVect)) < k60rad ) region = 2;    //forward
 878         if (TMath::Abs(leadVect.DeltaPhi(partVect)) > k120rad ) region = -2;  //backward
 879
 880         // skip leading particle
 881         if (leading == part)
 882           continue;
 883
 884         if (!region)continue;
 885         if (arrayMC && arrayMC->InheritsFrom("TClonesArray") && obj->InheritsFrom("AliAODEvent")){
 886                 Int_t label = ((AliAODTrack*)part)->GetLabel();
 887                 // re-define part as the matched MC particle
 888                 part = (AliAODMCParticle*)ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries);
 889                 if (!part)continue;
 890                 // skip leading particle
 891                 if (leading == part)
 892                   continue;
 893                 }
 894         if (arrayMC && arrayMC->InheritsFrom("AliMCEvent") && obj->InheritsFrom("AliESDEvent")){
 895                 Int_t label = ((AliESDtrack*)part)->GetLabel();
 896                 // look for the matched MC particle (but do not re-define part)
 897                 if (!ParticleWithCuts(arrayMC, TMath::Abs(label),onlyprimaries)) continue;
 898                 }
 899
 900         if ( region == 1 ) transverse1->Add(part);
 901         if ( region == -1 ) transverse2->Add(part);
 902         if ( region == 2 ) toward->Add(part);
 903         if ( region == -2 ) away->Add(part);
 904
 905         }//end loop on tracks
 906
 907   return regionParticles;
 908
 909 }
 910
 911
 912 //____________________________________________________________________
 913 Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
 914 {
 915   if (!obj) // MC
 916     return kFALSE;
 917
 918   // Use AliPhysicsSelection to select good events, works for ESD and AOD
 919   if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(fEventSelection)))
 920     return kFALSE;
 921
 922   return kTRUE;
 923 }
 924
 925 //____________________________________________________________________
 926 Bool_t  AliAnalyseLeadingTrackUE::VertexSelection(const TObject* obj, Int_t ntracks, Double_t zed)
 927 {
 928
 929   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 930
 931   if (obj->InheritsFrom("AliAODEvent")){
 932         Int_t nVertex = ((AliAODEvent*)obj)->GetNumberOfVertices();
 933         if( nVertex > 0 ) {
 934                 AliAODVertex* vertex = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
 935                 Int_t nTracksPrim = vertex->GetNContributors();
 936                 Double_t zVertex = vertex->GetZ();
 937                 if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 938                 // Reject TPC only vertex
 939                 TString name(vertex->GetName());
 940                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 941
 942                 // Select a quality vertex by number of tracks?
 943                 if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 944                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 945                         return kFALSE;
 946                         }
 947                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 948                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 949                 //  return kFALSE;
 950                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 951                 } else {
 952                         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 953                         return kFALSE;
 954                         }
 955         }
 956
 957   if (obj->InheritsFrom("AliMCEvent"))
 958   {
 959     if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
 960     {
 961       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 962       return kFALSE;
 963     }
 964   }
 965
 966   if (obj->InheritsFrom("AliAODMCHeader"))
 967   {
 968     if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
 969     {
 970       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
 971       return kFALSE;
 972     }
 973   }
 974
 975   // ESD case for DCA studies
 976   if (obj->InheritsFrom("AliESDEvent")){
 977        AliESDVertex* vertex = (AliESDVertex*)((AliESDEvent*)obj)->GetPrimaryVertex();
 978        if ( vertex){
 979                Int_t nTracksPrim = vertex->GetNContributors();
 980                Double_t zVertex = vertex->GetZ();
 981                if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 982                // Reject SPD or TPC only vertex
 983                TString name(vertex->GetName());
 984                if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
 985
 986                // Select a quality vertex by number of tracks?
 987                if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
 988                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 989                        return kFALSE;
 990                        }
 991                // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
 992                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
 993                 //  return kFALSE;
 994                if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 995                } else {
 996                        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 997                        return kFALSE;
 998                        }
 999        }
1000
1001   return kTRUE;
1002 }
1003
1004
1005 Int_t AliAnalyseLeadingTrackUE::GetParticleSpecies(AliVTrack      * trk)
1006 {
1007   // return PID according to detectors
1008   // Get PID response object, if needed
1009   // FIXME: THIS IS UGLY!! THE CODE IS COPIED AND ADAPTED FROM PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothPID.cxx
1010   // TO BE REPLACED BY THE PROPER CLASS ASAP
1011   // MF 10/12/12
1012
1013
1014   enum BothPIDType_t {kNSigmaTPC,kNSigmaTOF,kNSigmaTPCTOF};
1015
1016   enum BothParticleSpecies_t
1017   {
1018     kSpPion = 0,     kSpKaon,     kSpProton,
1019     kNSpecies,
1020     kSpUndefined,
1021   }; // Particle species used in plotting
1022
1023   UInt_t fPIDType = kNSigmaTPCTOF;
1024
1025   // Check TOF matching & PID
1026   // FIXME: this should go in the track selection!!
1027   UInt_t status;
1028   status=trk->GetStatus();
1029   if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0){//kTOFout and kTIME
1030     return kSpUndefined;
1031   }
1032
1033
1034
1035   // guess the particle based on the smaller nsigma
1036
1037   // rec[kSpPion]=false;
1038   // rec[kSpKaon]=false;
1039   // rec[kSpProton]=false;
1040
1041   static AliPIDResponse * fPIDResponse = 0;
1042
1043   if(!fPIDResponse)
1044     {
1045       AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
1046       AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
1047       fPIDResponse = inputHandler->GetPIDResponse();
1048     }
1049
1050   if(!fPIDResponse)
1051     {
1052       AliFatal("Cannot get pid response");
1053       return 0;
1054     }
1055
1056
1057   // Compute nsigma for each hypthesis
1058   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
1059
1060   const Double_t fshiftTPC = 0, fshiftTOF=0;
1061   const Double_t fNSigmaPID = 3; // FIXME
1062   // --- TPC
1063   Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
1064   Double_t nsigmaTPCkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
1065   Double_t nsigmaTPCkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
1066
1067
1068   Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
1069
1070   Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
1071   Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
1072   Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
1073   if(fPIDType==kNSigmaTOF)
1074     {
1075       nsigmaTPCTOFkProton = 100.0;
1076       nsigmaTPCTOFkKaon   = 100.0;
1077       nsigmaTPCTOFkPion   = 100.0;
1078
1079     }
1080
1081
1082   nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
1083   nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
1084   nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
1085   // the TOF info is used in combined
1086   nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
1087   nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
1088   nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
1089
1090   // select the nsigma to be used for the actual PID
1091   Double_t nsigmaPion=100;
1092   Double_t nsigmaKaon=100;
1093   Double_t nsigmaProton=100;
1094
1095   switch (fPIDType)
1096     {
1097     case kNSigmaTPC:
1098       nsigmaProton  =  nsigmaTPCkProton;
1099       nsigmaKaon    =  nsigmaTPCkKaon  ;
1100       nsigmaPion    =  nsigmaTPCkPion  ;
1101       break;
1102     case kNSigmaTOF:
1103       nsigmaProton  =  nsigmaTOFkProton;
1104       nsigmaKaon    =  nsigmaTOFkKaon  ;
1105       nsigmaPion    =  nsigmaTOFkPion  ;
1106       break;
1107     case kNSigmaTPCTOF:
1108       nsigmaProton  =  nsigmaTPCTOFkProton;
1109       nsigmaKaon    =  nsigmaTPCTOFkKaon  ;
1110       nsigmaPion    =  nsigmaTPCTOFkPion  ;
1111       break;
1112     }
1113
1114   // if(nsigmaPion   < fNSigmaPID)
1115   //   rec[kSpPion]=true;
1116   // if(nsigmaKaon   < fNSigmaPID)
1117   //   rec[kSpKaon]=true;
1118   // if(nsigmaProton   < fNSigmaPID)
1119   //   rec[kSpProton]=true;
1120
1121   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
1122     return kSpUndefined;//if is the default value for the three
1123   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID))
1124     {
1125       return kSpKaon;
1126     }
1127   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID))
1128     {
1129       return kSpPion;
1130     }
1131   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
1132     {
1133       return kSpProton;
1134     }
1135   // else, return undefined
1136   return kSpUndefined;
1137
1138 }