Updates in ITSsa pi,kp task: new multiplicity estimator, storeage of DCA cut, code...
authorfprino <prino@to.infn.it>
Wed, 6 Aug 2014 22:31:34 +0000 (00:31 +0200)
committerfprino <prino@to.infn.it>
Wed, 6 Aug 2014 22:59:23 +0000 (00:59 +0200)
PWGLF/CMakelibPWGLFspectra.pkg
PWGLF/PWGLFspectraLinkDef.h
PWGLF/SPECTRA/PiKaPr/ITSsa/AddTaskITSsaSpectra.C
PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx
PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.h
PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C [new file with mode: 0644]

index e83dea4..93d0399 100644 (file)
@@ -44,6 +44,7 @@ set ( SRCS SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAcceptanceCuts.cxx
        SPECTRA/ChargedHadrons/dNdPt/AlidNdPtTrackDumpTask.cxx
        SPECTRA/ChargedHadrons/dNdPt/AliPtResolAnalysis.cxx
        SPECTRA/ChargedHadrons/dNdPt/AliPtResolAnalysisPbPb.cxx
+       SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx
        SPECTRA/PiKaPr/TPCTOF/AliAnalysisCombinedHadronSpectra.cxx
        SPECTRA/PiKaPr/TPCTOFpA/AliAnalysisTPCTOFpA.cxx
        SPECTRA/PiKaPr/TOF/pp7/TOFSpectrappAnalysis.cxx
index b65c165..48a5235 100644 (file)
@@ -18,6 +18,7 @@
 #pragma link C++ class AlidNdPtTrackDumpTask+;
 #pragma link C++ class AliPtResolAnalysis+;
 #pragma link C++ class AliPtResolAnalysisPbPb+;
+#pragma link C++ class AliAnalysisTaskSEITSsaSpectra+;
 #pragma link C++ class AliAnalysisCombinedHadronSpectra+;
 #pragma link C++ class AliAnalysisTPCTOFpA+;
 #pragma link C++ class TOFSpectrappAnalysis+;
index 15eb068..4d09259 100644 (file)
@@ -55,7 +55,7 @@ AliAnalysisTaskSEITSsaSpectra *AddTaskITSsaSpectra(Int_t optNtuple=0,Int_t readM
   outputFileName += ":PWG2SpectraITSsa";
   
   AliAnalysisDataContainer *coutput =0x0;
-  
+  AliAnalysisDataContainer *coutputCuts = 0x0;
   if(hi)
     {
       coutput = mgr->CreateContainer(Form("clistITSsaCent%ito%i",lowcut,upcut),
@@ -67,11 +67,18 @@ AliAnalysisTaskSEITSsaSpectra *AddTaskITSsaSpectra(Int_t optNtuple=0,Int_t readM
     {
       coutput = mgr->CreateContainer(Form("clistITSsaMult%ito%i",lowcut,upcut),
                                     TList::Class(),
-                                     AliAnalysisManager::kOutputContainer,
-                                     outputFileName );    
+                                    AliAnalysisManager::kOutputContainer,
+                                    outputFileName );    
     }
+
+  coutputCuts = mgr->CreateContainer(Form("DCACutMult%ito%i",lowcut,upcut),
+                                    TList::Class(),
+                                    AliAnalysisManager::kParamContainer,
+                                    outputFileName );
+  
   
   mgr->ConnectInput(taskits, 0, mgr->GetCommonInputContainer());
   mgr->ConnectOutput(taskits, 1, coutput);
+  mgr->ConnectOutput(taskits, 2, coutputCuts);
   return taskits;
 }   
index def4a83..8b3ef29 100644 (file)
@@ -25,6 +25,7 @@
 ///////////////////////////////////////////////////////////////////////////
 
 #include <TH1F.h>
+#include <TF1.h>
 #include <TRandom3.h>
 #include <TH2F.h>
 #include <TChain.h>
@@ -53,10 +54,10 @@ ClassImp(AliAnalysisTaskSEITSsaSpectra)
 
 //________________________________________________________________________
 AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
-AliAnalysisTaskSE("Task CFit"),
+AliAnalysisTaskSE("taskITSsaSpectra"),
   fESD(0),
-  fesdTrackCutsMult(0),
   fOutput(0),
+  fListCuts(0),
   fHistNEvents(0),
   fHistMult(0),
   fHistCen(0),
@@ -67,6 +68,8 @@ AliAnalysisTaskSE("Task CFit"),
   fHistDEDXdouble(0),
   fHistBeforeEvSel(0),
   fHistAfterEvSel(0),
+  fDCAxyCutFunc(0x0),
+  fDCAzCutFunc(0x0),
   fITSPIDResponse(0),
   fMinSPDPts(1),
   fMinNdEdxSamples(3),
@@ -81,7 +84,7 @@ AliAnalysisTaskSE("Task CFit"),
   fUpMult(-1),
   fLowCentrality(-1.0),
   fUpCentrality(-1.0),
-  fSPD(0),
+  fMultEstimator(0),
   fHImode(0),
   fYear(2010),
   fMC(kFALSE), 
@@ -98,25 +101,9 @@ AliAnalysisTaskSE("Task CFit"),
   Double_t xbins[kNbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
   for(Int_t iBin=0; iBin<kNbins+1; iBin++) fPtBinLimits[iBin]=xbins[iBin];
   fRandGener=new TRandom3(0);
-  fesdTrackCutsMult = new AliESDtrackCuts;
-   // TPC  
-  fesdTrackCutsMult->SetMinNClustersTPC(70);
-  fesdTrackCutsMult->SetMaxChi2PerClusterTPC(4);
-  fesdTrackCutsMult->SetAcceptKinkDaughters(kFALSE);
-  fesdTrackCutsMult->SetRequireTPCRefit(kTRUE);
-  // ITS
-  fesdTrackCutsMult->SetRequireITSRefit(kTRUE);
-  fesdTrackCutsMult->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
-                                            AliESDtrackCuts::kAny);
-  fesdTrackCutsMult->SetDCAToVertex2D(kFALSE);
-  fesdTrackCutsMult->SetRequireSigmaToVertex(kFALSE);
-  fesdTrackCutsMult->SetEtaRange(-0.8,+0.8);
-  fesdTrackCutsMult->SetPtRange(0.15, 1e10);                                        
-  SetYear(2010);
-
   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
+  DefineOutput(2,TList::Class());
   Printf("end of AliAnalysisTaskSEITSsaSpectra");
 }
 
@@ -127,8 +114,14 @@ AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
     delete fOutput;
     fOutput = 0;
   }
+  if (fListCuts) {
+    delete fListCuts;
+    fListCuts = 0;
+  }
   if(fRandGener) delete fRandGener;
   if(fITSPIDResponse) delete fITSPIDResponse;
+  delete fDCAxyCutFunc;
+  delete fDCAzCutFunc;
 }
 
 //________________________________________________________________________
@@ -170,26 +163,37 @@ Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s) const {
 
 
 //________________________________________________________________________
-void AliAnalysisTaskSEITSsaSpectra::SetYear(Int_t year){
-  // Set year dependent quantities
-  fYear=year;
-  if(fYear==2009){
-    fesdTrackCutsMult->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/TMath::Power(pt,0.9)");  //2009 standard cut
-    fesdTrackCutsMult->SetMaxDCAToVertexZ(20); //2009 standard cut
-  }else{
-    fesdTrackCutsMult->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); //2010 standard cut
-    fesdTrackCutsMult->SetMaxDCAToVertexZ(2); //2010 standard cut
-  }
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC) const {
+Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const {
   // cut on transverse impact parameter updaated on 20-5-2010
   // from the study of L. Milano, F. Prino on the ITS standalone tracks
   // using the common binning of the TPC tracks
+
+  Double_t xyMax = fDCAxyCutFunc->Eval(pt); //in micron
+  if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;  
+  Double_t zMax = fDCAzCutFunc->Eval(pt); //in micron
+  if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
+  
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const {
+  // convert eta to y
+  Double_t mt = TMath::Sqrt(m*m + pt*pt);
+  return TMath::ASinH(pt/mt*TMath::SinH(eta));
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskSEITSsaSpectra::Init(){
+  //
+  // Initialization
+  //
+  fListCuts=new TList();
+  fListCuts->SetOwner();
   Double_t xyP[3];
   Double_t zP[3];
-  if(optMC){
+  if(fMC){
     if(fYear==2009){
       xyP[0]=88.63; //MC LHC10a12
       xyP[1]=19.57;
@@ -205,8 +209,7 @@ Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ
       zP[1]=59.8;
       zP[2]=1.2;
     }
-  }
-  else{
+  }else{
     if(fYear==2009){
       xyP[0]=85.28;//DATA 900 GeV pass6
       xyP[1]=25.78;
@@ -223,24 +226,21 @@ Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ
       zP[2]=1.2;
     }
   }
-  Double_t xySigma = xyP[0] + xyP[1]/TMath::Power(TMath::Abs(pt),xyP[2]);
-  Double_t xyMax = fNSigmaDCAxy*xySigma; //in micron
-  if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;
-  
-  Double_t zSigma = zP[0] + zP[1]/TMath::Power(TMath::Abs(pt),zP[2]);
-  Double_t zMax = fNSigmaDCAz*zSigma; //in micron
-  if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
-  
-  return kTRUE;
-}
+  fDCAxyCutFunc = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
+  for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutFunc->SetParameter(ipar,xyP[ipar]);
+  fDCAxyCutFunc->SetParameter(3,fNSigmaDCAxy);
 
-//________________________________________________________________________
-Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const {
-  // convert eta to y
-  Double_t mt = TMath::Sqrt(m*m + pt*pt);
-  return TMath::ASinH(pt/mt*TMath::SinH(eta));
-}
 
+                                               
+  fDCAzCutFunc = new TF1("fDCAzCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
+  for(Int_t ipar=0; ipar<3; ipar++) fDCAzCutFunc->SetParameter(ipar,zP[ipar]);
+  fDCAzCutFunc->SetParameter(3,fNSigmaDCAz);
+
+  fListCuts->Add(fDCAxyCutFunc);
+  fListCuts->Add(fDCAzCutFunc);
+
+  PostData(2,fListCuts);
+}
 
 //________________________________________________________________________
 void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
@@ -251,14 +251,26 @@ void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
   fOutput->SetOwner();
   fOutput->SetName("Spiderman");
   
-  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",8,-1.5,6.5);
+  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",9,-1.5,7.5);
   fHistNEvents->Sumw2();
   fHistNEvents->SetMinimum(0);
+  fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from ESD");
+  fHistNEvents->GetXaxis()->SetBinLabel(2,"Pass Phys. Sel. + Trig");
+  fHistNEvents->GetXaxis()->SetBinLabel(3,"SDD read out");
+  fHistNEvents->GetXaxis()->SetBinLabel(4,"In mult. range");
+  fHistNEvents->GetXaxis()->SetBinLabel(5,"With SPD vertex");
+  fHistNEvents->GetXaxis()->SetBinLabel(6,"Vertex contributors >0");
+  fHistNEvents->GetXaxis()->SetBinLabel(7,"|zVertex|<10");
+  fHistNEvents->GetXaxis()->SetBinLabel(8,"Error on zVertex<0.5");
+  fHistNEvents->GetXaxis()->SetBinLabel(9,"Good Z vertex");
   fOutput->Add(fHistNEvents);
   
   fHistMult = new TH1F("fHistMult", "Event Multiplicity",3000,-0.5,2999.5);
   fHistMult->Sumw2();
   fHistMult->SetMinimum(0);
+  if(fMultEstimator==0) fHistMult->GetXaxis()->SetTitle("Multiplicity |#eta|<0.8");
+  else if(fMultEstimator==1) fHistMult->GetXaxis()->SetTitle("Tracklets |#eta|<0.8");
+  else if(fMultEstimator==2) fHistMult->GetXaxis()->SetTitle("Clusters on SPD1");
   fOutput->Add(fHistMult);
   
   fHistCen = new TH1F("fHistCen", "Event Centrality",101,-0.5,100.5);
@@ -542,6 +554,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
   fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
   fOutput->Add(fNtupleMC);
 
+  
   // Post output data.
   PostData(1,fOutput);
 
@@ -577,20 +591,17 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
   } 
   fHistNEvents->Fill(-1);
   
-  if(fLowEnergypp){ // remove events without SDD in pp 2.6 TeV
-    Bool_t hasSDD=kFALSE;
-    for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {  
-      track = (AliESDtrack*)fESD->GetTrack(iTrack);      
-      if (!track) continue;
-      UInt_t clumap = track->GetITSClusterMap();
-      if(clumap&4) hasSDD=kTRUE;
-      if(clumap&8) hasSDD=kTRUE;
-      if(hasSDD) break;
-    }
-    if(!hasSDD) return;
-  }
+  UInt_t maskPhysSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+  TString firedTriggerClasses=fESD->GetFiredTriggerClasses();
+  //  if(!firedTriggerClasses.Contains("CINT1B")) return;
+  if((maskPhysSel & AliVEvent::kMB)==0) return;
   fHistNEvents->Fill(0);
 
+  if(fLowEnergypp && !fMC){ // remove events without SDD in pp 2.76 TeV
+    if(!firedTriggerClasses.Contains("ALL")) return;
+  }
+  fHistNEvents->Fill(1);
+
 
   if(fMC){
     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
@@ -618,7 +629,7 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
   if(stack) nTrackMC = stack->GetNtrack();     
   const AliESDVertex *vtx =  fESD->GetPrimaryVertexSPD();
 
-///////////selection of the centrality or multiplicity bin
+  ///////////selection of the centrality or multiplicity bin
 
   //selection on the event centrality
   if(fHImode){
@@ -634,74 +645,39 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
   }
   
   //selection on the event multiplicity based on global tracks
-  Int_t multiplicity = fesdTrackCutsMult->CountAcceptedTracks(fESD);
-  if(!fSPD){
-    if(fLowMult>-1)
-      {
-       if(multiplicity<fLowMult)
-         return;
-      }
-    if(fUpMult>-1)
-      {
-       if(multiplicity>fUpMult)
-         return;
-      }
-    
-    Printf("Multiplicity of the event (global tracks) : %i",multiplicity);
-    Printf("Multiplicity cut (global tracks) : %i to %i",fLowMult,fUpMult);
-    fHistMult->Fill(multiplicity);
-  }
-  
-  //multipicity selection based on SPD
-  if(fSPD){
-    Float_t spdCorr=-1.0;
+  Int_t multiplicity = -1;  
+  if(fMultEstimator==0){
+    // tracks+tracklets
+    multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTrackletsITSTPC, 0.8); 
+  }else if(fMultEstimator==1){
+    // tracklets
+    multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTracklets, 0.8);    
+  }else if(fMultEstimator==2){
+    // clusters in SPD1
     const AliMultiplicity *mult = fESD->GetMultiplicity();
-    Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0};
-    for(Int_t ilay=0; ilay<6; ilay++)
-      {
-       nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
-      } 
-    spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vtx->GetZ());
-    {
-      if(fLowMult>-1)
-       {
-         if(((Int_t)spdCorr)<fLowMult)
-           {
-             return;
-           }   
-       }
-      if(fUpMult>-1)
-       {
-         if(((Int_t)spdCorr)>fUpMult)
-           {
-             return;
-           }           
-       }
-    }
-    
-    Printf("Multiplicity of the event (SPD) : %i",(Int_t)spdCorr);
-    Printf("Multiplicity cut (SPD) : %i to %i",fLowMult,fUpMult);
-    fHistMult->Fill(spdCorr);
+    Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1);
+    multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5);
   }
-  
-  
+  if(fLowMult>-1 && multiplicity<fLowMult) return;
+  if(fUpMult>-1 && multiplicity>fUpMult) return;
+  fHistMult->Fill(multiplicity);
+  fHistNEvents->Fill(2);
   
   //event selection
-  fHistNEvents->Fill(1);
   if(!vtx)evSel=0;
   else{
-    fHistNEvents->Fill(2);
+    fHistNEvents->Fill(3);
     if(vtx->GetNContributors()<0) evSel=0;
     else{
-      fHistNEvents->Fill(3);
+      fHistNEvents->Fill(4);
       if(TMath::Abs(vtx->GetZv())>10) evSel=0;
       else{
-       fHistNEvents->Fill(4);
+       fHistNEvents->Fill(5);
        if(vtx->GetZRes()>0.5) evSel=0;
        else{
-         fHistNEvents->Fill(5);
+         fHistNEvents->Fill(6);
          if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
-         else fHistNEvents->Fill(6);
+         else fHistNEvents->Fill(7);
        }
       }
     }
@@ -757,26 +733,26 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
     
     if(jpart>=0){
       if(stack->IsPhysicalPrimary(imc)){
-       if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC));
-       else  fHistPrimMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC));
+       if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(ptMC);
+       else  fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC);
        if(evSel==1){
-         if(signMC>0) fHistPrimMCpos[jpart]->Fill(TMath::Abs(ptMC));
-         else  fHistPrimMCneg[jpart]->Fill(TMath::Abs(ptMC));
+         if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC);
+         else  fHistPrimMCneg[jpart]->Fill(ptMC);
        }
       }else{
        if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay
-         if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC));
-         else  fHistSecStrMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC));            
+         if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC);
+         else  fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC);        
          if(evSel==1){
-           if(signMC>0) fHistSecStrMCpos[jpart]->Fill(TMath::Abs(ptMC));
-           else  fHistSecStrMCneg[jpart]->Fill(TMath::Abs(ptMC));          
+           if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC);
+           else  fHistSecStrMCneg[jpart]->Fill(ptMC);      
          }
        }else{
-         if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(TMath::Abs(ptMC));
-         else  fHistSecMatMCnegBefEvSel[jpart]->Fill(TMath::Abs(ptMC));            
+         if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC);
+         else  fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC);        
          if(evSel==1){
-           if(signMC>0) fHistSecMatMCpos[jpart]->Fill(TMath::Abs(ptMC));
-           else  fHistSecMatMCneg[jpart]->Fill(TMath::Abs(ptMC));          
+           if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC);
+           else  fHistSecMatMCneg[jpart]->Fill(ptMC);      
          }
        }
       }
@@ -1086,7 +1062,7 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
     }
     
     //DCA cut on xy and z
-    if(!DCAcut(impactXY,impactZ,pt,fMC)) continue;
+    if(!DCAcut(impactXY,impactZ,pt)) continue;
     
     label="DCA";
     fHistNTracks->Fill(countBinTrk);
@@ -1119,8 +1095,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
          else signMC=-1;
          ptMC=part->Pt();
          if(stack->IsPhysicalPrimary(track->GetLabel())){
-           if(signMC>0) fHistPrimMCposReco[jpart]->Fill(TMath::Abs(ptMC));
-           else  fHistPrimMCnegReco[jpart]->Fill(TMath::Abs(ptMC));
+           if(signMC>0) fHistPrimMCposReco[jpart]->Fill(ptMC);
+           else  fHistPrimMCnegReco[jpart]->Fill(ptMC);
          }else{ 
            Int_t indexMoth=part->GetFirstMother();
            if(indexMoth>=0){
@@ -1130,11 +1106,11 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
            }
            uniqueID = part->GetUniqueID();
            if(mfl==3 && uniqueID == kPDecay){ // strangeness
-             if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(TMath::Abs(ptMC));
-             else  fHistSecStrMCnegReco[jpart]->Fill(TMath::Abs(ptMC));            
+             if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC);
+             else  fHistSecStrMCnegReco[jpart]->Fill(ptMC);        
            }else{
-             if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(TMath::Abs(ptMC));
-             else  fHistSecMatMCnegReco[jpart]->Fill(TMath::Abs(ptMC));            
+             if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC);
+             else  fHistSecMatMCnegReco[jpart]->Fill(ptMC);        
            }
          }
        }
@@ -1265,7 +1241,8 @@ void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
                
   // Post output data.
   PostData(1,fOutput);
-  Printf("............. end of Exec");
+  PostData(2,fListCuts);
+  //  Printf("............. end of Exec");
 }      
 
 //________________________________________________________________________
@@ -1279,147 +1256,7 @@ void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
     return;
   } 
   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
-  fHistMult = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMult"));
-  fHistCen = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCen"));
-  fHistNTracks = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNTracks"));
-  fHistNTracksPos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNTracksPos"));
-  fHistNTracksNeg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNTracksNeg"));
-  fHistDEDX = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDX"));
-  fHistDEDXdouble = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDXdouble"));
-
-  
-  fHistBeforeEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBeforeEvSel"));
-  fHistAfterEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistAfterEvSel"));
-
-       
-  for(Int_t j=0;j<3;j++){
-    fHistPosNSigmaSep[j] = dynamic_cast<TH2F*>(fOutput->FindObject(Form("fHistPosNSigmaSep%d",j)));
-    fHistNegNSigmaSep[j] = dynamic_cast<TH2F*>(fOutput->FindObject(Form("fHistNegNSigmaSep%d",j)));
-    if(fMC){
-    fHistPrimMCpos[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCpos%d",j)));
-    fHistPrimMCneg[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCneg%d",j)));
-    fHistSecStrMCpos[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCpos%d",j)));
-    fHistSecStrMCneg[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCneg%d",j)));
-    fHistSecMatMCpos[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCpos%d",j)));
-    fHistSecMatMCneg[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCneg%d",j)));
-    //
-    fHistPrimMCposBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCposBefEvSel%d",j)));
-    fHistPrimMCnegBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCnegBefEvSel%d",j)));
-    fHistSecStrMCposBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCposBefEvSel%d",j)));
-    fHistSecStrMCnegBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCnegBefEvSel%d",j)));
-    fHistSecMatMCposBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCposBefEvSel%d",j)));
-    fHistSecMatMCnegBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCnegBefEvSel%d",j)));
-    //
-    fHistPrimMCposReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCposReco%d",j)));
-    fHistPrimMCnegReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPrimMCnegReco%d",j)));
-    fHistSecStrMCposReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCposReco%d",j)));
-    fHistSecStrMCnegReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecStrMCnegReco%d",j)));
-    fHistSecMatMCposReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCposReco%d",j)));
-    fHistSecMatMCnegReco[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistSecMatMCnegReco%d",j)));
-    }
-  }
-  
-  for(Int_t i=0; i<4; i++){
-    fHistCharge[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistChargeLay%d",i)));
-  }
-
-  for(Int_t i=0; i<kNbins; i++){
-    fHistPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosPi%d",i)));
-    fHistPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosK%d",i)));
-    fHistPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosP%d",i)));
-    fHistNegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegPi%d",i)));
-    fHistNegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegK%d",i)));
-    fHistNegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegP%d",i)));
-    
-    fHistDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosPi%d",i)));
-    fHistDCAPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosK%d",i)));
-    fHistDCAPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosP%d",i)));
-    fHistDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegPi%d",i)));
-    fHistDCANegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegK%d",i)));
-    fHistDCANegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegP%d",i)));
-    
-    
-    if(fMC){   
-      
-      fHistMCPrimDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCAPosPi%d",i)));
-      fHistMCPrimDCAPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCAPosK%d",i)));
-      fHistMCPrimDCAPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCAPosP%d",i)));
-      fHistMCPrimDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCANegPi%d",i)));
-      fHistMCPrimDCANegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCANegK%d",i)));
-      fHistMCPrimDCANegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPrimDCANegP%d",i)));  
-      
-      fHistMCSecStDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCAPosPi%d",i)));
-      fHistMCSecStDCAPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCAPosK%d",i)));
-      fHistMCSecStDCAPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCAPosP%d",i)));
-      fHistMCSecStDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCANegPi%d",i)));
-      fHistMCSecStDCANegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCANegK%d",i)));
-      fHistMCSecStDCANegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecStDCANegP%d",i)));  
-      
-      fHistMCSecMatDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCAPosPi%d",i)));
-      fHistMCSecMatDCAPosK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCAPosK%d",i)));
-      fHistMCSecMatDCAPosP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCAPosP%d",i)));
-      fHistMCSecMatDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCANegPi%d",i)));
-      fHistMCSecMatDCANegK[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCANegK%d",i)));
-      fHistMCSecMatDCANegP[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCSecMatDCANegP%d",i)));  
-      
-      fHistMCPosOtherHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosOtherHypPion%d",i)));
-      fHistMCPosOtherHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosOtherHypKaon%d",i)));
-      fHistMCPosOtherHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosOtherHypProton%d",i)));
-      fHistMCPosElHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosElHypPion%d",i)));
-      fHistMCPosElHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosElHypKaon%d",i)));
-      fHistMCPosElHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosElHypProton%d",i)));
-      fHistMCPosPiHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPiHypPion%d",i)));
-      fHistMCPosPiHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPiHypKaon%d",i)));
-      fHistMCPosPiHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPiHypProton%d",i)));
-      fHistMCPosKHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosKHypPion%d",i)));
-      fHistMCPosKHypKaon[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosKHypKaon%d",i)));
-      fHistMCPosKHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosKHypProton%d",i)));
-      fHistMCPosPHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPHypPion%d",i)));
-      fHistMCPosPHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPHypKaon%d",i)));
-      fHistMCPosPHypProton[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPHypProton%d",i)));
-      
-      fHistMCNegOtherHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegOtherHypPion%d",i)));
-      fHistMCNegOtherHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegOtherHypKaon%d",i)));
-      fHistMCNegOtherHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegOtherHypProton%d",i)));
-      fHistMCNegElHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegElHypPion%d",i)));
-      fHistMCNegElHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegElHypKaon%d",i)));
-      fHistMCNegElHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegElHypProton%d",i)));
-      fHistMCNegPiHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPiHypPion%d",i)));
-      fHistMCNegPiHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPiHypKaon%d",i)));
-      fHistMCNegPiHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPiHypProton%d",i)));
-      fHistMCNegKHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegKHypPion%d",i)));
-      fHistMCNegKHypKaon[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegKHypKaon%d",i)));
-      fHistMCNegKHypProton[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegKHypProton%d",i)));
-      fHistMCNegPHypPion[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPHypPion%d",i)));
-      fHistMCNegPHypKaon[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPHypKaon%d",i)));
-      fHistMCNegPHypProton[i]  = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPHypProton%d",i)));
-      
-    }
-  }
-  
-  for(Int_t j=0;j<3;j++){
-    fHistPosNSigmaMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaMean%d",j)));
-    fHistNegNSigmaMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaMean%d",j)));
-    fHistPosNSigmaMCMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaMCMean%d",j)));
-    fHistNegNSigmaMCMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaMCMean%d",j)));
-    fHistPosNSigmaPrimMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrimMean%d",j)));
-    fHistNegNSigmaPrimMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrimMean%d",j)));
-    fHistPosNSigmaPrimMCMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrimMCMean%d",j)));
-    fHistNegNSigmaPrimMCMean[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrimMCMean%d",j)));
-  
-    fHistPosNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigma%d",j)));
-    fHistNegNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigma%d",j)));
-    fHistPosNSigmaMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaMC%d",j)));
-    fHistNegNSigmaMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaMC%d",j)));
-    fHistPosNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrim%d",j)));
-    fHistNegNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrim%d",j)));
-    fHistPosNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrimMC%d",j)));
-    fHistNegNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrimMC%d",j)));
-  
-  }
-  
-  fNtupleNSigma = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleNSigma"));
-  fNtupleMC = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleMC"));
+  printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1));
   
   Printf("end of Terminate");
   return;
index fc31f0f..b535ff8 100644 (file)
 class TString;
 class TTree;
 class TH1F;
+class TF1;
 class TH2F;
 class TRandom3;
 class AliESDEvent;
 class TNtuple;
-class AliESDtrackCuts;
 class AliITSPIDResponse;
 
 #include "AliAnalysisTaskSE.h"
@@ -32,6 +32,8 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
   virtual ~AliAnalysisTaskSEITSsaSpectra();
 
   virtual void   UserCreateOutputObjects();
+  virtual void Init();
+  virtual void LocalInit() {Init();}
   virtual void   UserExec(Option_t *);
   virtual void   Terminate(Option_t *);
 
@@ -66,14 +68,16 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
       fLowCentrality=low; fUpCentrality=up;
     }
   }
- void SetSPDMethodCut(){fSPD=kTRUE;}
  void SetHImode(){fHImode=kTRUE;}
+ void SetUseTrackMultiplicityEstimator(){fMultEstimator=0;}
+ void SetUseTrackletsMultiplicityEstimator(){fMultEstimator=1;}
+ void SetUseClustersSPD1MultiplicityEstimator(){fMultEstimator=2;}
   
   void SetEtaMax(Double_t maxeta){
     fEtaRange=maxeta;
   }
   
-  void SetYear(Int_t year);
   void SetReadMC(Bool_t flag = kTRUE) {fMC = flag;}
   void SetFillNtuple(Bool_t fill=kTRUE) {fFillNtuple=fill;}
   void SetLowEnergypp(Bool_t opt=kTRUE) {fLowEnergypp=opt;}
@@ -85,7 +89,7 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
   }
   Double_t CookdEdx(Double_t *s) const; 
   Double_t Eta2y(Double_t pt, Double_t m, Double_t eta) const;
-  Bool_t DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC) const;
+  Bool_t DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const;
 
  private:
   AliAnalysisTaskSEITSsaSpectra(const AliAnalysisTaskSEITSsaSpectra &source); 
@@ -94,9 +98,9 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
   enum {kNbins=22};
   
   AliESDEvent *fESD; //ESD object
-  AliESDtrackCuts *fesdTrackCutsMult;//cuts for multiplicity 
   
   TList *fOutput; //! tlist with output
+  TList *fListCuts; // list of functions storing DCA cut 
   TH1F *fHistNEvents; //! histo with number of events
   TH1F *fHistMult; //! histo with multiplicity of the events
   TH1F *fHistCen; //! histo with multiplicity of the events
@@ -218,6 +222,9 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
   TH1F *fHistNegNSigmaPrim[3];   //! NSigma histos for 6 species
   TH1F *fHistNegNSigmaPrimMC[3]; //! NSigma histos for 6 species
 
+  TF1* fDCAxyCutFunc;  // function with DCAz cut vs. pt
+  TF1* fDCAzCutFunc;   // function with DCAxy cut vs. pt
+
   Double_t fPtBinLimits[kNbins+1]; // limits of Pt Bins
   AliITSPIDResponse* fITSPIDResponse; //! class with BB parameterizations
   Int_t    fMinSPDPts;       // minimum number of SPD Points
@@ -233,7 +240,7 @@ class AliAnalysisTaskSEITSsaSpectra : public AliAnalysisTaskSE {
   Int_t fUpMult;      // Multiplicity bin
   Float_t fLowCentrality;//low Centrality cut
   Float_t fUpCentrality;//up  Centrality cut
-  Bool_t fSPD;//use spd2 as mulestimator 
+  Int_t fMultEstimator; // multiplicty estimator 
   Bool_t fHImode;//use spd2 as mulestimator 
   Int_t fYear;        // Year (2009, 2010)
   Bool_t   fMC;        //flag to switch on the MC analysis for the efficiency estimation
diff --git a/PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C b/PWGLF/SPECTRA/PiKaPr/ITSsa/FitDCADistributions.C
new file mode 100644 (file)
index 0000000..3a4e128
--- /dev/null
@@ -0,0 +1,303 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TFile.h>
+#include <TF1.h>
+#include <TH1F.h>
+#include <TCanvas.h>
+#include <TLatex.h>
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TFractionFitter.h>
+#endif
+
+enum species {kPion, kKaon, kProton};
+enum chargesign {kPositive,kNegative};
+
+Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax);
+
+void FitDCADistributions(TString period="LHC10d",
+                        TString MCname="LHC10f6aold",  
+                        Int_t iSpecies=2, 
+                        Int_t cSign=1
+                        )
+{
+  gROOT->SetStyle("Plain");
+  gStyle->SetFillColor(0);
+  gStyle->SetOptStat(0000);
+  //binning
+  const Int_t nptbins = 22;
+  Double_t ptbins[nptbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
+  //Final correction
+  TH1F *hPrimAllDATA=new TH1F("hPrimAllDATA","hPrimAllDATA",nptbins,ptbins);
+  TH1F *hPrimAllMC=new TH1F("hPrimAllMC","hPrimAllMC",nptbins,ptbins);
+  TH1F *hPrimAllDATAMC=new TH1F("hPrimAllDATAMC","hPrimAllDATAMC",nptbins,ptbins);
+  
+  
+  TString fnameMC=Form("/home/prino/alice/PID/pp7/%s/AnalysisResults.root",MCname.Data());
+  TString lnameMC="clistITSsaMult-1to-1";
+  TString lnameCutMC="DCACutMult-1to-1";
+  TString fnameDATA=Form("/home/prino/alice/PID/pp7/%s/AnalysisResults.root",period.Data());
+  TString lnameDATA="clistITSsaMult-1to-1";
+  TString lnameCutDATA="DCACutMult-1to-1";
+
+  const Int_t firstbin=1;
+  Int_t rebin=10;
+  Bool_t optSmooth=kTRUE;
+
+
+  TString pCharge="Pos";
+  TString pcharge="pos";
+  if(cSign==kNegative){
+    pCharge="Neg";
+    pcharge="neg";
+  }
+  TString species[3]={"Pi","K","P"};
+  TString speciesRoot[6]={"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"};
+  Int_t idPart=iSpecies*2+cSign;
+  TString partname=Form("%s%s",pCharge.Data(),species[iSpecies].Data());
+  printf("%s\n",partname.Data());
+
+  // MC histograms 
+  TFile *fMC=new TFile(fnameMC,"READ");  
+  TH1F *fHistDCAMC[nptbins]; 
+  TH1F *fHistDCAPrim[nptbins]; 
+  TH1F *fHistDCAsec[nptbins]; 
+  TH1F *fHistDCAsecSt[nptbins]; 
+  TDirectory *dirFileMC=(TDirectory*)fMC->Get("PWG2SpectraITSsa");
+  TList *liMC = (TList*)dirFileMC->Get(lnameMC.Data());
+  for(Int_t m=0;m<nptbins;m++){
+    fHistDCAMC[m] = (TH1F*)liMC->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
+    fHistDCAPrim[m] = (TH1F*)liMC->FindObject(Form("fHistMCPrimDCA%s%i",partname.Data(),m));
+    fHistDCAsec[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecMatDCA%s%i",partname.Data(),m));
+    fHistDCAsecSt[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecStDCA%s%i",partname.Data(),m));
+    fHistDCAMC[m]->Rebin(rebin);
+    fHistDCAPrim[m]->Rebin(rebin);
+    fHistDCAsec[m]->Rebin(rebin);
+    fHistDCAsecSt[m]->Rebin(rebin);
+  }
+  TList *liCutMC = (TList*)dirFileMC->Get(lnameCutMC.Data());
+  Bool_t setDCA=kFALSE;
+
+  Float_t DCAcut=7.;
+  Double_t xyPMC[3];   //parameters of DCA cut
+  xyPMC[0]=36.; //MC LHC10d1
+  xyPMC[1]=43.9;
+  xyPMC[2]=1.3;
+  TF1* fDCAxyCutMC=0x0;
+  TF1* fDCAzCutMC=0x0;
+  Float_t nSigmaDCAxyMC=0.;
+  Float_t nSigmaDCAzMC=0.;
+  if(liCutMC){
+    fDCAxyCutMC=(TF1*)liCutMC->FindObject("fDCAxyCutFunc");
+    fDCAzCutMC=(TF1*)liCutMC->FindObject("fDCAzCutFunc");
+    nSigmaDCAxyMC=fDCAxyCutMC->GetParameter(3);
+    nSigmaDCAzMC=fDCAzCutMC->GetParameter(3);
+    printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutMC->GetParameter(0),fDCAxyCutMC->GetParameter(1),fDCAxyCutMC->GetParameter(2),nSigmaDCAxyMC);
+    setDCA=kTRUE;
+  }
+  if(!setDCA){
+    printf("DCA cut values for MC not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
+    fDCAxyCutMC = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
+    for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutMC->SetParameter(ipar,xyPMC[ipar]);
+    fDCAxyCutMC->SetParameter(3,DCAcut);
+    nSigmaDCAxyMC=DCAcut;
+ }
+
+  //DATA histograms
+  TFile *fDATA=new TFile(fnameDATA,"READ");
+  TDirectory *dirFileDATA=(TDirectory*)fDATA->Get("PWG2SpectraITSsa");
+  TList *liDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
+  TH1F *fHistDCADATA[nptbins]; 
+  for(Int_t m=0;m<nptbins;m++){
+    fHistDCADATA[m] = (TH1F*)liDATA->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
+    fHistDCADATA[m]->Rebin(rebin);
+  }
+  TList *liCutDATA = (TList*)dirFileDATA->Get(lnameCutDATA.Data());
+  setDCA=kFALSE;
+  Float_t nSigmaDCAxyDATA=0.;
+  Float_t nSigmaDCAzDATA=0.;
+  Double_t xyP[3];   // parameters of DCA cut
+  xyP[0]=32.7;  //DATA 7 TeV pass2
+  xyP[1]=44.8;
+  xyP[2]=1.3;
+  TF1* fDCAxyCutDATA=0x0;
+  TF1* fDCAzCutDATA=0x0;
+  if(liCutDATA){
+    fDCAxyCutDATA=(TF1*)liCutDATA->FindObject("fDCAxyCutFunc");
+    fDCAzCutDATA=(TF1*)liCutDATA->FindObject("fDCAzCutFunc");
+    nSigmaDCAxyDATA=fDCAxyCutDATA->GetParameter(3);
+    nSigmaDCAzDATA=fDCAzCutDATA->GetParameter(3);
+    printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutDATA->GetParameter(0),fDCAxyCutDATA->GetParameter(1),fDCAxyCutDATA->GetParameter(2),nSigmaDCAxyDATA);
+    setDCA=kTRUE;
+  }
+  if(!setDCA){
+    printf("DCA cut values for DATA not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
+    fDCAxyCutDATA = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
+    for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutDATA->SetParameter(ipar,xyP[ipar]);
+    fDCAxyCutDATA->SetParameter(3,DCAcut);
+    nSigmaDCAxyDATA=DCAcut;
+  }
+
+  if(TMath::Abs(nSigmaDCAxyDATA-nSigmaDCAxyMC)<0.01) DCAcut=nSigmaDCAxyDATA;
+  else{
+    printf("ERROR: DCAxy cuts do not match between data and MC\n");
+    return;
+  }
+  if(TMath::Abs(nSigmaDCAzDATA-nSigmaDCAzMC)<0.01){
+    printf("ERROR: DCAz cuts do not match between data and MC\n");
+    return;
+  }
+  TCanvas *cfitDATA=new TCanvas("cfitDATA","cfitDATA");
+  cfitDATA->Divide(6,4);
+  TCanvas *cfitMC=new TCanvas("cfitMC","cfitMC");
+  cfitMC->Divide(6,4);
+
+  for(Int_t ibin=firstbin;ibin<nptbins;ibin++){    
+    
+    Double_t xyMax =fDCAxyCutDATA->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;
+    Double_t xyMaxMC =fDCAxyCutMC->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;
+
+    TObjArray *mcTemplates = 0x0;
+    if(partname=="PosP"){
+      mcTemplates = new TObjArray(3);        // MC histograms are put in this array
+      mcTemplates->Add(fHistDCAPrim[ibin]);
+      mcTemplates->Add(fHistDCAsecSt[ibin]);
+      mcTemplates->Add(fHistDCAsec[ibin]);
+    }else{
+      mcTemplates = new TObjArray(2);        // MC histograms are put in this array
+      mcTemplates->Add(fHistDCAPrim[ibin]);
+      mcTemplates->Add(fHistDCAsecSt[ibin]);
+    }
+    if(fHistDCADATA[ibin]->Integral()==0) continue;
+
+    ////////////////// Fit on DATA
+    cfitDATA->cd(ibin);
+    gPad->SetLogy();
+    Double_t primAllDATA=PerformFit(fHistDCADATA[ibin],mcTemplates,partname,xyMax);
+    
+    ////////////////// Fit on MC
+    cfitMC->cd(ibin);
+    gPad->SetLogy();
+    Double_t primAllMC=PerformFit(fHistDCAMC[ibin],mcTemplates,partname,xyMaxMC);
+    
+    if(hPrimAllMC!=0 && hPrimAllDATA!=0){
+      hPrimAllDATA->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA);
+      hPrimAllMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllMC);
+      hPrimAllDATAMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA/primAllMC);
+    }
+  }
+  
+  //Adding MC truth
+  TH1F* hAll=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaMean%d",pCharge.Data(),iSpecies));
+  TH1F* hPrim=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaPrimMean%d",pCharge.Data(),iSpecies));
+  TH1F* hPrimRecMC=(TH1F*)liMC->FindObject(Form("fHistPrimMC%sReco%d",pcharge.Data(),iSpecies));
+  TH1F* hSecMatRecMC=(TH1F*)liMC->FindObject(Form("fHistSecMatMC%sReco%d",pcharge.Data(),iSpecies));
+  TH1F* hSecStrRecMC=(TH1F*)liMC->FindObject(Form("fHistSecStrMC%sReco%d",pcharge.Data(),iSpecies));
+  hSecStrRecMC->Add(hSecMatRecMC);
+  hSecStrRecMC->Add(hPrimRecMC);
+  hPrim->SetLineColor(8);
+  hPrim->SetLineWidth(2);
+  hPrim->SetTitle("Prim/All from MC Truth");
+  hPrim->Divide(hAll);
+  hPrimRecMC->SetLineColor(11);
+  hPrimRecMC->SetLineWidth(2);
+  hPrimRecMC->SetTitle("Prim/All from MC Truth Reco");
+  hPrimRecMC->Divide(hSecStrRecMC);
+
+  hPrimAllDATAMC->GetXaxis()->SetRangeUser(0.09,0.79);
+  hPrimAllMC->GetXaxis()->SetRangeUser(0.09,0.79);
+  hPrimAllDATA->GetXaxis()->SetRangeUser(0.09,0.79);
+  if(optSmooth) hPrimAllDATAMC->Smooth(1,"R");
+  
+  hPrimAllDATA->SetTitle("Prim/all data");
+  hPrimAllDATA->SetLineColor(2);
+  hPrimAllDATA->SetLineWidth(2);
+  hPrimAllMC->SetTitle("Prim/all mc");
+  hPrimAllMC->SetLineColor(4);
+  hPrimAllMC->SetLineWidth(2);
+  hPrimAllDATAMC->SetTitle("Prim/all data/mc");
+  hPrimAllDATAMC->SetLineColor(1);
+  hPrimAllDATAMC->SetLineWidth(2);
+
+  TCanvas *cPrimAll=new TCanvas("cPrimAll","cPrimAll");
+  hPrimAllDATA->SetMinimum(0.7);
+  hPrimAllDATA->SetMaximum(1.1);
+  hPrimAllDATA->Draw("l");
+  hPrimAllMC->Draw("lsame");
+  hPrimAllDATAMC->Draw("lsame");
+  hPrim->Draw("lsame");
+  hPrimRecMC->Draw("lsame");
+  gPad->BuildLegend();
+  TLatex* tsp=new TLatex(0.4,0.75,speciesRoot[idPart].Data());
+  tsp->SetTextFont(43);
+  tsp->SetTextSize(28);
+  tsp->SetNDC();
+  tsp->Draw();
+  cPrimAll->Update();
+
+  TString fout;
+  fout=Form("DCACorr%s_%s_%.0fDCA_%s_TFraction.root",period.Data(),MCname.Data(),DCAcut,partname.Data());
+  TFile *out=new TFile(fout.Data(),"recreate");
+  hPrimAllDATA->Write();
+  hPrimAllMC->Write();
+  hPrimAllDATAMC->Write();
+  out->Close();
+  delete out;
+  
+}
+
+Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax){
+  
+  Double_t prim=0,secSt=0,sec=0;
+  TFractionFitter* fit = new TFractionFitter(fHistDCA, mc); // initialise
+  fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
+  fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
+  if(partname=="PosP")fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
+  
+  //fit->SetRangeX(200,1000);                    // use only the first 15 bins in the fit
+  Int_t status = fit->Fit();               // perform the fit
+  cout << "fit status: " << status << endl;
+  
+  if (status == 0) {                       // check on fit status
+    TH1F* result = (TH1F*) fit->GetPlot();
+    TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
+    TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
+    TH1F* secMCPred=0x0;
+    if(partname=="PosP"){
+      secMCPred=(TH1F*)fit->GetMCPrediction(2);
+      secMCPred->SetLineColor(4);
+    }
+    //Drawing section
+    PrimMCPred->SetLineColor(2);
+    secStMCPred->SetLineColor(6);
+    fHistDCA->SetTitle("DCA distr - data");
+    fHistDCA->DrawCopy("Ep");
+    result->SetTitle("Fit result");
+    result->DrawCopy("same");
+    Double_t value,error;
+    
+    fit->GetResult(0,value,error);
+    PrimMCPred->Scale(fHistDCA->GetSumOfWeights()*value/PrimMCPred->GetSumOfWeights());
+    PrimMCPred->SetTitle("Primaries");
+    PrimMCPred->DrawCopy("same");
+    fit->GetResult(1,value,error);
+    secStMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secStMCPred->GetSumOfWeights());
+    secStMCPred->SetTitle("Sec from strangeness");
+    secStMCPred->DrawCopy("same");
+    if(partname=="PosP" && secMCPred){
+      fit->GetResult(2,value,error);
+      secMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secMCPred->GetSumOfWeights());
+      secMCPred->SetTitle("Sec from material");
+      secMCPred->DrawCopy("same");
+    }
+    prim=PrimMCPred->Integral(PrimMCPred->FindBin(-xyMax),PrimMCPred->FindBin(xyMax));
+    secSt=secStMCPred->Integral(secStMCPred->FindBin(-xyMax),secStMCPred->FindBin(xyMax));
+    if(partname=="PosP")sec=secMCPred->Integral(secMCPred->FindBin(-xyMax),secMCPred->FindBin(xyMax));
+  }
+  else{
+    prim=1;
+    secSt=1;
+    if(partname=="PosP")sec=1;
+  }
+  Printf("Yields:  primary=%f material=%f strange=%f\n",prim,sec,secSt);
+  return prim/(prim+secSt+sec);
+}