PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothPID.cxx

   1 #include "AliSpectraBothPID.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 "AliSpectraBothTrackCuts.h"
  12
  13 ClassImp(AliSpectraBothPID)
  14
  15 AliSpectraBothPID::AliSpectraBothPID() : TNamed("PID", "PID object"), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0),fshiftTPC(0),fshiftTOF(0)
  16 {
  17
  18 }
  19
  20 AliSpectraBothPID::AliSpectraBothPID(BothPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0), fshiftTPC(0),fshiftTOF(0)
  21 {
  22
  23
  24
  25 }
  26
  27
  28
  29 void AliSpectraBothPID::FillQAHistos(AliSpectraBothHistoManager * hman, AliVTrack * track, AliSpectraBothTrackCuts * trackCuts)
  30 {
  31
  32   // fill a bunch of QA histos
  33
  34
  35   // Get PID response object, if needed
  36         if(!fPIDResponse)
  37         {
  38                 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
  39                  AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
  40                 fPIDResponse = inputHandler->GetPIDResponse();
  41         }
  42
  43   //Response
  44         AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);
  45
  46         hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
  47         hman->GetPIDHistogram(kHistPIDTPCPion)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kPion)*track->Charge()); // Expected PIDPion histo
  48         hman->GetPIDHistogram(kHistPIDTPCKaon)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kKaon)*track->Charge()); // Expected PIDKaon histo
  49         hman->GetPIDHistogram(kHistPIDTPCProton)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kProton)*track->Charge()); // Expected PIDProton histo
  50
  51         Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
  52         Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
  53         Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
  54         Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
  55
  56         Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
  57         Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
  58         Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
  59         if(fPIDType==kNSigmaTOF)
  60         {
  61                 nsigmaTPCTOFkProton = 100.0;
  62                 nsigmaTPCTOFkKaon   = 100.0;
  63                 nsigmaTPCTOFkPion   =100.0;
  64
  65         }
  66
  67         if(track->Pt()>trackCuts->GetPtTOFMatching())
  68         {
  69
  70                  hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
  71
  72                  nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
  73                  nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
  74                  nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
  75
  76     //TOF
  77                  hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
  78                 hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
  79                 hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
  80                 hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
  81                 hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
  82                 hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
  83
  84                  nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2);
  85                 nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2);
  86                 nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2);
  87
  88         }
  89         else
  90         {
  91                 if (fPIDType==kNSigmaTOF)
  92                         return;
  93         }
  94
  95   //TPC
  96         hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
  97         hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
  98         hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
  99          hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
 100         hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
 101          hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
 102   //TPCTOF
 103         hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
 104         hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
 105         hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
 106         hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
 107         hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
 108         hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
 109
 110 }
 111
 112 Int_t AliSpectraBothPID::GetParticleSpecie(AliAODMCParticle * part)
 113 {
 114   // return PID according to MC truth
 115         switch(TMath::Abs(part->PdgCode()))
 116         {
 117                 case 2212:
 118                 return kSpProton;
 119                  break;
 120                 case 321:
 121                 return kSpKaon;
 122                 break;
 123                 case 211:
 124                 return kSpPion;
 125                 break;
 126                 default:
 127                 return kSpUndefined;
 128         }
 129 }
 130
 131 Int_t AliSpectraBothPID::GetParticleSpecie(TParticle* part)
 132 {
 133   // return PID according to MC truth
 134   switch(TMath::Abs(part->GetPdgCode()))
 135   {
 136         case 2212:
 137         return kSpProton;
 138          break;
 139         case 321:
 140         return kSpKaon;
 141         break;
 142         case 211:
 143         return kSpPion;
 144         break;
 145         default:
 146         return kSpUndefined;
 147   }
 148 }
 149
 150 Int_t AliSpectraBothPID::GetParticleSpecie(AliSpectraBothHistoManager * hman,AliVTrack      * trk, AliSpectraBothTrackCuts * trackCuts, Bool_t* rec)
 151 {
 152   // return PID according to detectors
 153   // Get PID response object, if needed
 154
 155           // guess the particle based on the smaller nsigma
 156
 157         rec[kSpPion]=false;
 158          rec[kSpKaon]=false;
 159         rec[kSpProton]=false;
 160
 161         if(!fPIDResponse)
 162          {
 163                 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 164                 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 165                  fPIDResponse = inputHandler->GetPIDResponse();
 166         }
 167
 168         if(!fPIDResponse)
 169         {
 170                 AliFatal("Cannot get pid response");
 171         return 0;
 172         }
 173
 174
 175   // Compute nsigma for each hypthesis
 176         AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
 177
 178   // --- TPC
 179         Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
 180         Double_t nsigmaTPCkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
 181          Double_t nsigmaTPCkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
 182
 183
 184
 185         Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
 186
 187         Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 188         Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 189         Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 190         if(fPIDType==kNSigmaTOF)
 191         {
 192                 nsigmaTPCTOFkProton = 100.0;
 193                 nsigmaTPCTOFkKaon   = 100.0;
 194                 nsigmaTPCTOFkPion   =100.0;
 195
 196         }
 197
 198
 199
 200         if(trk->Pt()>trackCuts->GetPtTOFMatching())
 201         {
 202                 nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
 203                 nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
 204                 nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
 205                 // the TOF info is used in combined
 206                 nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 207                 nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 208                 nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 209         }
 210          else
 211         {
 212                 if (fPIDType==kNSigmaTOF)
 213                         return kSpUndefined;
 214
 215         }
 216
 217         // select the nsigma to be used for the actual PID
 218         Double_t nsigmaPion=100;
 219         Double_t nsigmaKaon=100;
 220         Double_t nsigmaProton=100;
 221
 222         switch (fPIDType)
 223         {
 224                 case kNSigmaTPC:
 225                 nsigmaProton  =  nsigmaTPCkProton;
 226                 nsigmaKaon        =  nsigmaTPCkKaon  ;
 227                 nsigmaPion    =  nsigmaTPCkPion  ;
 228                 break;
 229                 case kNSigmaTOF:
 230                 nsigmaProton  =  nsigmaTOFkProton;
 231                 nsigmaKaon        =  nsigmaTOFkKaon  ;
 232                 nsigmaPion    =  nsigmaTOFkPion  ;
 233                 break;
 234                 case kNSigmaTPCTOF:
 235                 nsigmaProton  =  nsigmaTPCTOFkProton;
 236                 nsigmaKaon        =  nsigmaTPCTOFkKaon  ;
 237                 nsigmaPion    =  nsigmaTPCTOFkPion  ;
 238                 break;
 239         }
 240
 241         if(nsigmaPion   < fNSigmaPID)
 242                 rec[kSpPion]=true;
 243         if(nsigmaKaon   < fNSigmaPID)
 244                 rec[kSpKaon]=true;
 245          if(nsigmaProton   < fNSigmaPID)
 246                 rec[kSpProton]=true;
 247
 248         if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
 249                 return kSpUndefined;//if is the default value for the three
 250         if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID))
 251         {
 252                 hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
 253                 return kSpKaon;
 254         }
 255         if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID))
 256         {
 257                 hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
 258                 return kSpPion;
 259         }
 260         if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
 261         {
 262                 hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
 263                 return kSpProton;
 264         }
 265         // else, return undefined
 266         return kSpUndefined;
 267
 268 }
 269
 270 Long64_t AliSpectraBothPID::Merge(TCollection* list)
 271 {
 272   // Merging interface.
 273   // Returns the number of merged objects (including this).
 274
 275   Printf("Merging");
 276
 277   if (!list)
 278     return 0;
 279
 280   if (list->IsEmpty())
 281     return 1;
 282
 283   TIterator* iter = list->MakeIterator();
 284   TObject* obj;
 285
 286   // Actually, we don't do anything here...
 287   // collections of all histograms
 288   //  TList collections;
 289
 290   Int_t count = 0;
 291
 292   while ((obj = iter->Next())) {
 293     AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
 294     if (entry == 0)
 295       continue;
 296
 297     // TH1I * histo = entry->GetHistoCuts();
 298     // collections.Add(histo);
 299     count++;
 300   }
 301
 302   //  fHistoCuts->Merge(&collections);
 303
 304   delete iter;
 305   Printf("OK");
 306   return count+1;
 307 }
 308