PWGCF/FEMTOSCOPY/AliFemtoUser/AliAnalysisTaskParticleEfficiency.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TCanvas.h"
   6 #include "TList.h"
   7 #include "TObjArray.h"
   8 #include "TString.h"
   9 #include "TParticle.h"
  10
  11 #include "AliAnalysisTask.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliAnalysisTaskSE.h"
  14 #include "AliCentrality.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliESDtrackCuts.h"
  18 #include "AliESDInputHandler.h"
  19
  20 #include "AliAODEvent.h"
  21 #include "AliAODTrack.h"
  22 #include "AliAODHandler.h"
  23 #include "AliAODInputHandler.h"
  24 #include "AliAODMCParticle.h"
  25 #include "AliAODpidUtil.h"
  26 #include "AliAnalysisUtils.h"
  27
  28 #include "AliMCEvent.h"
  29 #include "AliStack.h"
  30 #include "AliInputEventHandler.h"
  31
  32 #include "AliAnalysisTaskParticleEfficiency.h"
  33
  34
  35 ClassImp(AliAnalysisTaskParticleEfficiency)
  36 //ClassImp(AliAnalysisTaskParticleEfficiency)
  37
  38 using std::cout;
  39 using std::endl;
  40
  41 //_______________________________________________________
  42
  43 AliAnalysisTaskParticleEfficiency::AliAnalysisTaskParticleEfficiency(const Char_t *partName) :
  44 AliAnalysisTaskSE(partName), centrality(0), fHistoList(0),  fHistEv(0), fpidResponse(0)
  45 {
  46
  47   for(Int_t i = 0; i < MULTBINS*PARTTYPES; i++)  {
  48     fGeneratedMCPrimaries[i] = NULL;
  49     fMCPrimariesThatAreReconstructed[i] = NULL;
  50     fReconstructedAfterCuts[i] = NULL;
  51     fReconstructedNotPrimaries[i] = NULL;
  52     fReconstructedPrimaries[i] = NULL;
  53     fContamination[i] = NULL;
  54   }
  55   for ( Int_t i = 0; i < 11; i++) {
  56     fHistQA[i] = NULL;
  57     if(i<3) fHistQA2D[i] = NULL;
  58   }
  59
  60   /* init track cuts */
  61   //fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
  62   /*fTrackCuts =  AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  63     if( !fTrackCuts ) return;
  64     fTrackCuts->SetMinNClustersTPC(70);*/
  65
  66
  67   /* create output */
  68   fHistoList = new TList();
  69   fHistoList->SetOwner(kTRUE);
  70
  71   DefineInput(0, TChain::Class());
  72   // DefineOutput(0, TTree::Class());
  73   DefineOutput(1, TList::Class());
  74
  75
  76 }
  77
  78 //_______________________________________________________
  79
  80 AliAnalysisTaskParticleEfficiency::~AliAnalysisTaskParticleEfficiency()
  81 {
  82   /* if(centrality) delete centrality;
  83 s128     if(fHistoList) delete fHistoList;
  84      if(vertex) delete vertex;
  85      if(vtxSPD) delete vtxSPD;*/
  86 }
  87
  88 //_______________________________________________________
  89
  90 void AliAnalysisTaskParticleEfficiency::UserCreateOutputObjects()
  91 {
  92
  93   TString hname1, hname2, hname3, hname4, hname5;
  94
  95   TString htitle1, htitle2, htitle3, htitle4,htitle5;
  96
  97   TString hname1M, hname2M, hname3M, hname4M, hname5M, hname;
  98
  99   TString htitle1M, htitle2M, htitle3M, htitle4M, htitle5M, htitle;
 100
 101   TString parttypename = "None";
 102
 103   for(Int_t j = 0; j < PARTTYPES; j++)  {
 104     if (j==0) parttypename="All";
 105     else if (j==1) parttypename="Pion";
 106     else if (j==2) parttypename="Kaon";
 107     else if (j==3) parttypename="Proton";
 108
 109     for(Int_t i = 0; i < MULTBINS; i++)  {
 110       hname1  = "hGeneratedMCPrimariesEffM"; hname1+=i; hname1+=parttypename;
 111       htitle1 = "Kinematic level eta_pT (prim only) M"; htitle1+=i; htitle1+=parttypename;
 112       fGeneratedMCPrimaries[i*PARTTYPES+j] = new TH2F(hname1.Data(),htitle1.Data(),50, -1.5, 1.5,1000,0.,10.0);
 113
 114       hname3  = "hMCPrimariesThatAreReconstructedM"; hname3+=i; hname3+=parttypename;
 115       htitle3 = "Reconstructed level eta_pT (prim only) M"; htitle3+=i; htitle3+=parttypename;
 116       fMCPrimariesThatAreReconstructed[i*PARTTYPES+j] = new TH2F(hname3.Data(),htitle3.Data(),50, -1.5, 1.5,1000,0.,10.0);
 117
 118       hname2  = "hHistoReconstructedAfterCutsM"; hname2+=i; hname2+=parttypename;
 119       htitle2 = "Total Reconstructed tracks M "; htitle2+=i; htitle2+=parttypename;
 120       fReconstructedAfterCuts[i*PARTTYPES+j] = new TH2F(hname2.Data(),htitle2.Data(),50, -1.5, 1.5,1000,0.,10.0);
 121
 122       hname4  = "hHistoReconstructedNotPrimariesM"; hname4+=i; hname4+=parttypename;
 123       htitle4 = "Reconstructed level eta_pT (not primaries) M"; htitle4+=i; htitle4+=parttypename;
 124       fReconstructedNotPrimaries[i*PARTTYPES+j] = new TH2F(hname4.Data(),htitle4.Data(),50, -1.5, 1.5,1000,0.,10.0);
 125
 126       hname4  = "hHistoReconstructedPrimariesM"; hname4+=i; hname4+=parttypename;
 127       htitle4 = "Reconstructed level eta_pT (primaries) M"; htitle4+=i; htitle4+=parttypename;
 128       fReconstructedPrimaries[i*PARTTYPES+j] = new TH2F(hname4.Data(),htitle4.Data(),50, -1.5, 1.5,1000,0.,10.0);
 129
 130       hname5  = "hContaminationM"; hname5+=i; hname5+=parttypename;
 131       htitle5 = "Contamination M"; htitle5+=i; htitle5+=parttypename;
 132       fContamination[i*PARTTYPES+j] = new TH2F(hname5.Data(),htitle5.Data(),6000, -3000, 3000.,1000,0.,10.0);
 133
 134
 135       fReconstructedAfterCuts[i*PARTTYPES+j]->Sumw2();
 136       fReconstructedNotPrimaries[i*PARTTYPES+j]->Sumw2();
 137       fReconstructedPrimaries[i*PARTTYPES+j]->Sumw2();
 138     }
 139     hname  = "pidTPCdEdx";  hname+=parttypename;
 140     htitle = parttypename + " TPC dEdx vs. momentum";
 141     fHistQAPID[0][j] = new TH2F(hname, htitle, 100, 0.0, 5.0, 250, 0.0, 500.0);
 142     hname  = "pidTOFTime";  hname+=parttypename;
 143     htitle = parttypename + " TOF Time vs. momentum";
 144     fHistQAPID[1][j] = new TH2F(hname, htitle, 100, 0.1, 5.0, 400, -4000.0, 4000.0);
 145     hname  = "pidTOFNSigma";  hname+=parttypename;
 146     htitle = parttypename + " TOF NSigma vs. momentum";
 147     fHistQAPID[2][j] = new TH2F(hname,htitle, 100, 0.0, 5.0, 100, -5.0, 5.0);
 148     hname  = "pidTPCNSigma";  hname+=parttypename;
 149     htitle = parttypename + " TPC NSigma vs. momentum";
 150     fHistQAPID[3][j] = new TH2F(hname,htitle, 100, 0.0, 5.0, 100, -5.0, 5.0);
 151     hname  = "pidTPCTOFNSigma";  hname+=parttypename;
 152     htitle = parttypename + " TPC vs TOF NSigma";
 153     fHistQAPID[4][j] = new TH2F(hname,htitle, 200, -10.0, 10.0, 200, -10.0, 10.0);
 154
 155   }
 156
 157   fHistEv = new TH1F("fHistEv", "Multiplicity", 100, 0, 100);
 158   fHistoList->Add(fHistEv);
 159
 160   for(Int_t i = 0; i < MULTBINS; i++)  {
 161     hname = "fHistEventCutsM";
 162     hname+= i;
 163
 164     fHistEvCuts[i] = new TH1F(hname,Form("Event Cuts M%d",i) , 4, 0, 5);
 165     fHistEvCuts[i]->GetXaxis()->SetBinLabel(1,"All");
 166     fHistEvCuts[i]->GetXaxis()->SetBinLabel(2,"NoVertex");
 167     fHistEvCuts[i]->GetXaxis()->SetBinLabel(3,"PileUp");
 168     fHistEvCuts[i]->GetXaxis()->SetBinLabel(4,"z-vertex>10");
 169     fHistoList->Add(fHistEvCuts[i]);
 170
 171
 172     hname  = "hMisidentificationM"; hname+=i;
 173     htitle = "Misidentification Fraction M"; htitle+=i;
 174     fMisidentification[i] = new TH2F(hname.Data(),htitle.Data(), 3, 0.5, 3.5, 4 , 0, 4);
 175     fMisidentification[i]->GetXaxis()->SetBinLabel(1,"Pions, MC");
 176     fMisidentification[i]->GetXaxis()->SetBinLabel(2,"Kaons, MC");
 177     fMisidentification[i]->GetXaxis()->SetBinLabel(3,"Protons, MC");
 178     fMisidentification[i]->GetYaxis()->SetBinLabel(1,"Pions, Data");
 179     fMisidentification[i]->GetYaxis()->SetBinLabel(2,"Kaons, Data");
 180     fMisidentification[i]->GetYaxis()->SetBinLabel(3,"Protons, Data");
 181     fMisidentification[i]->GetYaxis()->SetBinLabel(4,"Other, Data");
 182     fHistoList->Add(fMisidentification[i]);
 183
 184   }
 185
 186   fHistQA[0] = new TH1F("fHistVtx", "Z vertex distribution", 100, -15., 15.);
 187   fHistQA[1] = new TH1F("fHistnTpcCluster", "n TPC Cluster", 100, 0., 200.);
 188   fHistQA[2] = new TH1F("fHistnTpcClusterF", "n TPC Cluster findable", 100, 0., 200.);
 189   fHistQA[3] = new TH1F("dcaHistDcaXY1D", "DCA XY", 210, -2.1, 2.1);
 190   fHistQA[4] = new TH1F("dcaHistDcaZ1D", "DCA Z", 210, -2.1, 2.1);
 191   fHistQA[5] = new TH1F("fHistChi2Tpc", "Chi2 TPC", 100, 0., 8.);
 192   fHistQA[6] = new TH1F("fHistpT", "pT distribution",1000,0.,10.0);
 193   fHistQA[7] = new TH1F("fHistPhi", "Phi distribution" , 100, -TMath::Pi(), TMath::Pi());
 194   fHistQA[8] = new TH1F("fHistEta", "Eta distribution" , 100, -2, 2);
 195
 196   fHistQA[9] = new TH1F("fHistEventCuts", "Event Cuts" , 4, 0, 5);
 197   fHistQA[9]->GetXaxis()->SetBinLabel(1,"All");
 198   fHistQA[9]->GetXaxis()->SetBinLabel(2,"NoVertex");
 199   fHistQA[9]->GetXaxis()->SetBinLabel(3,"PileUp");
 200   fHistQA[9]->GetXaxis()->SetBinLabel(4,"z-vertex>10");
 201
 202
 203   fHistQA[10] = new TH1F("fHistTrackCuts", "Track Cuts" , 7, 0.5, 7.5);
 204   fHistQA[10]->GetXaxis()->SetBinLabel(1,"AllTracksInEvents");
 205   fHistQA[10]->GetXaxis()->SetBinLabel(2,"GetTrack");
 206   fHistQA[10]->GetXaxis()->SetBinLabel(3,"Filter bit");
 207   fHistQA[10]->GetXaxis()->SetBinLabel(4,"Eta");
 208   fHistQA[10]->GetXaxis()->SetBinLabel(5,"Pt");
 209   fHistQA[10]->GetXaxis()->SetBinLabel(6,"DCA");
 210   fHistQA[10]->GetXaxis()->SetBinLabel(7,"Electron Rejection");
 211
 212   fHistQA2D[0] = new TH2F("dcaHistDcaXY","DCA XY",50, 0, 5,210, -2.1, 2.1);
 213   fHistQA2D[1] = new TH2F("dcaHistDcaZ","DCA Z", 50, 0, 5, 210, -2.1, 2.1);
 214   fHistQA2D[2] = new TH2F("fPhiEta","Eta-Phi",100, -2, 2, 100, -TMath::Pi(), TMath::Pi());
 215
 216
 217
 218   for ( Int_t i = 0; i < 11; i++)
 219     {
 220       fHistoList->Add(fHistQA[i]);
 221       if(i<3) fHistoList->Add(fHistQA2D[i]);
 222       if(i<5) {
 223         for(Int_t j = 0 ; j<PARTTYPES; j++)
 224           fHistoList->Add(fHistQAPID[i][j]);
 225       }
 226     }
 227
 228   for (Int_t i = 0; i < MULTBINS*PARTTYPES; i++){
 229     fHistoList->Add(fGeneratedMCPrimaries[i]);
 230     fHistoList->Add(fMCPrimariesThatAreReconstructed[i]);
 231     fHistoList->Add(fReconstructedAfterCuts[i]);
 232     fHistoList->Add(fReconstructedNotPrimaries[i]);
 233     fHistoList->Add(fReconstructedPrimaries[i]);
 234     fHistoList->Add(fContamination[i]);
 235
 236   }
 237
 238   //********** PID ****************
 239
 240   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 241   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 242   fpidResponse = inputHandler->GetPIDResponse();
 243   cout<<"*******"<< fpidResponse<<endl;
 244
 245   // ************************
 246
 247   PostData(1, fHistoList);
 248 }
 249
 250
 251 //_____________________________________________________________________
 252
 253 bool IsPionNSigma(float mom, float nsigmaTPCPi, float nsigmaTOFPi)
 254 {
 255
 256     if (mom > 0.5) {
 257         if (TMath::Hypot( nsigmaTOFPi, nsigmaTPCPi ) < 3)
 258             return true;
 259         }
 260     else {
 261         if (TMath::Abs(nsigmaTPCPi) < 3)
 262             return true;
 263     }
 264 /*
 265   if(mom<0.65)
 266     {
 267       if(nsigmaTOFPi<-999.)
 268         {
 269           if(mom<0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
 270           else if(mom<0.5 && mom>=0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
 271           else if(mom>=0.5 && TMath::Abs(nsigmaTPCPi)<2.0) return true;
 272           else return false;
 273         }
 274       else if(TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
 275     }
 276   else if(nsigmaTOFPi<-999.)
 277     {
 278       return false;
 279     }
 280   else if(mom<1.5 && TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
 281   else if(mom>=1.5 && TMath::Abs(nsigmaTOFPi)<2.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
 282 */
 283   return false;
 284 }
 285
 286 bool IsKaonNSigma(float mom, float nsigmaTPCK, float nsigmaTOFK)
 287 {
 288     if (mom > 0.5) {
 289         if (TMath::Hypot( nsigmaTOFK, nsigmaTPCK ) < 3)
 290             return true;
 291         }
 292     else {
 293         if (TMath::Abs(nsigmaTPCK) < 3)
 294             return true;
 295     }
 296
 297  /* if(mom<0.4)
 298     {
 299       if(nsigmaTOFK<-999.)
 300         {
 301           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
 302         }
 303       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
 304     }
 305   else if(mom>=0.4 && mom<=0.6)
 306     {
 307       if(nsigmaTOFK<-999.)
 308         {
 309           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
 310         }
 311       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
 312     }
 313   else if(nsigmaTOFK<-999.)
 314     {
 315       return false;
 316     }
 317   else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
 318 */
 319   return false;
 320 }
 321
 322 bool IsProtonNSigma(float mom, float nsigmaTPCP, float nsigmaTOFP)
 323 {
 324
 325     if (mom > 0.5) {
 326         if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP ) < 3)
 327             return true;
 328         }
 329     else {
 330         if (TMath::Abs(nsigmaTPCP) < 3)
 331             return true;
 332     }
 333
 334   // if (mom > 0.8 && mom < 2.5) {
 335   //   if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 3.0)
 336   //     return true;
 337   // }
 338   // else if (mom > 2.5) {
 339   //   if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 2.0)
 340   //     return true;
 341   // }
 342   // else {
 343   //   if (TMath::Abs(nsigmaTPCP) < 3.0)
 344   //     return true;
 345   // }
 346
 347
 348   return false;
 349 }
 350
 351
 352 bool IsElectron(float nsigmaTPCE, float nsigmaTPCPi,float nsigmaTPCK, float nsigmaTPCP)
 353 {
 354   if(TMath::Abs(nsigmaTPCE)<3 && TMath::Abs(nsigmaTPCPi)>3 && TMath::Abs(nsigmaTPCK)>3 && TMath::Abs(nsigmaTPCP)>3)
 355       return true;
 356    else
 357      return false;
 358 }
 359
 360 //_______________________________________________________
 361
 362 void AliAnalysisTaskParticleEfficiency::UserExec(Option_t *)
 363 {
 364
 365   /***Get Event****/
 366   //AliESDEvent *esdEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 367   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 368   if (!aodEvent) return;
 369   AliAODHeader *fAODheader = dynamic_cast<AliAODHeader*>(aodEvent->GetHeader());
 370   if(!fAODheader) AliFatal("Not a standard AOD");
 371   Double_t mult = fAODheader->GetRefMultiplicity();
 372 // AliCentrality* alicent= aodEvent->GetCentrality(); //in PbPb and pPb
 373 //  Double_t mult = alicent->GetCentralityPercentile("V0A"); //in pPb
 374 //  Double_t mult = alicent->GetCentralityPercentile("V0A"); //in PbPb
 375   fHistEv->Fill(mult);
 376
 377
 378   //******************
 379   // load MC array
 380   // arrayMC =  (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 381   //  if(!arrayMC) {
 382   //  printf("AliAnalysisTaskParticleEficiency::UserExec: MC particles branch not found!\n");
 383   //return;
 384   //}
 385
 386   // load MC header
 387   //mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 388   //if(!mcHeader) {
 389   //printf("AliAnalysisTaskSEDplusCorrelations::UserExec: MC header branch not found!\n");
 390   //return;
 391   // }
 392   //*********************
 393
 394
 395
 396   // EVENT SELECTION ********************
 397
 398   fHistQA[9]->Fill(1);
 399
 400   // collision candidate
 401   // if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
 402
 403
 404   //****** Multiplicity selection *********
 405   Int_t fcent = -999;
 406   //if(mult >= 0 && mult <=20)  fcent = 0;
 407   //else if(mult >= 20 && mult <=39) fcent = 1;
 408   //else if(mult >= 40 && mult <=59) fcent = 2;
 409   //else if(mult >= 60 && mult <=90) fcent = 3;
 410   //else if(mult >= 99990 && mult <=99936) fcent = 4;
 411   //else if(mult >= 999937 && mult <=99944) fcent = 5;
 412   //else if(mult >= 999945 && mult <=99957) fcent = 6;
 413   //else if(mult >= 999958 && mult <=99149) fcent = 6;
 414   //else fcent = 7;
 415   //if (fcent == 7) return;
 416
 417   if(mult >= 2&& mult <=150)  fcent = 0;
 418   else if(mult >= 2 && mult <=19) fcent = 1;
 419   else if(mult >= 20 && mult <=49) fcent = 2;
 420   else if(mult >= 50 && mult <=150) fcent = 3;
 421   else return;
 422
 423   if(fcent==0)fHistEvCuts[0]->Fill(1);
 424   else if(fcent==1)fHistEvCuts[1]->Fill(1);
 425   else if(fcent==2)fHistEvCuts[2]->Fill(1);
 426   else if(fcent==3)fHistEvCuts[3]->Fill(1);
 427
 428   //"ESDs/pass2/AOD049/*AliAOD.root");
 429   const AliAODVertex* vertex =(AliAODVertex*) aodEvent->GetPrimaryVertex();
 430   if (!vertex || vertex->GetNContributors()<=0) return;
 431
 432   fHistQA[9]->Fill(2);
 433  if(fcent==0)fHistEvCuts[0]->Fill(2);
 434   else if(fcent==1)fHistEvCuts[1]->Fill(2);
 435   else if(fcent==2)fHistEvCuts[2]->Fill(2);
 436   else if(fcent==3)fHistEvCuts[3]->Fill(2);
 437
 438   AliAnalysisUtils *anaUtil=new AliAnalysisUtils();
 439
 440   Bool_t fpA2013 = kFALSE;
 441   Bool_t fMVPlp = kFALSE;
 442   Bool_t fisPileUp = kTRUE;
 443   Int_t fMinPlpContribMV = 0;
 444   Int_t fMinPlpContribSPD = 3;
 445
 446   if(fpA2013)
 447     if(anaUtil->IsVertexSelected2013pA(aodEvent)==kFALSE) return;
 448
 449   if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE);
 450   else anaUtil->SetUseMVPlpSelection(kFALSE);
 451
 452   if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV);
 453   if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD);
 454
 455   if(fisPileUp)
 456     if(anaUtil->IsPileUpEvent(aodEvent)) return;
 457
 458   delete anaUtil;
 459
 460  fHistQA[9]->Fill(3);
 461   if(fcent==0)fHistEvCuts[0]->Fill(3);
 462   else if(fcent==1)fHistEvCuts[1]->Fill(3);
 463   else if(fcent==2)fHistEvCuts[2]->Fill(3);
 464   else if(fcent==3)fHistEvCuts[3]->Fill(3);
 465
 466   //TString vtxTtl = vertex->GetTitle();
 467   //if (!vtxTtl.Contains("VertexerTracks")) return;
 468   Float_t zvtx = vertex->GetZ();
 469   if (TMath::Abs(zvtx) > 10) return;
 470   fHistQA[0]->Fill(zvtx);
 471   fHistQA[9]->Fill(4);
 472   if(fcent==0)fHistEvCuts[0]->Fill(4);
 473   else if(fcent==1)fHistEvCuts[1]->Fill(4);
 474   else if(fcent==2)fHistEvCuts[2]->Fill(4);
 475   else if(fcent==3)fHistEvCuts[3]->Fill(4);
 476
 477  //**** getting MC array ******
 478   TClonesArray  *arrayMC;
 479
 480   arrayMC = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 481
 482
 483   //get stack
 484   //AliStack *mcStack = mcEvent->Stack();
 485   //if (!mcStack) return;
 486   //***********************
 487
 488
 489   // old vertex selection
 490   /*const AliESDVertex *vertex = esdEvent->GetPrimaryVertex();
 491     if (vertex->GetNContributors() < 1) return;
 492
 493     //z-vertex cut
 494     if (TMath::Abs(vertex->GetZ()) > 10.) return;
 495
 496
 497     const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();*/
 498   // Double_t zVertex = vtxSPD->GetZ();
 499
 500   //cout << "Event  Z vtx ==========> " << vertex->GetZ() <<endl;
 501   // centrality selection
 502   //AliCentrality *centrality = aodEvent->GetHeader()->GetCentralityP();
 503   //if (centrality->GetQuality() != 0) return;
 504   //Double_t cent = centrality->GetCentralityPercentileUnchecked("V0M");
 505   //if(cent < 0 || cent > 100.) return;
 506
 507
 508 //copying pid information for FB 128
 509   int labels[20000];
 510   for (int il=0; il<20000; il++) labels[il] = -1;
 511
 512   // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
 513   for (int i=0;i<aodEvent->GetNumberOfTracks();i++) {
 514     const AliAODTrack *aodtrack=dynamic_cast<const AliAODTrack*>(aodEvent->GetTrack(i));
 515     if(!aodtrack) AliFatal("Not a standard AOD");
 516     if (!aodtrack->TestFilterBit(128)) {
 517       if(aodtrack->GetID() < 0) continue;
 518       labels[aodtrack->GetID()] = i;
 519     }
 520   }
 521
 522   //RECONSTRUCTED TRACKS
 523
 524   TObjArray recoParticleArray[PARTTYPES];
 525
 526   fHistQA[10]->Fill(1,aodEvent->GetNumberOfTracks());
 527   //loop over AOD tracks
 528   for (Int_t iTracks = 0; iTracks < aodEvent->GetNumberOfTracks(); iTracks++) {
 529     //get track
 530
 531     //AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esdEvent),iTracks);
 532     AliAODTrack *track = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(iTracks));
 533     if(!track) AliFatal("Not a standard AOD");
 534     if (!track)continue;
 535     fHistQA[10]->Fill(2);
 536
 537     //UInt_t filterBit = (1 << (0));
 538     UInt_t filterBit = 96;
 539     if(!track->TestFilterBit(filterBit))continue;
 540
 541     //if(!track->IsHybridGlobalConstrainedGlobal())continue;
 542     //if((track->IsHybridGlobalConstrainedGlobal())==false)continue;
 543     // if(!track->IsHybridTPCConstrainedGlobal())continue;
 544     // if(!track->IsTPCConstrained())continue;
 545     //if(!track->IsGlobalConstrained())continue;
 546     //if((track->TestFilterMask(AliAODTrack::kTrkTPCOnly)==false))continue;//cut0_BIT(0)
 547
 548     //   if((track->IsHybridGlobalConstrainedGlobal())==false)
 549     //  continue;//def_BIT(272)
 550
 551     //if((track->TestFilterMask(AliAODTrack::kTrkGlobal)==false))continue;//cut1_BIT(5)
 552
 553     fHistQA[10]->Fill(3);
 554
 555     if(track->Eta() < -0.8 || track->Eta() > 0.8)
 556       continue;
 557     fHistQA[10]->Fill(4);
 558
 559     if (track->Pt() < 0.2 || track->Pt() > 20)
 560       continue;
 561     fHistQA[10]->Fill(5);
 562
 563     //single track cuts
 564     // if(track->Chi2perNDF() > 4.0) continue;
 565     // if(track->GetTPCNcls() < 70) continue;
 566
 567     //DCA
 568
 569     Double_t DCAXY;
 570     Double_t DCAZ;
 571     //  if(filterBit==(1 << (7))){
 572     DCAXY = -TMath::Abs(track->DCA());
 573     DCAZ = -TMath::Abs(track->ZAtDCA());
 574
 575       if(!(DCAXY==-999 || DCAZ==-999)){
 576         //if(TMath::Abs(DCAXY) > 0.0182 + 0.035*TMath::Power(track->Pt(), -1.01)) continue; //XY, Pt dep
 577         //no DCA cut
 578         //if(TMath::Abs(DCAXY) > 1000.0) {continue;} //XY
 579         //if(TMath::Abs(DCAZ) > 1000.0) {continue;} //Z
 580       }
 581     else {
 582       // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
 583       // const AliAODVertex* vertex = (AliAODVertex*) aodEvent->GetPrimaryVertex(); (already defined above)
 584       float vertexX  = -999.;
 585       float vertexY  = -999.;
 586       float vertexZ  = -999.;
 587
 588       if(vertex) {
 589         Double32_t fCov[6];
 590         vertex->GetCovarianceMatrix(fCov);
 591         if(vertex->GetNContributors() > 0) {
 592           if(fCov[5] != 0) {
 593             vertexX = vertex->GetX();
 594             vertexY = vertex->GetY();
 595             vertexZ = vertex->GetZ();
 596
 597           }
 598         }
 599       }
 600
 601       Double_t pos[3];
 602       track->GetXYZ(pos);
 603
 604       Double_t DCAX = pos[0] - vertexX;
 605       Double_t DCAY = pos[1] - vertexY;
 606       DCAZ = pos[2] - vertexZ;
 607       DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
 608
 609       //if(TMath::Abs(DCAXY) > 0.0182 + 0.035*TMath::Power(track->Pt(), -1.01)) continue; //XY, Pt dep
 610       //if(TMath::Abs(impactD) > 0.44 + 0.07*TMath::Power(tPt, -1.94)) continue; //XY, Pt dep
 611       //no DCA cut
 612       //if(TMath::Abs(DCAXY) > 1000.0) continue;
 613       //if(TMath::Abs(DCAZ) > 1000.0) continue;
 614     }
 615
 616     fHistQA[10]->Fill(6);
 617
 618     AliAODTrack* aodtrackpid;
 619
 620     //for FB 128 - tpc only tracks
 621     if(filterBit==(1 << (7))) {
 622       aodtrackpid = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(labels[-1-aodEvent->GetTrack(iTracks)->GetID()]));
 623       if(!aodtrackpid) AliFatal("Not a standard AOD");
 624     }
 625     else
 626       aodtrackpid = track;
 627
 628
 629  //PANOS--------------------------
 630
 631     AliMCEvent* mcEvent = MCEvent();
 632     if (!mcEvent) {
 633       AliError("ERROR: Could not retrieve MC event");
 634       return;
 635     }
 636
 637     Int_t nMCParticles = mcEvent->GetNumberOfTracks();
 638     Int_t labelp = TMath::Abs(track->GetLabel());
 639     if(labelp > nMCParticles) continue;
 640
 641     AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(labelp);
 642
 643     /*
 644     //contamination from secondaries without electron rejection
 645     if (!AODmcTrack->IsPhysicalPrimary()) {
 646       fReconstructedNotPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 647     }
 648     else{
 649       fReconstructedPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 650     }
 651     */
 652     //-------------------------
 653
 654     //Electron rejection
 655     double nSigmaTPCPi = fpidResponse->NumberOfSigmasTPC(aodtrackpid,AliPID::kPion);
 656     double nSigmaTPCK = fpidResponse->NumberOfSigmasTPC(aodtrackpid,AliPID::kKaon);
 657     double nSigmaTPCP = fpidResponse->NumberOfSigmasTPC(aodtrackpid,AliPID::kProton);
 658     double nSigmaTPCe = fpidResponse->NumberOfSigmasTPC(aodtrackpid,AliPID::kElectron);
 659     if(IsElectron(nSigmaTPCe,nSigmaTPCPi,nSigmaTPCK,nSigmaTPCP))
 660       continue;
 661     fHistQA[10]->Fill(7);
 662
 663     fHistQA[1]->Fill(track->GetTPCClusterInfo(2,1));
 664     //fHistQA[2]->Fill(track->GetTPCNclsF());
 665      fHistQA[3]->Fill(DCAXY);
 666      fHistQA[4]->Fill(DCAZ);
 667     Float_t chi2Tpc = track->Chi2perNDF();
 668     fHistQA[5]->Fill(chi2Tpc);
 669     fHistQA[6]->Fill(track->Pt());
 670
 671     float px=track->Px(); float py=track->Py();  float ph=atan2(py,px); //track->Phi()
 672     float tPt = track->Pt();
 673
 674     fHistQA[7]->Fill(ph);
 675     fHistQA[8]->Fill(track->Eta());
 676     fHistQA2D[2]->Fill(track->Eta(),ph);
 677
 678     fHistQA2D[0]->Fill(tPt,DCAXY);
 679     fHistQA2D[1]->Fill(tPt,DCAZ);
 680
 681
 682  //PANOS
 683     //contamination from secondaries with electron rejection
 684     // if (!AODmcTrack->IsPhysicalPrimary()) {
 685     //   fReconstructedNotPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 686     // }
 687     // else{
 688     //   fReconstructedPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 689     // }
 690
 691
 692     //PID monitors
 693     double nSigmaTOFPi = fpidResponse->NumberOfSigmasTOF(aodtrackpid,AliPID::kPion);
 694     double nSigmaTOFK = fpidResponse->NumberOfSigmasTOF(aodtrackpid,AliPID::kKaon);
 695     double nSigmaTOFP = fpidResponse->NumberOfSigmasTOF(aodtrackpid,AliPID::kProton);
 696
 697     float tdEdx = aodtrackpid->GetTPCsignal();
 698     float tTofSig = aodtrackpid->GetTOFsignal();
 699     double pidTime[5]; aodtrackpid->GetIntegratedTimes(pidTime);
 700
 701
 702     fHistQAPID[0][0]->Fill(tPt,tdEdx);
 703     fHistQAPID[1][0]->Fill(tPt,tTofSig-pidTime[2]);//pion
 704     fHistQAPID[2][0]->Fill(tPt,nSigmaTOFPi);
 705     fHistQAPID[3][0]->Fill(tPt,nSigmaTPCPi);
 706     fHistQAPID[4][0]->Fill(nSigmaTPCPi,nSigmaTOFPi);
 707
 708     if (IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi)){
 709       fHistQAPID[0][1]->Fill(tPt,tdEdx);
 710       fHistQAPID[1][1]->Fill(tPt,tTofSig-pidTime[2]);//pion
 711       fHistQAPID[2][1]->Fill(tPt,nSigmaTOFPi);
 712       fHistQAPID[3][1]->Fill(tPt,nSigmaTPCPi);
 713       fHistQAPID[4][1]->Fill(nSigmaTPCPi,nSigmaTOFPi);
 714     }
 715     if (IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK)){
 716       fHistQAPID[0][2]->Fill(tPt,tdEdx);
 717       fHistQAPID[1][2]->Fill(tPt,tTofSig-pidTime[3]);//kaon
 718       fHistQAPID[2][2]->Fill(tPt,nSigmaTOFK);
 719       fHistQAPID[3][2]->Fill(tPt,nSigmaTPCK);
 720       fHistQAPID[4][2]->Fill(nSigmaTPCK,nSigmaTOFK);
 721     }
 722     if (IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP)){
 723       fHistQAPID[0][3]->Fill(tPt,tdEdx);
 724       fHistQAPID[1][3]->Fill(tPt,tTofSig-pidTime[4]);//proton
 725       fHistQAPID[2][3]->Fill(tPt,nSigmaTOFP);
 726       fHistQAPID[3][3]->Fill(tPt,nSigmaTPCP);
 727       fHistQAPID[4][3]->Fill(nSigmaTPCP,nSigmaTOFP);
 728     }
 729
 730     fReconstructedAfterCuts[fcent]->Fill(track->Eta(), track->Pt());//Fills hist. for all reconstructed particles after cuts
 731
 732     if(!arrayMC){
 733       continue;
 734     }
 735     //get coresponding MC particle
 736     Int_t label = TMath::Abs(track->GetLabel());
 737     AliAODMCParticle *MCtrk = (AliAODMCParticle*)arrayMC->At(label);
 738
 739    //getting no. of tracks for each particle species after all the cuts:
 740
 741     //********* PID - pions ********
 742      if (IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi)){
 743        fReconstructedAfterCuts[PARTTYPES*fcent+1]->Fill(track->Eta(), track->Pt());
 744        if (!MCtrk) continue;
 745        recoParticleArray[1].Add(MCtrk);
 746        }
 747        //Fills for all identified pions found after cuts (reconstructed) - numerator for Efficiency
 748
 749      //********* PID - kaons ********
 750      if (IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK)){
 751        fReconstructedAfterCuts[PARTTYPES*fcent+2]->Fill(track->Eta(), track->Pt());
 752        if (!MCtrk) continue;
 753        recoParticleArray[2].Add(MCtrk);
 754        }
 755        //Fills for all identified kaons found after cuts (reconstructed) - numerator for Efficiency
 756
 757     //********* PID - protons ********
 758      if (IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP)){
 759        fReconstructedAfterCuts[PARTTYPES*fcent+3]->Fill(track->Eta(), track->Pt());
 760        if (!MCtrk) continue;
 761        recoParticleArray[3].Add(MCtrk);
 762        }
 763
 764       //Fills for all identified protos found after cuts (reconstructed) - numerator for Efficiency
 765    //******************************
 766
 767      //get coresponding MC particle
 768      // Int_t label = TMath::Abs(track->GetLabel()); //moved up
 769      // if(!label) cout<<"no label"<<endl;
 770      //if(label) cout<<"label = "<<label<<endl;
 771
 772     //AliAODMCParticle *MCtrk = (AliAODMCParticle*)arrayMC->At(label); //moved up
 773     if (!MCtrk) continue;
 774     if(MCtrk->Charge()==0){cout<<"!!!"<<endl; continue;}
 775     recoParticleArray[0].Add(MCtrk);
 776
 777
 778     //Fills histogram for particles that are contamination from secondaries:
 779     if (!AODmcTrack->IsPhysicalPrimary()) {
 780       fReconstructedNotPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 781     }
 782     else{
 783       fReconstructedPrimaries[PARTTYPES*fcent]->Fill(track->Eta(), track->Pt());
 784     }
 785
 786     int PDGcode = MCtrk->GetPdgCode();
 787
 788    //And secondaries for different particle species:
 789     if (!AODmcTrack->IsPhysicalPrimary() && (IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi) && PDGcode==211)) { //secondaries in pions
 790       fReconstructedNotPrimaries[PARTTYPES*fcent+1]->Fill(track->Eta(), track->Pt());
 791     }
 792     else if(AODmcTrack->IsPhysicalPrimary() && (IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi) && PDGcode==211)) {
 793       fReconstructedPrimaries[PARTTYPES*fcent+1]->Fill(track->Eta(), track->Pt());
 794     }
 795
 796     if (!AODmcTrack->IsPhysicalPrimary() && (IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK) && PDGcode==321)) { //secondaries in kaons
 797       fReconstructedNotPrimaries[PARTTYPES*fcent+2]->Fill(track->Eta(), track->Pt());
 798     }
 799     else if(AODmcTrack->IsPhysicalPrimary() && (IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK) && PDGcode==321)) {
 800       fReconstructedPrimaries[PARTTYPES*fcent+2]->Fill(track->Eta(), track->Pt());
 801       }
 802
 803     if (!AODmcTrack->IsPhysicalPrimary() && (IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP) && PDGcode==2212)) { //secondaries in protons
 804       fReconstructedNotPrimaries[PARTTYPES*fcent+3]->Fill(track->Eta(), track->Pt());
 805     }
 806     else if(AODmcTrack->IsPhysicalPrimary() && (IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP) && PDGcode==2212)) {
 807       fReconstructedPrimaries[PARTTYPES*fcent+3]->Fill(track->Eta(), track->Pt());
 808     }
 809
 810     //Misidentification fraction
 811     if(PDGcode==211)
 812       {
 813         if(IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi))
 814           fMisidentification[fcent]-> Fill(1,0.5);
 815         if(IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK))
 816           fMisidentification[fcent]-> Fill(1,1.5);
 817         if(IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 818           fMisidentification[fcent]-> Fill(1,2.5);
 819         if(!IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi) && !IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK) && !IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 820           fMisidentification[fcent]-> Fill(1,3.5);
 821
 822
 823       }
 824     else if(PDGcode==321)
 825       {
 826         if(IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi))
 827           fMisidentification[fcent]-> Fill(2,0.5);
 828         if(IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK))
 829           fMisidentification[fcent]-> Fill(2,1.5);
 830         if(IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 831           fMisidentification[fcent]-> Fill(2,2.5);
 832         if(!IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi) && !IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK) && !IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 833           fMisidentification[fcent]-> Fill(2,3.5);
 834
 835
 836       }
 837     else if(PDGcode == 2212)
 838       {
 839         if(IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi))
 840           fMisidentification[fcent]-> Fill(3,0.5);
 841         if(IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK))
 842           fMisidentification[fcent]-> Fill(3,1.5);
 843         if(IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 844           {
 845           fMisidentification[fcent]-> Fill(3,2.5);
 846           }
 847         if(!IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi) && !IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK) && !IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 848           fMisidentification[fcent]-> Fill(3,3.5);
 849       }
 850
 851
 852     //Contaminations: "how many pions are in the kaons sample"? etc.
 853     //Do not use for corrections: using those values will be dependant on i.e. Pi/K ratio in MC
 854     //Use misidentification fraction instead
 855     if(IsPionNSigma(track->P(),nSigmaTPCPi, nSigmaTOFPi))
 856       {
 857         fContamination[PARTTYPES*fcent+1]-> Fill(PDGcode,track->Pt()); // filling contamination histogram for pions
 858       }
 859     if(IsKaonNSigma(track->P(),nSigmaTPCK, nSigmaTOFK))
 860       {
 861         fContamination[PARTTYPES*fcent+2]-> Fill(PDGcode,track->Pt()); // filling contamination histogram for kaons
 862       }
 863     if(IsProtonNSigma(track->P(),nSigmaTPCP, nSigmaTOFP))
 864       {
 865         fContamination[PARTTYPES*fcent+3]-> Fill(PDGcode,track->Pt()); // filling contamination histogram for protons
 866       }
 867
 868
 869
 870   }
 871
 872
 873   // MONTECARLO PARTICLES
 874   if(!arrayMC){
 875     AliError("Array of MC particles not found");
 876     return;
 877   }
 878   // loop over MC stack
 879   for (Int_t ipart = 0; ipart < arrayMC->GetEntries(); ipart++) {
 880     //cout<<"Entered MC loop"<<endl;
 881
 882     AliAODMCParticle *MCtrk = (AliAODMCParticle*)arrayMC->At(ipart);
 883
 884     if (!MCtrk) continue;
 885     //cout<<"particle obtained"<<endl;
 886
 887     Int_t PDGcode = TMath::Abs(MCtrk->GetPdgCode());
 888
 889
 890     if(MCtrk->Charge() == 0) continue;
 891
 892     //*** PID - check if pion ***
 893     //if(PDGcode!=211) continue; //(PDGcode==11 || PDGcode==321 || PDGcode==2212 || PDGcode==13)
 894
 895       if(MCtrk->Eta() < -0.8 || MCtrk->Eta() > 0.8){
 896         continue; }
 897
 898       if (MCtrk->Pt() < 0.2 || MCtrk->Pt() > 20){
 899         continue;}
 900
 901       // check physical primary
 902       if (MCtrk->IsPhysicalPrimary()){
 903
 904         // Filling histograms for MC truth particles
 905         fGeneratedMCPrimaries[fcent*PARTTYPES]->Fill(MCtrk->Eta(), MCtrk->Pt());
 906         if(PDGcode==211)
 907           fGeneratedMCPrimaries[fcent*PARTTYPES+1]->Fill(MCtrk->Eta(), MCtrk->Pt());
 908         else if(PDGcode==321)
 909           fGeneratedMCPrimaries[fcent*PARTTYPES+2]->Fill(MCtrk->Eta(), MCtrk->Pt());
 910         else if(PDGcode==2212)
 911           fGeneratedMCPrimaries[fcent*PARTTYPES+3]->Fill(MCtrk->Eta(), MCtrk->Pt());
 912
 913           //Filling data from MC truth particles only for particles that were reconstruced
 914         if (recoParticleArray[0].Contains(MCtrk)){ //All
 915           fMCPrimariesThatAreReconstructed[fcent*PARTTYPES]->Fill(MCtrk->Eta(), MCtrk->Pt());
 916         }
 917         if (recoParticleArray[1].Contains(MCtrk)){ //Pions
 918           if(PDGcode==211)
 919             fMCPrimariesThatAreReconstructed[fcent*PARTTYPES+1]->Fill(MCtrk->Eta(), MCtrk->Pt());
 920         }
 921         if (recoParticleArray[2].Contains(MCtrk)){ //Kaons
 922           if(PDGcode==321)
 923             fMCPrimariesThatAreReconstructed[fcent*PARTTYPES+2]->Fill(MCtrk->Eta(), MCtrk->Pt());
 924         }
 925         if (recoParticleArray[3].Contains(MCtrk)){ //Protons
 926           if(PDGcode==2212)
 927             fMCPrimariesThatAreReconstructed[fcent*PARTTYPES+3]->Fill(MCtrk->Eta(), MCtrk->Pt());
 928         }
 929
 930       }
 931
 932   }
 933   PostData(1, fHistoList);
 934 }
 935 //-----------------------------------------------------------------
 936
 937 //void AliAnalysisTaskParticleEfficiency::Terminate(Option_t *)
 938 //{}