Update of cut variation function

diff --git a/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx b/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx
index 2078cad739fd7b70949887224492ee3c82e85ddf..0279b92c6923808821a130103a68a0c3e5cbbd4c 100644 (file)
--- a/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx
+++ b/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx
@@ -77,7 +77,7 @@ AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
-    fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0)
+    fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
 
 {
 
@@ -100,7 +100,7 @@ AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
-    fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0)
+    fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
 
 {
 
@@ -342,7 +342,7 @@ void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
   fListHist->Add(fHistPsi2sMassVsPt);
   
-  fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
+  fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",50,2.5,5.5);
   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
   fListHist->Add(fHistPsi2sMassCoherent);
   
@@ -381,6 +381,10 @@ fListJPsiLoose->Add(fHistJPsiDCAzLoose);
 TH1D *fHistJPsiDCAxyLoose = new TH1D("JPsiDCAxyLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 fListJPsiLoose->Add(fHistJPsiDCAxyLoose);
 
+TH1D *fHistJPsiITShitsLoose = new TH1D("JPsiITShitsLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
+fListJPsiLoose->Add(fHistJPsiITShitsLoose);
+
+
 fListJPsiTight = new TList();
 fListJPsiTight->SetOwner();
 fListJPsiTight->SetName("JPsiTight");
@@ -399,6 +403,45 @@ TH1D *fHistJPsiDCAxyTight = new TH1D("JPsiDCAxyTight","Invariant mass of J/#psi
 fListJPsiTight->Add(fHistJPsiDCAxyTight);
 
 
+fListPsi2sLoose = new TList();
+fListPsi2sLoose->SetOwner();
+fListPsi2sLoose->SetName("Psi2sLoose");
+fListSystematics->Add(fListPsi2sLoose);
+
+TH1D *fHistPsi2sNClusLoose = new TH1D("Psi2sNClusLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sLoose->Add(fHistPsi2sNClusLoose);
+
+TH1D *fHistPsi2sChi2Loose = new TH1D("Psi2sChi2Loose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sLoose->Add(fHistPsi2sChi2Loose);
+
+TH1D *fHistPsi2sDCAzLoose = new TH1D("Psi2sDCAzLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sLoose->Add(fHistPsi2sDCAzLoose);
+
+TH1D *fHistPsi2sDCAxyLoose = new TH1D("Psi2sDCAxyLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sLoose->Add(fHistPsi2sDCAxyLoose);
+
+TH1D *fHistPsi2sITShitsLoose = new TH1D("Psi2sITShitsLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sLoose->Add(fHistPsi2sITShitsLoose);
+
+
+fListPsi2sTight = new TList();
+fListPsi2sTight->SetOwner();
+fListPsi2sTight->SetName("Psi2sTight");
+fListSystematics->Add(fListPsi2sTight);
+
+TH1D *fHistPsi2sNClusTight = new TH1D("Psi2sNClusTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sTight->Add(fHistPsi2sNClusTight);
+
+TH1D *fHistPsi2sChi2Tight = new TH1D("Psi2sChi2Tight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sTight->Add(fHistPsi2sChi2Tight);
+
+TH1D *fHistPsi2sDCAzTight = new TH1D("Psi2sDCAzTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sTight->Add(fHistPsi2sDCAzTight);
+
+TH1D *fHistPsi2sDCAxyTight = new TH1D("Psi2sDCAxyTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
+fListPsi2sTight->Add(fHistPsi2sDCAxyTight);
+
+
 }
 
 //_____________________________________________________________________________
@@ -477,7 +520,7 @@ void AliAnalysisTaskUpcPsi2s::RunAODhist()
   AliAODEvent *aod = (AliAODEvent*) InputEvent();
   if(!aod) return;
   
-  cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
+  //cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
 
   fHistNeventsJPsi->Fill(1);
   fHistNeventsPsi2s->Fill(1);
@@ -526,10 +569,12 @@ void AliAnalysisTaskUpcPsi2s::RunAODhist()
   
   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
   Short_t qLepton[4], qPion[4];
-  UInt_t nLepton=0, nPion=0, nHighPt=0;
+  UInt_t nLepton=0, nPion=0, nHighPt=0, nSpdHits=0;
   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
   Int_t fChannel = 0;
   Int_t mass[3]={-1,-1,-1};
+  Double_t TrackPt[5]={0,0,0,0,0};
+  Double_t MeanPt = -1;
   
    
   //Four track loop
@@ -546,10 +591,13 @@ void AliAnalysisTaskUpcPsi2s::RunAODhist()
       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
       delete trk_clone;
-
       if(TMath::Abs(dca[1]) > 2) continue;
+      Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
+      if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
+      if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
      
       TrackIndex[nGoodTracks] = itr;
+      TrackPt[nGoodTracks] = trk->Pt();
       nGoodTracks++;
                                  
       if(nGoodTracks > 4) break;  
@@ -558,12 +606,13 @@ void AliAnalysisTaskUpcPsi2s::RunAODhist()
   nLepton=0; nPion=0; nHighPt=0;
   mass[0]= -1; mass[1]= -1, mass[2]= -1;
   
-  if(nGoodTracks == 4){
+  if(nGoodTracks == 4 && nSpdHits>1){
+         MeanPt = GetMedian(TrackPt);
          fHistNeventsPsi2s->Fill(6);
          for(Int_t i=0; i<4; i++){
                AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
                
-               if(trk->Pt() > 1){   
+               if(trk->Pt() > MeanPt){   
                        fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
                        qLepton[nLepton] = trk->Charge();
                        if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
@@ -591,14 +640,26 @@ void AliAnalysisTaskUpcPsi2s::RunAODhist()
                if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
                if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
                        fHistNeventsPsi2s->Fill(11);
-                       if(mass[0] == mass[1]) {
+                       if(mass[0] != -1 && mass[1] != -1) {
                                fHistNeventsPsi2s->Fill(12); 
                                vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
                                vDilepton = vLepton[0]+vLepton[1];
                                fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
-                               if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
-                               if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);   
-                               if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
+                               fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
+                               if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
+                               else { 
+                                       fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
+                                       if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1; 
+                                       }
+                               
+                               if(fChannel == -1) {
+                                       fHistNeventsPsi2s->Fill(13);
+                                       if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
+                                       }       
+                               if(fChannel == 1){ 
+                                       fHistNeventsPsi2s->Fill(14);
+                                       if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) fHistPsi2sMassCoherent->Fill(vCandidate.M());
+                                       }
                                }
                        }
                }
@@ -724,9 +785,9 @@ void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
   Int_t nGoodTracks = 0;
   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
   
-  TLorentzVector vLepton[4], vCandidate, vDilepton;
-  Short_t qLepton[4];
-  UInt_t nLepton=0, nHighPt=0;
+  TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
+  Short_t qLepton[4],qPion[4];
+  UInt_t nLepton=0, nPion=0, nHighPt=0;
   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
   Int_t fChannel = 0;
 
@@ -740,12 +801,12 @@ void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
   Double_t electronMass = partElectron->Mass();
   
- // TParticlePDG *partPion = pdgdat->GetParticle( 211 );
- // Double_t pionMass = partPion->Mass();
+  TParticlePDG *partPion = pdgdat->GetParticle( 211 );
+  Double_t pionMass = partPion->Mass();
 
   
-for(Int_t i=0; i<4; i++){
-         cout<<"Loose sytematics, cut"<<i<<endl;
+for(Int_t i=0; i<5; i++){
+         //cout<<"Loose sytematics, cut"<<i<<endl;
          for(Int_t j=0; j<4; j++){
                  if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
                  else fJPsiSels[j] = fJPsiSelsMid[j];
@@ -759,7 +820,7 @@ for(Int_t i=0; i<4; i++){
 
       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
-      if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
+      if(i!=4){ if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;}
       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
@@ -782,7 +843,6 @@ for(Int_t i=0; i<4; i++){
   nLepton=0; nHighPt=0;
   
   if(nGoodTracks == 2){
-         fHistNeventsJPsi->Fill(6);
          for(Int_t k=0; k<2; k++){
                AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);                
                if(trk->Pt() > 1) nHighPt++;     
@@ -815,7 +875,7 @@ for(Int_t i=0; i<4; i++){
 }//loose cuts
 
 for(Int_t i=0; i<4; i++){
-         cout<<"Tight sytematics, cut"<<i<<endl;
+         //cout<<"Tight sytematics, cut"<<i<<endl;
          for(Int_t j=0; j<4; j++){
                  if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
                  else fJPsiSels[j] = fJPsiSelsMid[j];
@@ -853,7 +913,6 @@ for(Int_t i=0; i<4; i++){
   nLepton=0; nHighPt=0;
   
   if(nGoodTracks == 2){
-         fHistNeventsJPsi->Fill(6);
          for(Int_t k=0; k<2; k++){
                AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);                
                if(trk->Pt() > 1) nHighPt++;     
@@ -885,6 +944,197 @@ for(Int_t i=0; i<4; i++){
   }
 }//tight cuts
 
+//---------------------------------------------Psi2s------------------------------------------------------------------------
+
+  Double_t fPsi2sSels[4];
+
+  fPsi2sSels[0] =   50; //min number of TPC clusters
+  fPsi2sSels[1] =   4; //chi2
+  fPsi2sSels[2] =   2; //DCAz
+  fPsi2sSels[3] =   4; // DCAxy 1x 
+
+  Double_t fPsi2sSelsMid[4];
+
+  fPsi2sSelsMid[0] =   50; //min number of TPC clusters
+  fPsi2sSelsMid[1] =   4; //chi2
+  fPsi2sSelsMid[2] =   2; //DCAz
+  fPsi2sSelsMid[3] =   4; // DCAxy 1x 
+  
+  Double_t fPsi2sSelsLoose[4];
+
+  fPsi2sSelsLoose[0] =   60; //min number of TPC clusters
+  fPsi2sSelsLoose[1] =   5; //chi2
+  fPsi2sSelsLoose[2] =   3; //DCAz
+  fPsi2sSelsLoose[3] =   6; // DCAxy 2x 
+
+  Double_t fPsi2sSelsTight[4];
+
+  fPsi2sSelsTight[0] =   70; //min number of TPC clusters
+  fPsi2sSelsTight[1] =   3.5; //chi2
+  fPsi2sSelsTight[2] =   1; //DCAz
+  fPsi2sSelsTight[3] =   2; // DCAxy 0.5x 
+
+  nGoodTracks = 0; nLepton=0; nHighPt=0; fChannel = 0;
+  Int_t nSpdHits = 0;
+  Double_t TrackPt[5]={0,0,0,0,0};
+  Double_t MeanPt = -1;
+
+for(Int_t i=0; i<5; i++){
+         //cout<<"Loose systematics psi2s, cut"<<i<<endl;
+         for(Int_t j=0; j<4; j++){
+                 if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
+                 else fJPsiSels[j] = fJPsiSelsMid[j];
+         }
+ 
+  //Two track loop
+  for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
+    AliAODTrack *trk = aod->GetTrack(itr);
+    if( !trk ) continue;
+    if(!(trk->TestFilterBit(1<<0))) continue;
+
+      if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
+      if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
+      if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
+      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
+      AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
+      if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+      delete trk_clone;
+      Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
+      
+      if(trk->GetTPCNcls() < fJPsiSels[0])continue;
+      if(trk->Chi2perNDF() > fJPsiSels[1])continue;
+      if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;      
+      if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
+      if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
+     
+      TrackIndex[nGoodTracks] = itr;
+      TrackPt[nGoodTracks] = trk->Pt();
+      nGoodTracks++;
+                                 
+      if(nGoodTracks > 4) break;  
+  }//Track loop
+    
+  Int_t mass[3]={-1,-1,-1};
+  fChannel = 0;
+  nLepton=0; nHighPt=0;
+  
+  if(nGoodTracks == 4){
+         if(i!=4){ if(nSpdHits<2) continue;} 
+         MeanPt = GetMedian(TrackPt);
+         for(Int_t k=0; k<4; k++){
+               AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
+               if(trk->Pt() > MeanPt){   
+                       fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
+                       qLepton[nLepton] = trk->Charge();
+                       if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
+                                       vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
+                                       mass[nLepton] = 0;
+                                       }
+                       if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
+                                       vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
+                                       mass[nLepton] = 1;
+                                       }
+                       nLepton++;
+                       }
+               else{
+                       qPion[nPion] = trk->Charge();
+                       vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
+                       nPion++;
+                       }             
+               }
+       if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
+               vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
+               vDilepton = vLepton[0]+vLepton[1];
+               fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
+               if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
+               else { 
+                       fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
+                       if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1; 
+                       }                       
+               if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());              
+               if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
+       }
+  }   
+}//loose cuts
+
+for(Int_t i=0; i<4; i++){
+         //cout<<"Tight systematics psi2s, cut"<<i<<endl;
+         for(Int_t j=0; j<4; j++){
+                 if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
+                 else fJPsiSels[j] = fJPsiSelsMid[j];
+         }
+ 
+  //Two track loop
+  for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
+    AliAODTrack *trk = aod->GetTrack(itr);
+    if( !trk ) continue;
+    if(!(trk->TestFilterBit(1<<0))) continue;
+
+      if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
+      if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
+      if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
+      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
+      AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
+      if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+      delete trk_clone;
+      Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
+      
+      if(trk->GetTPCNcls() < fJPsiSels[0])continue;
+      if(trk->Chi2perNDF() > fJPsiSels[1])continue;
+      if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;      
+      if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
+      if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
+     
+      TrackIndex[nGoodTracks] = itr;
+      TrackPt[nGoodTracks] = trk->Pt();
+      nGoodTracks++;
+                                 
+      if(nGoodTracks > 4) break;  
+  }//Track loop
+    
+  Int_t mass[3]={-1,-1,-1};
+  fChannel = 0;
+  nLepton=0; nHighPt=0;
+  
+  if(nGoodTracks == 4){
+         if(nSpdHits<2) continue; 
+         MeanPt = GetMedian(TrackPt);
+         for(Int_t k=0; k<4; k++){
+               AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
+               if(trk->Pt() > MeanPt){   
+                       fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
+                       qLepton[nLepton] = trk->Charge();
+                       if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
+                                       vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
+                                       mass[nLepton] = 0;
+                                       }
+                       if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
+                                       vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
+                                       mass[nLepton] = 1;
+                                       }
+                       nLepton++;
+                       }
+               else{
+                       qPion[nPion] = trk->Charge();
+                       vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
+                       nPion++;
+                       }             
+               }
+       if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
+               vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
+               vDilepton = vLepton[0]+vLepton[1];
+               fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
+               if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
+               else { 
+                       fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
+                       if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1; 
+                       }                       
+               if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());              
+               if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
+       }
+  }   
+}//Tight cuts
+
 
 }
 //_____________________________________________________________________________
@@ -1079,7 +1329,6 @@ void AliAnalysisTaskUpcPsi2s::RunAODtree()
       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
       delete trk_clone;
-      if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
       if(TMath::Abs(dca[1]) > 2) continue;
       Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
@@ -1532,7 +1781,7 @@ void AliAnalysisTaskUpcPsi2s::RunESDtree()
       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
       trk->GetImpactParameters(dca[0],dca[1]);
       if(TMath::Abs(dca[1]) > 2) continue;
-      if(TMath::Abs(dca[1]) > 0.2) continue;
+      if(TMath::Abs(dca[0]) > 0.2) continue;
       
       TrackIndex[nGoodTracks] = itr;
       nGoodTracks++;
@@ -1662,3 +1911,26 @@ void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
   cout<<"Analysis complete."<<endl;
 }//Terminate
 
+//_____________________________________________________________________________
+Double_t AliAnalysisTaskUpcPsi2s::GetMedian(Double_t *daArray) {
+    // Allocate an array of the same size and sort it.
+    Double_t dpSorted[4];
+    for (Int_t i = 0; i < 4; ++i) {
+        dpSorted[i] = daArray[i];
+    }
+    for (Int_t i = 3; i > 0; --i) {
+        for (Int_t j = 0; j < i; ++j) {
+            if (dpSorted[j] > dpSorted[j+1]) {
+                Double_t dTemp = dpSorted[j];
+                dpSorted[j] = dpSorted[j+1];
+                dpSorted[j+1] = dTemp;
+            }
+        }
+    }
+
+    // Middle or average of middle values in the sorted array.
+    Double_t dMedian = 0.0;
+    dMedian = (dpSorted[2] + dpSorted[1])/2.0;
+    
+    return dMedian;
+}

diff --git a/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.h b/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.h
index fa85e3c0b7ae7f2e2645b0c5f270ba1ce78c1b54..f0b4737bb16be69b1e306fc026e4346a0d934f69 100644 (file)
--- a/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.h
+++ b/PWGUD/UPC/AliAnalysisTaskUpcPsi2s.h
@@ -41,6 +41,7 @@ class AliAnalysisTaskUpcPsi2s : public AliAnalysisTaskSE {
   void SetRunSyst(Bool_t runSyst){fRunSystematics = runSyst;}
   void SetIsMC(Bool_t MC){isMC = MC;}
   void InitSystematics();
+  Double_t GetMedian(Double_t *daArray);
 
  private:
   Int_t fType; // 0 - ESD, 1 - AOD
@@ -112,6 +113,8 @@ class AliAnalysisTaskUpcPsi2s : public AliAnalysisTaskSE {
   TList *fListSystematics;
   TList *fListJPsiLoose;
   TList *fListJPsiTight;
+  TList *fListPsi2sLoose;
+  TList *fListPsi2sTight;
   
   AliAnalysisTaskUpcPsi2s(const AliAnalysisTaskUpcPsi2s&); //not implemented
   AliAnalysisTaskUpcPsi2s& operator =(const AliAnalysisTaskUpcPsi2s&); //not implemented