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