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         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);
 100                 nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2);
 101                 nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2);
 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=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
 203
 204         Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 205         Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 206         Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 207         if(fPIDType==kNSigmaTOF)
 208         {
 209                 nsigmaTPCTOFkProton = 100.0;
 210                 nsigmaTPCTOFkKaon   = 100.0;
 211                 nsigmaTPCTOFkPion   =100.0;
 212
 213         }
 214
 215
 216
 217         if(trk->Pt()>trackCuts->GetPtTOFMatching())
 218         {
 219                 nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
 220                 nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
 221                 nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
 222                 // the TOF info is used in combined
 223                 nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 224                 nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 225                 nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 226         }
 227          else
 228         {
 229                 if (fPIDType==kNSigmaTOF)
 230                         return kSpUndefined;
 231
 232         }
 233
 234         // select the nsigma to be used for the actual PID
 235         Double_t nsigmaPion=100;
 236         Double_t nsigmaKaon=100;
 237         Double_t nsigmaProton=100;
 238
 239         switch (fPIDType)
 240         {
 241                 case kNSigmaTPC:
 242                 nsigmaProton  =  nsigmaTPCkProton;
 243                 nsigmaKaon        =  nsigmaTPCkKaon  ;
 244                 nsigmaPion    =  nsigmaTPCkPion  ;
 245                 break;
 246                 case kNSigmaTOF:
 247                 nsigmaProton  =  nsigmaTOFkProton;
 248                 nsigmaKaon        =  nsigmaTOFkKaon  ;
 249                 nsigmaPion    =  nsigmaTOFkPion  ;
 250                 break;
 251                 case kNSigmaTPCTOF:
 252                 nsigmaProton  =  nsigmaTPCTOFkProton;
 253                 nsigmaKaon        =  nsigmaTPCTOFkKaon  ;
 254                 nsigmaPion    =  nsigmaTPCTOFkPion  ;
 255                 break;
 256         }
 257
 258         if(nsigmaPion   < fNSigmaPID)
 259                 rec[kSpPion]=true;
 260         if(nsigmaKaon   < fNSigmaPID)
 261                 rec[kSpKaon]=true;
 262          if(nsigmaProton   < fNSigmaPID)
 263                 rec[kSpProton]=true;
 264
 265         if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
 266                 return kSpUndefined;//if is the default value for the three
 267         if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID))
 268         {
 269                 hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
 270                 return kSpKaon;
 271         }
 272         if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID))
 273         {
 274                 hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
 275                 return kSpPion;
 276         }
 277         if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
 278         {
 279                 hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
 280                 return kSpProton;
 281         }
 282         // else, return undefined
 283         return kSpUndefined;
 284
 285 }
 286
 287 Long64_t AliSpectraBothPID::Merge(TCollection* list)
 288 {
 289   // Merging interface.
 290   // Returns the number of merged objects (including this).
 291
 292   Printf("Merging");
 293
 294   if (!list)
 295     return 0;
 296
 297   if (list->IsEmpty())
 298     return 1;
 299
 300   TIterator* iter = list->MakeIterator();
 301   TObject* obj;
 302
 303   // Actually, we don't do anything here...
 304   // collections of all histograms
 305   //  TList collections;
 306
 307   Int_t count = 0;
 308
 309   while ((obj = iter->Next())) {
 310     AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
 311     if (entry == 0)
 312       continue;
 313
 314     // TH1I * histo = entry->GetHistoCuts();
 315     // collections.Add(histo);
 316     count++;
 317   }
 318
 319   //  fHistoCuts->Merge(&collections);
 320
 321   delete iter;
 322   Printf("OK");
 323   return count+1;
 324 }
 325