Added analysis of fakes (Francesco)

diff --git a/PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx b/PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx
index 0e457c5960e052c0958956adb81c947c9dbd68cb..4d1f636dfcd99eb9336f5ae8101df76ee346e31f 100644 (file)
--- a/PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx
+++ b/PWG1/ITS/AliAnalysisTaskITSsaTracks.cxx
@@ -12,6 +12,7 @@
 #include <TParticle.h>
 #include <TSystem.h>
 #include <TTree.h>
+#include <TNtuple.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TChain.h>
@@ -48,25 +49,16 @@ ClassImp(AliAnalysisTaskITSsaTracks)
 AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("ITSsa resolution"), 
   fOutput(0),
   fHistNEvents(0),
-  fHistPtTPCITSAll(0),
-  fHistPtTPCITSGood(0),
-  fHistPtTPCITSFake(0),
-  fHistPtITSsaAll(0),
-  fHistPtITSsaGood(0),
-  fHistPtITSsaFake(0),
-  fHistPtITSpureSAAll(0),
-  fHistPtITSpureSAGood(0),
-  fHistPtITSpureSAFake(0),
-  fHistdedxvsPtITSpureSA3cls(0),
-  fHistdedxvsPITSpureSA3cls(0),
-  fHistdedxvsPtITSpureSA4cls(0),
-  fHistdedxvsPITSpureSA4cls(0),
+  fHistMCPhiResid(0),
+  fHistPhiResid(0),
+  fNtupleTracks(0),
   fNPtBins(kMaxPtBins),
   fMinITSpts(4),
   fMinSPDpts(1),
-  fMinITSptsForPid(3),
-  fMinITSptsForMatch(6),
+  fMinPtsforPid(3),
   fMinTPCpts(50),
+  fMaxITSChi2Clu(100.),
+  fFillNtuple(kFALSE),
   fReadMC(kFALSE),
   fUseMCId(kFALSE)
 {
@@ -85,6 +77,52 @@ AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("IT
 //___________________________________________________________________________
 AliAnalysisTaskITSsaTracks::~AliAnalysisTaskITSsaTracks(){
   //
+  delete fHistNEvents;
+  for(Int_t iType=0; iType<kNtrackTypes; iType++){
+    delete fHistPt[iType];
+    delete fHistPtGood[iType];
+    delete fHistPtFake[iType];
+    delete fHistEtaPhi[iType];
+    delete fHistEtaPhiGood[iType];
+    delete fHistEtaPhiFake[iType];
+    delete fHistChi2[iType];
+    delete fHistChi2Good[iType];
+    delete fHistChi2Fake[iType];
+    delete fHistNclu[iType];
+    delete fHistNcluGood[iType];
+    delete fHistNcluFake[iType];
+    delete fHistdedxvsPt3cls[iType];
+    delete fHistdedxvsPt4cls[iType];
+  }
+  for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
+    delete fHistPtTPCITS[iSpec];
+    delete fHistPtITSsa[iSpec];
+    delete fHistPtITSpureSA[iSpec];
+    delete fHistEtaPhiTPCITS[iSpec];
+    delete fHistEtaPhiITSsa[iSpec];
+    delete fHistEtaPhiITSpureSA[iSpec];
+    delete fHistNcluTPCITS[iSpec];
+    delete fHistNcluITSsa[iSpec];
+    delete fHistNcluITSpureSA[iSpec];
+    delete fHistd0rphiITSpureSA[iSpec];
+    delete fHistd0zITSpureSA[iSpec];
+    delete fHistCluInLayTPCITS[iSpec];
+    delete fHistCluInLayITSsa[iSpec];
+    delete fHistCluInLayITSpureSA[iSpec];
+    delete fHistOuterLayITSpureSA[iSpec];
+    delete fHistPtResid[iSpec];
+    delete fHistPtRelResid[iSpec];
+    delete fHistInvPtResid[iSpec];
+    delete fHistInvPtRelResid[iSpec];
+    delete fHistMCPtResid[iSpec];
+    delete fHistMCPtRelResid[iSpec];
+    delete fHistMCInvPtResid[iSpec];
+    delete fHistMCInvPtRelResid[iSpec];
+  } 
+
+  delete fHistMCPhiResid;
+  delete fHistPhiResid;
+  delete fNtupleTracks;
   if (fOutput) {
     delete fOutput;
     fOutput = 0;
@@ -111,6 +149,17 @@ void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
   fOutput->SetOwner();
   fOutput->SetName("OutputHistos");
 
+  fNtupleTracks= new TNtuple("ntupleTracks","ntupleTracks","pt:px:py:pz:eta:phi:d0xy:d0z:dedx3:dedx4:dedx5:dedx6:truncmean:chi2:status:ITSrefit:TPCin:TPCrefit:ITSpureSA:nITSclu:nTPCclu:clumap:spdin:spdout:sddin:sddout:ssdin:ssdout:label:ptgen:pxgen:pygen:pzgen:etagen:phigen:ntracks");
+  // kinematics: pt, p eta,phi        ->  6 variables
+  // impact parameters: d0xy, d0z     ->  2 variables
+  // dE/dx: 4 Layers + trunc mean     ->  5 variables
+  // chi2 and track status:           ->  6 variables
+  // cluster infos:                   ->  9 variables
+  // MC info:                         ->  7 variables
+  // Multiplicity                     ->  1 variable
+  // Total:                              36
+  fOutput->Add(fNtupleTracks);
+
   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
   fHistNEvents->Sumw2();
   fHistNEvents->SetMinimum(0);
@@ -129,59 +178,65 @@ void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
     hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
   }
   
-  fHistdedxvsPtITSpureSA3cls = new TH2F("hdedxvsPtITSpureSA3cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
-  fHistdedxvsPtITSpureSA3cls->Sumw2();
-  fOutput->Add(fHistdedxvsPtITSpureSA3cls);
-  
-  fHistdedxvsPITSpureSA3cls = new TH2F("hdedxvsPITSpureSA3cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
-  fHistdedxvsPITSpureSA3cls->Sumw2();
-  fOutput->Add(fHistdedxvsPITSpureSA3cls);
-  
-  fHistdedxvsPtITSpureSA4cls = new TH2F("hdedxvsPtITSpureSA4cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
-  fHistdedxvsPtITSpureSA4cls->Sumw2();
-  fOutput->Add(fHistdedxvsPtITSpureSA4cls);
-  
-  fHistdedxvsPITSpureSA4cls = new TH2F("hdedxvsPITSpureSA4cls","",hnbinsdedx,hxbinsdedx,900,0,1000);
-  fHistdedxvsPITSpureSA4cls->Sumw2();
-  fOutput->Add(fHistdedxvsPITSpureSA4cls);
-  
+  TString tit[kNtrackTypes]={"TPCITS","ITSsa","ITSpureSA"};
+  for(Int_t iType=0; iType<kNtrackTypes; iType++){
+    fHistPt[iType]=new TH1F(Form("hPt%s",tit[iType].Data()),"",100,0.,2.);
+    fHistPt[iType]->Sumw2();
+    fOutput->Add(fHistPt[iType]);
+    fHistPtGood[iType]=new TH1F(Form("hPtGood%s",tit[iType].Data()),"",100,0.,2.);
+    fHistPtGood[iType]->Sumw2();
+    fOutput->Add(fHistPtGood[iType]);
+    fHistPtFake[iType]=new TH1F(Form("hPtFake%s",tit[iType].Data()),"",100,0.,2.);
+    fHistPtFake[iType]->Sumw2();
+    fOutput->Add(fHistPtFake[iType]);
+
+    fHistEtaPhi[iType] = new TH2F(Form("hEtaPhi%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
+    fHistEtaPhi[iType]->Sumw2();
+    fOutput->Add(fHistEtaPhi[iType]);
+    fHistEtaPhiGood[iType] = new TH2F(Form("hEtaPhiGood%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
+    fHistEtaPhiGood[iType]->Sumw2();
+    fOutput->Add(fHistEtaPhiGood[iType]);
+    fHistEtaPhiFake[iType] = new TH2F(Form("hEtaPhiFake%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
+    fHistEtaPhiFake[iType]->Sumw2();
+    fOutput->Add(fHistEtaPhiFake[iType]);
+
+    fHistChi2[iType]=new TH1F(Form("hChi2%s",tit[iType].Data()),"",100,0.,10.);
+    fHistChi2[iType]->Sumw2();
+    fOutput->Add(fHistChi2[iType]);
+    fHistChi2Good[iType]=new TH1F(Form("hChi2Good%s",tit[iType].Data()),"",100,0.,10.);
+    fHistChi2Good[iType]->Sumw2();
+    fOutput->Add(fHistChi2Good[iType]);
+    fHistChi2Fake[iType]=new TH1F(Form("hChi2Fake%s",tit[iType].Data()),"",100,0.,10.);
+    fHistChi2Fake[iType]->Sumw2();
+    fOutput->Add(fHistChi2Fake[iType]);
+
+    fHistNclu[iType]=new TH1F(Form("hNclu%s",tit[iType].Data()),"",7,-0.5,6.5);
+    fHistNclu[iType]->Sumw2();
+    fOutput->Add(fHistNclu[iType]);
+    fHistNcluGood[iType]=new TH1F(Form("hNcluGood%s",tit[iType].Data()),"",7,-0.5,6.5);
+    fHistNcluGood[iType]->Sumw2();
+    fOutput->Add(fHistNcluGood[iType]);
+    fHistNcluFake[iType]=new TH1F(Form("hNcluFake%s",tit[iType].Data()),"",7,-0.5,6.5);
+    fHistNcluFake[iType]->Sumw2();
+    fOutput->Add(fHistNcluFake[iType]);
+
+    fHistdedxvsPt3cls[iType] = new TH2F(Form("hdedxvsPt3cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
+    fHistdedxvsPt3cls[iType]->Sumw2();
+    fOutput->Add(fHistdedxvsPt3cls[iType]);
+
+    fHistdedxvsPt4cls[iType] = new TH2F(Form("hdedxvsPt4cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
+    fHistdedxvsPt4cls[iType]->Sumw2();
+    fOutput->Add(fHistdedxvsPt4cls[iType]);
+  }
+
+  //-----------------------------------------------------------
+
   TString spname[3]={"Pion","Kaon","Proton"};
   TString hisname;
   const Int_t nbins = fNPtBins;
   Double_t xbins[nbins+1];
   for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
 
-
-  fHistPtTPCITSAll = new TH1F("hPtTPCITSAll","",100,0.,2.);
-  fHistPtTPCITSAll->Sumw2();
-  fOutput->Add(fHistPtTPCITSAll);
-  fHistPtTPCITSGood = new TH1F("hPtTPCITSGood","",100,0.,2.);
-  fHistPtTPCITSGood->Sumw2();
-  fOutput->Add(fHistPtTPCITSGood);
-  fHistPtTPCITSFake = new TH1F("hPtTPCITSFake","",100,0.,2.);
-  fHistPtTPCITSFake->Sumw2();
-  fOutput->Add(fHistPtTPCITSFake);
-
-  fHistPtITSsaAll  = new TH1F("hPtITSsaAll","",100,0.,2.);
-  fHistPtITSsaAll->Sumw2();
-  fOutput->Add(fHistPtITSsaAll);
-  fHistPtITSsaGood  = new TH1F("hPtITSsaGood","",100,0.,2.);
-  fHistPtITSsaGood->Sumw2();
-  fOutput->Add(fHistPtITSsaGood);
-  fHistPtITSsaFake  = new TH1F("hPtITSsaFake","",100,0.,2.);
-  fHistPtITSsaFake->Sumw2();
-  fOutput->Add(fHistPtITSsaFake);
-
-  fHistPtITSpureSAAll = new TH1F("hPtITSpureSAAll","",100,0.,2.);
-  fHistPtITSpureSAAll->Sumw2();
-  fOutput->Add(fHistPtITSpureSAAll);
-  fHistPtITSpureSAGood = new TH1F("hPtITSpureSAGood","",100,0.,2.);
-  fHistPtITSpureSAGood->Sumw2();
-  fOutput->Add(fHistPtITSpureSAGood);
-  fHistPtITSpureSAFake = new TH1F("hPtITSpureSAFake","",100,0.,2.);
-  fHistPtITSpureSAFake->Sumw2();
-  fOutput->Add(fHistPtITSpureSAFake);
-
   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
 
     hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
@@ -307,6 +362,15 @@ void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
     
   }
 
+  fHistMCPhiResid=new TH2F("hMCPhiResid","",nbins,xbins,100,-0.5,0.5);
+  fHistMCPhiResid->Sumw2();
+  fOutput->Add(fHistMCPhiResid);
+  fHistPhiResid=new TH2F("hPhiResid","",nbins,xbins,100,-0.5,0.5);
+  fHistPhiResid->Sumw2();
+  fOutput->Add(fHistPhiResid);
+
+  PostData(1,fOutput);
+
 }
 //______________________________________________________________________________
 void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
@@ -347,26 +411,104 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
     }
   }
 
-  PostData(1, fOutput);
 
   fHistNEvents->Fill(0);
   const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
   if(spdv->GetNContributors()<=0) return;
   fHistNEvents->Fill(1);
-
+  
+  const Int_t ntSize=36;
+  Float_t xnt[ntSize];
+  
   Int_t ntracks = esd->GetNumberOfTracks();
   for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
     AliESDtrack * track = esd->GetTrack(iTrack);
     if (!track) continue;
+    Float_t pttrack=track->Pt();
+    Float_t ptrack=track->P();
+    Float_t pxtrack=track->Px();
+    Float_t pytrack=track->Py();
+    Float_t pztrack=track->Pz();
+    Float_t impactXY=-999, impactZ=-999;
+    track->GetImpactParameters(impactXY, impactZ);
+    Double_t dedxlay[4];
+    track->GetITSdEdxSamples(dedxlay);
+    Float_t dedx=track->GetITSsignal();
+    Float_t chi2=track->GetITSchi2();
     Int_t status=track->GetStatus();
-    if(!(status&AliESDtrack::kITSrefit)) continue;
+    Bool_t isITSrefit=kFALSE;
+    Bool_t isTPCin=kFALSE;
+    Bool_t isTPCrefit=kFALSE;
+    Bool_t isPureSA=kFALSE;
+    if(status&AliESDtrack::kITSrefit) isITSrefit=kTRUE; 
+    if(status&AliESDtrack::kTPCin) isTPCin=kTRUE; 
+    if(status&AliESDtrack::kTPCrefit) isTPCrefit=kTRUE; 
+    if(status&AliESDtrack::kITSpureSA) isPureSA=kTRUE;
     Bool_t isSA=kTRUE;
     if(status&AliESDtrack::kTPCin) isSA=kFALSE;
     Int_t nTPCclus=track->GetNcls(1);
     Int_t nITSclus=track->GetNcls(0);
+    UChar_t clumap=track->GetITSClusterMap();
+    Int_t idet;
+    Float_t xloc,zloc;
+    Int_t statusLay[6];
+    for(Int_t iLay=0; iLay<6;iLay++) track->GetITSModuleIndexInfo(iLay,idet,statusLay[iLay],xloc,zloc);
+    Int_t trlabel=track->GetLabel();
+    Float_t ptgen=-999.;
+    Float_t pgen=-999.;
+    Float_t pxgen=-999.;
+    Float_t pygen=-999.;
+    Float_t pzgen=-999.;
+    Float_t etagen=-999.;
+    Float_t phigen=-999.;
+    if(fReadMC){
+      TParticle* part = stack->Particle(TMath::Abs(trlabel));
+      ptgen=part->Pt();
+      pgen=part->P();
+      pxgen=part->Px();
+      pygen=part->Py();
+      pzgen=part->Pz();
+      etagen=part->Eta();
+      phigen=part->Phi();
+    }
+
+    Int_t indexn=0;
+    xnt[indexn++]=pttrack;
+    xnt[indexn++]=pxtrack;
+    xnt[indexn++]=pytrack;
+    xnt[indexn++]=pztrack;
+    xnt[indexn++]=track->Eta();
+    xnt[indexn++]=track->Phi();
+    xnt[indexn++]=impactXY;
+    xnt[indexn++]=impactZ;
+    for(Int_t iLay=0; iLay<4; iLay++) xnt[indexn++]=dedxlay[iLay];
+    xnt[indexn++]=dedx;
+    xnt[indexn++]=chi2;
+    xnt[indexn++]=status;
+    xnt[indexn++]=isITSrefit;
+    xnt[indexn++]=isTPCin;
+    xnt[indexn++]=isTPCrefit;
+    xnt[indexn++]=isPureSA;
+    xnt[indexn++]=nITSclus;
+    xnt[indexn++]=nTPCclus;
+    xnt[indexn++]=clumap;
+    for(Int_t iLay=0; iLay<6; iLay++) xnt[indexn++]=(Float_t)statusLay[iLay];
+    xnt[indexn++]=trlabel;
+    xnt[indexn++]=ptgen;
+    xnt[indexn++]=pxgen;
+    xnt[indexn++]=pygen;
+    xnt[indexn++]=pzgen;
+    xnt[indexn++]=etagen;
+    xnt[indexn++]=phigen;
+    xnt[indexn++]=ntracks;
+
+    if(indexn>ntSize) printf("AliAnalysisTaskITSsaTracks: ERROR ntuple insexout of range\n");
+    if(fFillNtuple) fNtupleTracks->Fill(xnt);
+
+    if(!(status&AliESDtrack::kITSrefit)) continue;
     if(nITSclus<fMinITSpts)continue;
+    if((chi2/nITSclus) > fMaxITSChi2Clu) continue;
 
-    UChar_t clumap=track->GetITSClusterMap();
     Int_t nSPDPoints=0;
     for(Int_t i=0; i<2; i++){
       if(clumap&(1<<i)) ++nSPDPoints;
@@ -377,52 +519,50 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
     for(Int_t i=2; i<6; i++){
       if(clumap&(1<<i)) ++nPointsForPid;
     }
-    
-    Float_t pttrack=track->Pt();
-    Float_t ptrack=track->P();
-    Float_t dedx=track->GetITSsignal();
-    Int_t hadronSpecie=-1;
-    if(status&AliESDtrack::kTPCin){ 
-      fHistPtTPCITSAll->Fill(pttrack);
+
+    Int_t iTrackType=-1;
+    if(status&AliESDtrack::kTPCin){
+      iTrackType=kTypeTPCITS;
     }else{
-      if(status&AliESDtrack::kITSpureSA){
-       fHistPtITSpureSAAll->Fill(pttrack);
+      if(status&AliESDtrack::kITSpureSA) iTrackType=kTypeITSpureSA;
+      else iTrackType=kTypeITSsa;
+    }
+    if(iTrackType<0 || iTrackType>=kNtrackTypes) continue;
+    fHistPt[iTrackType]->Fill(pttrack);
+    fHistEtaPhi[iTrackType]->Fill(track->Eta(),track->Phi());
+    fHistChi2[iTrackType]->Fill(chi2/nITSclus);
+    fHistNclu[iTrackType]->Fill(nITSclus);
+    if(nPointsForPid==3){
+      fHistdedxvsPt3cls[iTrackType]->Fill(pttrack,dedx);
+    }
+    if(nPointsForPid==4){
+      fHistdedxvsPt4cls[iTrackType]->Fill(pttrack,dedx);
+    }
+    if(fReadMC){
+      if(trlabel<0){
+       fHistPtFake[iTrackType]->Fill(pttrack);
+       fHistEtaPhiFake[iTrackType]->Fill(track->Eta(),track->Phi());
+       fHistChi2Fake[iTrackType]->Fill(chi2/nITSclus);
+       fHistNcluFake[iTrackType]->Fill(nITSclus);
       }else{
-       fHistPtITSsaAll->Fill(pttrack);
+       fHistPtGood[iTrackType]->Fill(pttrack);
+       fHistEtaPhiGood[iTrackType]->Fill(track->Eta(),track->Phi());
+       fHistChi2Good[iTrackType]->Fill(chi2/nITSclus);
+       fHistNcluGood[iTrackType]->Fill(nITSclus);
       }
     }
 
+    Int_t hadronSpecie=-1;
     if(fReadMC && fUseMCId){
-      Int_t trlabel=track->GetLabel();
       if(trlabel>=0){
-       if(status&AliESDtrack::kTPCin){ 
-         fHistPtTPCITSGood->Fill(pttrack);
-       }else{
-         if(status&AliESDtrack::kITSpureSA){
-           fHistPtITSpureSAGood->Fill(pttrack);
-         }else{
-           fHistPtITSsaGood->Fill(pttrack);
-         }
-       }
-      }else{
-       if(status&AliESDtrack::kTPCin){ 
-         fHistPtTPCITSFake->Fill(pttrack);
-       }else{
-         if(status&AliESDtrack::kITSpureSA){
-           fHistPtITSpureSAFake->Fill(pttrack);
-         }else{
-           fHistPtITSsaFake->Fill(pttrack);
-         }
-       }
-       continue; // fake track
+       TParticle* part = stack->Particle(trlabel);
+       Int_t pdg=TMath::Abs(part->GetPdgCode());
+       if(pdg==211) hadronSpecie=kPion;
+       else if(pdg==321) hadronSpecie=kKaon;
+       else if(pdg==2212) hadronSpecie=kProton;
       }
-      TParticle* part = stack->Particle(trlabel);
-      Int_t pdg=TMath::Abs(part->GetPdgCode());
-      if(pdg==211) hadronSpecie=kPion;
-      else if(pdg==321) hadronSpecie=kKaon;
-      else if(pdg==2212) hadronSpecie=kProton;
     }else{
-      if(nPointsForPid<fMinITSptsForPid)continue;
+      if(nPointsForPid<fMinPtsforPid)continue;
       AliITSPIDResponse pidits(fReadMC);
       Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
       if(nSigmaPion>-2. && nSigmaPion<2.){
@@ -440,7 +580,7 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
       } 
     }
     if(hadronSpecie<0) continue;
-    if(status&AliESDtrack::kTPCin){ // TPC+ITS tracks
+    if(iTrackType==kTypeTPCITS){ // TPC+ITS tracks
       fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
       fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
       fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
@@ -448,39 +588,27 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
       for(Int_t iBit=0; iBit<6; iBit++){
        if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
       }
-    }else{ // ITS standalone and pureSA tracks
-      if(status&AliESDtrack::kITSpureSA){
-       Float_t impactXY=-999, impactZ=-999;
-       track->GetImpactParameters(impactXY, impactZ);
-       fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
-       fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
-       fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
-       fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
-       fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
-       if(nPointsForPid==3){
-         fHistdedxvsPtITSpureSA3cls->Fill(pttrack,dedx);
-         fHistdedxvsPITSpureSA3cls->Fill(ptrack,dedx);
-       }
-       if(nPointsForPid==4){
-         fHistdedxvsPtITSpureSA4cls->Fill(pttrack,dedx);
-         fHistdedxvsPITSpureSA4cls->Fill(ptrack,dedx);
+    }else if(iTrackType==kTypeITSpureSA){ // TPC+ITS tracks
+      fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
+      fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
+      fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
+      fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
+      fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
+      fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
+      Int_t outerLay=-1;
+      for(Int_t iBit=0; iBit<6; iBit++){
+       if(clumap&(1<<iBit)){
+         fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
+         if(iBit>outerLay) outerLay=iBit;
        }
-       fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
-       Int_t outerLay=-1;
-       for(Int_t iBit=0; iBit<6; iBit++){
-         if(clumap&(1<<iBit)){
-           fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
-           if(iBit>outerLay) outerLay=iBit;
-         }
        }
-       fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);   
+      fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);    
+      
        
-       
-       if(fReadMC){  
-         Int_t trlabel=track->GetLabel();
-         if(trlabel<0)continue; // fake track
+      if(fReadMC){  
+       if(trlabel>=0){
          TParticle* part = stack->Particle(trlabel);
-         Float_t ptgen=part->Pt();
+         ptgen=part->Pt();
          Float_t invpttrack=track->OneOverPt();
          Float_t invptgen=0.;
          if(ptgen>0.) invptgen=1./ptgen;
@@ -488,19 +616,24 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
          fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
          fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
          fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);       
+         Float_t deltaphi=track->Phi()-part->Phi();
+         if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
+         if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
+         fHistMCPhiResid->Fill(pttrack,deltaphi);
        }
-      }else{
-       fHistPtITSsa[hadronSpecie]->Fill(pttrack);
-       fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
-       fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
-       fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
-       for(Int_t iBit=0; iBit<6; iBit++){
-         if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
-       }
+      }
+    }else if(iTrackType==kTypeITSsa){ // TPC+ITS tracks
+      fHistPtITSsa[hadronSpecie]->Fill(pttrack);
+      fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
+      fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
+      fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
+      for(Int_t iBit=0; iBit<6; iBit++){
+       if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
       }
     }
-
-    if(nITSclus<fMinITSptsForMatch || nTPCclus<fMinTPCpts) continue;      
+  
+    // match TPCITS with ITSpureSA
+    if(nITSclus<fMinITSpts || nTPCclus<fMinTPCpts) continue;      
     Int_t idxMI[12],idxSA[12];
     for(Int_t icl=0; icl<12; icl++){ 
       idxMI[icl]=-1; 
@@ -587,7 +720,11 @@ void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
        fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
        fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
        fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
-       fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);   
+       fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
+       Float_t deltaphi=track->Phi()-track2->Phi();
+       if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
+       if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
+       fHistPhiResid->Fill(pt1,deltaphi);
       }
     }
   }
@@ -604,25 +741,10 @@ void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
     printf("ERROR: fOutput not available\n");
     return;
   }
+  fNtupleTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("ntupleTracks"));
   fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
-  fHistdedxvsPtITSpureSA3cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPtITSpureSA3cls"));
-  fHistdedxvsPITSpureSA3cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPITSpureSA3cls"));
-  fHistdedxvsPtITSpureSA4cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPtITSpureSA4cls"));
-  fHistdedxvsPITSpureSA4cls= dynamic_cast<TH2F*>(fOutput->FindObject("hdedxvsPITSpureSA4cls"));
-    
-  TString spname[3]={"Pion","Kaon","Proton"};
-
-  fHistPtTPCITSAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSAll"));
-  fHistPtTPCITSGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSGood"));
-  fHistPtTPCITSFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtTPCITSFake"));
 
-  fHistPtITSsaAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaAll"));
-  fHistPtITSsaGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaGood"));
-  fHistPtITSsaFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSsaFake"));
-
-  fHistPtITSpureSAAll =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAAll"));
-  fHistPtITSpureSAGood =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAGood"));
-  fHistPtITSpureSAFake =dynamic_cast<TH1F*>(fOutput->FindObject("hPtITSpureSAFake"));
+  TString spname[3]={"Pion","Kaon","Proton"};
 
   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
     fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
@@ -656,6 +778,8 @@ void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
     fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
   
   }
+  fHistMCPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hMCPhiResid"));
+  fHistPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hPhiResid"));
   return;
 }
 
@@ -663,4 +787,3 @@ void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
 
 
 
-

diff --git a/PWG1/ITS/AliAnalysisTaskITSsaTracks.h b/PWG1/ITS/AliAnalysisTaskITSsaTracks.h
index 06b13c63a6fa0c513950f6d5f4e99a2263f0a9ec..7e7607c67490ce2568384f316e3801d1147bb81c 100644 (file)
--- a/PWG1/ITS/AliAnalysisTaskITSsaTracks.h
+++ b/PWG1/ITS/AliAnalysisTaskITSsaTracks.h
@@ -15,6 +15,7 @@
 //*************************************************************************
 
 class TList;
+class TNtuple;
 class TH1F;
 class TH2F;
 class TTree;
@@ -45,10 +46,10 @@ class AliAnalysisTaskITSsaTracks : public AliAnalysisTaskSE {
     fMinSPDpts=minp;
   }
   void SetMinPointsForITSPid(Int_t minp=3){
-    fMinITSptsForPid=minp;
+    fMinPtsforPid=minp;
   }
-  void SetMinPointsForMatchToTPC(Int_t minp=6){
-    fMinITSptsForMatch=minp;
+  void SetITChi2Cut(Float_t maxchi2=2.5){
+    fMaxITSChi2Clu=maxchi2;
   }
 
   void SetPtBins(Int_t n, Double_t* lim);
@@ -64,6 +65,9 @@ class AliAnalysisTaskITSsaTracks : public AliAnalysisTaskSE {
     fRequirePoint[1]=kTRUE;
   }
 
+  void SetFillNtuple(Bool_t fill=kTRUE){
+    fFillNtuple=fill;
+  }  
   void SetReadMC(Bool_t optMC=kTRUE){
     fReadMC=optMC;
   }
@@ -74,6 +78,7 @@ class AliAnalysisTaskITSsaTracks : public AliAnalysisTaskSE {
 
  private:
   enum {kPion=0,kKaon,kProton,kNspecies};
+  enum {kTypeTPCITS=0, kTypeITSsa, kTypeITSpureSA, kNtrackTypes};
   enum {kMaxPtBins=40};
 
   AliAnalysisTaskITSsaTracks(const AliAnalysisTaskITSsaTracks &source);
@@ -82,29 +87,35 @@ class AliAnalysisTaskITSsaTracks : public AliAnalysisTaskSE {
   TList*  fOutput;          //! list of output histos
   TH1F*   fHistNEvents;     //! histo with N of events  
 
+  
+  TH1F*   fHistPt[kNtrackTypes];          //! pt distr., no PID
+  TH1F*   fHistPtGood[kNtrackTypes];      //! pt distr. good tracks, no PID
+  TH1F*   fHistPtFake[kNtrackTypes];      //! pt distr. fake tracks, no PID
+
+  TH2F*   fHistEtaPhi[kNtrackTypes];      //! etaphi distr., no PID
+  TH2F*   fHistEtaPhiGood[kNtrackTypes];  //! etaphi distr. good tracks, no PID
+  TH2F*   fHistEtaPhiFake[kNtrackTypes];  //! etaphi distr. fake tracks, no PID
+
+  TH1F*   fHistChi2[kNtrackTypes];        //! chi2 distr., no PID
+  TH1F*   fHistChi2Good[kNtrackTypes];    //! chi2 distr., good tracks, no PID
+  TH1F*   fHistChi2Fake[kNtrackTypes];    //! chi2 distr., fake tracks, no PID
+
+  TH1F*   fHistNclu[kNtrackTypes];        //! ITS clu distr., no PID
+  TH1F*   fHistNcluGood[kNtrackTypes];    //! ITS clu distr., good tracks, no PID
+  TH1F*   fHistNcluFake[kNtrackTypes];    //! ITS clu distr., fake tracks, no PID
+
+  TH2F*   fHistdedxvsPt3cls[kNtrackTypes]; //! dedx vs. pt for tracks with 3 clus in SDD+SSD
+  TH2F*   fHistdedxvsPt4cls[kNtrackTypes]; //! dedx vs. pt for tracks with 4 clus in SDD+SSD
+
+
   TH1F*   fHistPtTPCITS[kNspecies];    //! pt distribution of TPC+ITS tracks
   TH1F*   fHistPtITSsa[kNspecies];     //! pt distribution of ITSsa tracks
   TH1F*   fHistPtITSpureSA[kNspecies]; //! pt distribution of ITS pure SA tracks
 
-  TH1F*   fHistPtTPCITSAll;     //! pt distribution of all TPC+ITS tracks 
-  TH1F*   fHistPtTPCITSGood;    //! pt distribution of good TPC+ITS tracks 
-  TH1F*   fHistPtTPCITSFake;    //! pt distribution of fake TPC+ITS tracks
-  TH1F*   fHistPtITSsaAll;      //! pt distribution of all ITSsa tracks
-  TH1F*   fHistPtITSsaGood;     //! pt distribution of good ITSsa tracks
-  TH1F*   fHistPtITSsaFake;     //! pt distribution of fake ITSsa tracks
-  TH1F*   fHistPtITSpureSAAll;  //! pt distribution of all ITS pure SA tracks
-  TH1F*   fHistPtITSpureSAGood; //! pt distribution of good ITS pure SA tracks
-  TH1F*   fHistPtITSpureSAFake; //! pt distribution of fake ITS pure SA tracks
-
   TH2F*   fHistEtaPhiTPCITS[kNspecies];    //! etaphi distr. of TPC+ITS tracks
   TH2F*   fHistEtaPhiITSsa[kNspecies];     //! etaphi distr. of ITSsa tracks
   TH2F*   fHistEtaPhiITSpureSA[kNspecies]; //! etaphi distr. of ITSpureSA tracks
 
-  TH2F*   fHistdedxvsPtITSpureSA3cls; //! dedx for ITSpureSA tracks vs. pt
-  TH2F*   fHistdedxvsPITSpureSA3cls; //! dedx for ITSpureSA tracks vs. p
-  TH2F*   fHistdedxvsPtITSpureSA4cls; //! dedx for ITSpureSA tracks vs. pt
-  TH2F*   fHistdedxvsPITSpureSA4cls; //! dedx for ITSpureSA tracks vs. p
-
   TH2F*   fHistNcluTPCITS[kNspecies];    //! n. of clusters for TPC+ITS tracks vs. pt
   TH2F*   fHistNcluITSsa[kNspecies];     //! n. of clusters for ITSsa tracks vs. pt
   TH2F*   fHistNcluITSpureSA[kNspecies]; //! n. of clusters for ITSpureSA tracks vs. pt
@@ -125,20 +136,24 @@ class AliAnalysisTaskITSsaTracks : public AliAnalysisTaskSE {
   TH2F*   fHistMCInvPtResid[kNspecies];    //! 1/pt residuals (MC) vs. pt
   TH2F*   fHistMCInvPtRelResid[kNspecies]; //! 1/pt relative residulas (MC) vs. pt
 
-  Int_t fNPtBins;                  // number of Pt bins
+  TH2F*   fHistMCPhiResid; //! phi residuals in pt bins
+  TH2F*   fHistPhiResid;   //! phi residuals in pt bins
+  TNtuple* fNtupleTracks;  //! output ntuple
+
+  Int_t   fNPtBins;                  // number of Pt bins
   Float_t fPtLimits[kMaxPtBins+1]; // Pt bin limits
-  Int_t   fMinITSpts;          // Minimum number of ITS points per track
-  Int_t   fMinSPDpts;          // Minimum number of SPD points per track
-  Int_t   fMinITSptsForPid;    // Minimum number of SDD+SSD points per track
-  Int_t   fMinITSptsForMatch;  // Minimum number of ITS points to macth to TPC
-  Int_t   fMinTPCpts;          // Minimum number of TPC points per track
+  Int_t   fMinITSpts;       // Minimum number of ITS points per track
+  Int_t   fMinSPDpts;       // Minimum number of SPD points per track
+  Int_t   fMinPtsforPid;    // Minimum number of SDD+SSD points per track
+  Int_t   fMinTPCpts;       // Minimum number of TPC points per track
+  Float_t fMaxITSChi2Clu;   // Maximum value of ITS chi2 per cluster
   Bool_t  fRequirePoint[6]; // require point in given layer
+  Bool_t  fFillNtuple;      // flag to control fill of ntuple  
   Bool_t  fReadMC;          // flag read/not-read MC truth info
   Bool_t  fUseMCId;         // flag use/not-use MC identity for PID
 
-  ClassDef(AliAnalysisTaskITSsaTracks,1);  
+  ClassDef(AliAnalysisTaskITSsaTracks,2);  
 };
 
 
 #endif
-