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         Float_t tmppt=TMath::Min(trackCuts->GetPtTOFMatchingPion(),TMath::Min(trackCuts->GetPtTOFMatchingKaon(),trackCuts->GetPtTOFMatchingProton()));
  72         if(track->Pt()>tmppt&&trackCuts->CheckTOFMatchingParticleType(kSpProton)&&trackCuts->CheckTOFMatchingParticleType(kSpKaon)&&trackCuts->CheckTOFMatchingParticleType(kSpPion))
  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;
 203         Double_t nsigmaTOFkKaon=0.0;
 204         Double_t nsigmaTOFkPion=0.0;
 205
 206
 207
 208         Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
 209         Double_t nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
 210         Double_t nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
 211         if(fPIDType==kNSigmaTOF)
 212         {
 213                 nsigmaTPCTOFkProton = 10000.0;
 214                 nsigmaTPCTOFkKaon   = 10000.0;
 215                 nsigmaTPCTOFkPion   =10000.0;
 216                 nsigmaTOFkProton = 10000.0;
 217                 nsigmaTOFkKaon   = 10000.0;
 218                 nsigmaTOFkPion   =10000.0;
 219
 220         }
 221
 222
 223         if(trackCuts->GetUseTypeDependedTOFCut())
 224         {
 225                 if(trackCuts->CheckTOFMatchingParticleType(kSpPion)&&trk->Pt()>trackCuts->GetPtTOFMatchingPion())
 226                 {
 227                         nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
 228                         nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 229                 }
 230                 else if (trk->Pt()>trackCuts->GetPtTOFMatchingPion()) // tracks without proper flags
 231                 {
 232                         nsigmaTOFkPion   = 10000.0;
 233                         nsigmaTPCTOFkPion   =10000.0;
 234
 235                 }
 236
 237                 if(trackCuts->CheckTOFMatchingParticleType(kSpKaon)&&trk->Pt()>trackCuts->GetPtTOFMatchingKaon())
 238                 {
 239                         nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
 240                         nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 241                 }
 242                 else if (trk->Pt()>trackCuts->GetPtTOFMatchingKaon()) // tracks without proper flags
 243                 {
 244                         nsigmaTOFkKaon   = 10000.0;
 245                         nsigmaTPCTOFkKaon   =10000.0;
 246
 247                 }
 248
 249                 if(trackCuts->CheckTOFMatchingParticleType(kSpProton)&&trk->Pt()>trackCuts->GetPtTOFMatchingProton())
 250                 {
 251                         nsigmaTOFkProton   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
 252                         nsigmaTPCTOFkProton   = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 253                 }
 254                 else if (trk->Pt()>trackCuts->GetPtTOFMatchingProton()) // tracks without proper flags
 255                 {
 256                         nsigmaTOFkProton   = 10000.0;
 257                         nsigmaTPCTOFkProton   =10000.0;
 258
 259                 }
 260
 261
 262         }
 263         else if(trk->Pt()>trackCuts->GetPtTOFMatching())
 264         {
 265                 nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
 266                 nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
 267                 nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
 268                 // the TOF info is used in combined
 269                 nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
 270                 nsigmaTPCTOFkKaon   = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
 271                 nsigmaTPCTOFkPion   = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
 272         }
 273          else
 274         {
 275                 if (fPIDType==kNSigmaTOF)
 276                         return kSpUndefined;
 277
 278         }
 279
 280         // select the nsigma to be used for the actual PID
 281         Double_t nsigmaPion=10000.0;
 282         Double_t nsigmaKaon=10000.0;
 283         Double_t nsigmaProton=10000.0;
 284
 285         switch (fPIDType)
 286         {
 287                 case kNSigmaTPC:
 288                 nsigmaProton  =  nsigmaTPCkProton;
 289                 nsigmaKaon        =  nsigmaTPCkKaon  ;
 290                 nsigmaPion    =  nsigmaTPCkPion  ;
 291                 break;
 292                 case kNSigmaTOF:
 293                 nsigmaProton  =  nsigmaTOFkProton;
 294                 nsigmaKaon        =  nsigmaTOFkKaon  ;
 295                 nsigmaPion    =  nsigmaTOFkPion  ;
 296                 break;
 297                 case kNSigmacircleTPCTOF:
 298                 nsigmaProton  =  nsigmaTPCTOFkProton;
 299                 nsigmaKaon        =  nsigmaTPCTOFkKaon  ;
 300                 nsigmaPion    =  nsigmaTPCTOFkPion  ;
 301                 break;
 302                 case kNSigmasquareTPCTOF:
 303                 nsigmaProton  =  TMath::Max(nsigmaTPCkProton,nsigmaTOFkProton);
 304                 nsigmaKaon        =  TMath::Max(nsigmaTPCkKaon,nsigmaTOFkKaon);
 305                 nsigmaPion    =  TMath::Max(nsigmaTPCkPion,nsigmaTOFkPion)  ;
 306                 break;
 307
 308         }
 309
 310         if(nsigmaPion   < fNSigmaPID)
 311                 rec[kSpPion]=true;
 312         if(nsigmaKaon   < fNSigmaPID)
 313                 rec[kSpKaon]=true;
 314          if(nsigmaProton   < fNSigmaPID)
 315                 rec[kSpProton]=true;
 316
 317         if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
 318                 return kSpUndefined;//if is the default value for the three
 319         if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID))
 320         {
 321                 if(hman->GetPIDHistogram(kHistPIDTPCKaonRec))
 322                         hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
 323                 return kSpKaon;
 324         }
 325         if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID))
 326         {
 327                 if(hman->GetPIDHistogram(kHistPIDTPCPionRec))
 328                         hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
 329                 return kSpPion;
 330         }
 331         if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
 332         {
 333                 if(hman->GetPIDHistogram(kHistPIDTPCProtonRec))
 334                         hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
 335                 return kSpProton;
 336         }
 337         // else, return undefined
 338         return kSpUndefined;
 339
 340 }
 341
 342 Long64_t AliSpectraBothPID::Merge(TCollection* list)
 343 {
 344   // Merging interface.
 345   // Returns the number of merged objects (including this).
 346
 347   Printf("Merging");
 348
 349   if (!list)
 350     return 0;
 351
 352   if (list->IsEmpty())
 353     return 1;
 354
 355   TIterator* iter = list->MakeIterator();
 356   TObject* obj;
 357
 358   // Actually, we don't do anything here...
 359   // collections of all histograms
 360   //  TList collections;
 361
 362   Int_t count = 0;
 363
 364   while ((obj = iter->Next())) {
 365     AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
 366     if (entry == 0)
 367       continue;
 368
 369     // TH1I * histo = entry->GetHistoCuts();
 370     // collections.Add(histo);
 371     count++;
 372   }
 373
 374   //  fHistoCuts->Merge(&collections);
 375
 376   delete iter;
 377   Printf("OK");
 378   return count+1;
 379 }
 380