PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODPID.cxx

   1 #include "AliSpectraAODPID.h"
   2 #include "AliAODEvent.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TList.h"
   6 #include "AliAODTrack.h"
   7 #include "AliAODMCParticle.h"
   8 #include "AliPIDResponse.h"
   9 #include "AliAnalysisManager.h"
  10 #include "AliInputEventHandler.h"
  11 #include "AliSpectraAODTrackCuts.h"
  12
  13 ClassImp(AliSpectraAODPID)
  14
  15 AliSpectraAODPID::AliSpectraAODPID() : TNamed("PID", "PID object"), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0) {
  16
  17 }
  18
  19 AliSpectraAODPID::AliSpectraAODPID(AODPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0) {
  20
  21
  22
  23 }
  24
  25
  26
  27 void AliSpectraAODPID::FillQAHistos(AliSpectraAODHistoManager * hman, AliAODTrack * track, AliSpectraAODTrackCuts * trackCuts) {
  28
  29   // fill a bunch of QA histos
  30
  31
  32   // Get PID response object, if needed
  33   if(!fPIDResponse) {
  34     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
  35     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
  36     fPIDResponse = inputHandler->GetPIDResponse();
  37   }
  38
  39   //Response
  40   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);
  41
  42   hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
  43   hman->GetPIDHistogram(kHistPIDTPCPion)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kPion)*track->Charge()); // Expected PIDPion histo
  44   hman->GetPIDHistogram(kHistPIDTPCKaon)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kKaon)*track->Charge()); // Expected PIDKaon histo
  45   hman->GetPIDHistogram(kHistPIDTPCProton)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kProton)*track->Charge()); // Expected PIDProton histo
  46
  47   Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
  48   Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
  49   Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
  50   Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;
  51
  52   if(track->Pt()>trackCuts->GetPtTOFMatching()){
  53     hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
  54
  55     nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
  56     nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
  57     nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
  58
  59     //TOF
  60     hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
  61     hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
  62     hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
  63     hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
  64     hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
  65     hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
  66   }
  67
  68   Double_t nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2);
  69   Double_t nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2);
  70   Double_t nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2);
  71
  72   //TPC
  73   hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
  74   hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
  75   hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
  76   hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
  77   hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
  78   hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
  79   //TPCTOF
  80   hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
  81   hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
  82   hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
  83   hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
  84   hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
  85   hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
  86
  87 }
  88
  89 Int_t AliSpectraAODPID::GetParticleSpecie(AliAODMCParticle * part) {
  90   // return PID according to MC truth
  91   switch(TMath::Abs(part->PdgCode())){
  92   case 2212:
  93     return kSpProton;
  94     break;
  95   case 321:
  96     return kSpKaon;
  97     break;
  98   case 211:
  99     return kSpPion;
 100     break;
 101   default:
 102     return kSpUndefined;
 103   }
 104 }
 105
 106
 107 Int_t AliSpectraAODPID::GetParticleSpecie(AliSpectraAODHistoManager * hman,AliAODTrack      * trk, AliSpectraAODTrackCuts * trackCuts) {
 108   // return PID according to detectors
 109
 110   // Get PID response object, if needed
 111   if(!fPIDResponse) {
 112     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 113     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 114     fPIDResponse = inputHandler->GetPIDResponse();
 115   }
 116
 117   if(!fPIDResponse) {
 118     AliFatal("Cannot get pid response");
 119     return 0;
 120   }
 121
 122
 123   // Compute nsigma for each hypthesis
 124   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
 125
 126   // --- TPC
 127   Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
 128   Double_t nsigmaTPCkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
 129   Double_t nsigmaTPCkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
 130   // --- TOF
 131   Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;
 132   if(trk->Pt()>trackCuts->GetPtTOFMatching()){
 133     nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
 134     nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
 135     nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
 136   }
 137
 138   // --- combined
 139   Double_t nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 140   Double_t nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 141   Double_t nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 142
 143
 144   // select the nsigma to be used for the actual PID
 145   Double_t nsigmaPion=0, nsigmaKaon=0, nsigmaProton=0;
 146
 147   switch (fPIDType) {
 148   case kNSigmaTPC:
 149     nsigmaProton  =  nsigmaTPCkProton;
 150     nsigmaKaon    =  nsigmaTPCkKaon  ;
 151     nsigmaPion    =  nsigmaTPCkPion  ;
 152     break;
 153   case kNSigmaTOF:
 154     nsigmaProton  =  nsigmaTOFkProton;
 155     nsigmaKaon    =  nsigmaTOFkKaon  ;
 156     nsigmaPion    =  nsigmaTOFkPion  ;
 157     break;
 158   case kNSigmaTPCTOF:
 159     nsigmaProton  =  nsigmaTPCTOFkProton;
 160     nsigmaKaon    =  nsigmaTPCTOFkKaon  ;
 161     nsigmaPion    =  nsigmaTPCTOFkPion  ;
 162     break;
 163   }
 164
 165   // guess the particle based on the smaller nsigma
 166   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
 167   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)){
 168     hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
 169     return kSpKaon;
 170   }
 171   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)){
 172     hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
 173     return kSpPion;
 174   }
 175   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)){
 176     hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
 177     return kSpProton;
 178   }
 179   // else, return undefined
 180   return kSpUndefined;
 181
 182 }
 183
 184 Long64_t AliSpectraAODPID::Merge(TCollection* list)
 185 {
 186   // Merging interface.
 187   // Returns the number of merged objects (including this).
 188
 189   Printf("Merging");
 190
 191   if (!list)
 192     return 0;
 193
 194   if (list->IsEmpty())
 195     return 1;
 196
 197   TIterator* iter = list->MakeIterator();
 198   TObject* obj;
 199
 200   // Actually, we don't do anything here...
 201   // collections of all histograms
 202   //  TList collections;
 203
 204   Int_t count = 0;
 205
 206   while ((obj = iter->Next())) {
 207     AliSpectraAODPID* entry = dynamic_cast<AliSpectraAODPID*> (obj);
 208     if (entry == 0)
 209       continue;
 210
 211     // TH1I * histo = entry->GetHistoCuts();
 212     // collections.Add(histo);
 213     count++;
 214   }
 215
 216   //  fHistoCuts->Merge(&collections);
 217
 218   delete iter;
 219   Printf("OK");
 220   return count+1;
 221 }
 222