classification of lost clusters detailed

diff --git a/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.cxx b/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.cxx
index a2e2cfc9a412f5b9ee29185a78157a4846db26b5..55a71e415b1067e4ef2bfd96945fe1cbe1aef1b4 100644 (file)
--- a/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.cxx
+++ b/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.cxx
@@ -375,7 +375,8 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
     
          //MC tagging: reasons of partner loss etc.
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_%s_cent%d",cPID[iPID],cen),"Decay photons with partner not in Stack", nPt,0.,ptMax )) ;
-         fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartn_%s_cent%d",cPID[iPID],cen),"Decay photon with found partner", nPt,0.,ptMax )) ;
+         for(Int_t iType=0; iType<9; iType++)
+          fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartnType%d_%s_cent%d",iType,cPID[iPID],cen),"Decay photon with found partner", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_%s_cent%d",cPID[iPID],cen),"Decay photon with wrong mass", nPt,0.,ptMax )) ;       
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_%s_cent%d",cPID[iPID],cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax )) ;
@@ -387,13 +388,13 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_%s_cent%d",cPID[iPID],cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax )) ;
 
-       
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax )) ;
          fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
+         fOutputContainer->Add(new TH1F(Form("hMCDecWRecUniqPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
 
         fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
          fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
@@ -717,7 +718,7 @@ void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
     p->SetIsolationTag(isolation) ;
     
     p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
-    
+    p->SetBC(i) ; //reference to CaloCluster
     p->SetTagInfo(0); //Strict PID pi0 partner not found
     p->SetTagged(kFALSE);   //Reconstructed pairs found
     
@@ -785,6 +786,8 @@ void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
   const Double_t rcut=1. ; //cut on vertex to consider particle as "primary" 
   const Double_t phiMin=260.*TMath::Pi()/180. ;
   const Double_t phiMax=320.*TMath::Pi()/180. ;
+
+  AliVEvent* event = (AliVEvent*)InputEvent();
   
   Int_t nPrim = fStack->GetEntriesFast() ;
   //Fill Primary particl yields
@@ -983,16 +986,51 @@ void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
                break ;
              }
            }
+           if(pp){
+             //Partner reconstructed, but did not pass cuts
+                FillPIDHistograms("hMCDecWRecUniqPartn",p) ;   
+           }
+           //Partner not found. Check if it is not dominant contributor?
+           if(!pp){
+             for(Int_t ii=0;(ii<n) && (!pp);ii++){
+               if(i==ii) continue; 
+               AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
+               Int_t iCaloCluster=tmp->GetBC();//index of AODCaloCluster
+                AliVCluster * clu = event->GetCaloCluster(iCaloCluster);
+               Int_t nCluPrimaries = clu->GetNLabels() ;
+               for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!pp); iAODLabel++){
+                 Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
+                 while(ipartnPrim>-1){
+                    if(ipartnPrim==ipartner)
+                     break ;
+                   ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
+                 }
+                 if(ipartnPrim==ipartner){
+                   pp=tmp ;
+                   break ;
+                 }
+               }
+             }
+           }
 
            if(pp){
              //Partner reconstructed, but did not pass cuts
                 FillPIDHistograms("hMCDecWRecPartn",p) ;       
                Double_t invMass=(*p+ *pp).M() ;
                FillHistogram("hMCmass",invMass,p->Pt(),p->GetWeight()) ;
-               if(IsInPi0Band(invMass,p->Pt())){
-                 FillPIDHistograms("hMCDecWithFoundPartn",p) ;
+               Double_t nSigma=InPi0Band(invMass,p->Pt()) ;
+               // analog to Tag
+                for(Int_t eminType=0; eminType<3; eminType++){
+                  if(pp->E()>0.1*(eminType+1)){
+                   for(Int_t isigma=0; isigma<3; isigma++){
+                     if(nSigma<1.+isigma){
+                        Int_t iType=3*eminType+isigma ;
+                        FillPIDHistograms(Form("hMCDecWithFoundPartnType%d",iType),p) ;
+                     }
+                   }
+                 }
                }
-               else{
+               if(nSigma>3.){
                  FillPIDHistograms("hMCDecWithWrongMass",p) ;
                }
            }
@@ -1046,6 +1084,28 @@ void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
                }
                if(!isPartnerLost){ //Reason not found!!!!!                               
                  FillPIDHistograms("hMCDecWMisPartnOther",p);
+                  Int_t multClust = event->GetNumberOfCaloClusters();
+                  for (Int_t iclu=0; (iclu<multClust) && (!isPartnerLost); iclu++) {
+                    AliVCluster * clu = event->GetCaloCluster(iclu);
+                    if(!clu->IsPHOS())
+                      continue ; 
+                   Int_t nCluPrimaries = clu->GetNLabels() ;
+                   for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!isPartnerLost); iAODLabel++){
+                     Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
+                     while(ipartnPrim>-1){
+                        if(ipartnPrim==ipartner)
+                         break ;
+                       ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
+                     }
+                     if(ipartnPrim==ipartner){
+                       isPartnerLost=kTRUE;
+                       break ;
+                     }
+                   }
+                 }
+                 if(isPartnerLost){//Did not pass default cuts                                   
+                   FillPIDHistograms("hMCDecWMisPartnDefCuts",p);
+                 }               
                }
                else{//Sum of all missed partners
                  FillPIDHistograms("hMCDecWMisPartnAll",p);
@@ -1086,8 +1146,8 @@ void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
             FillPIDHistograms("hInvM_Re_Emin3",p1,p2,invMass,kTRUE) ;
      
            //Fill izolated pi0s
-           Double_t nsigma1 = IsInPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
-           Double_t nsigma2 = IsInPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
+           Double_t nsigma1 = InPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
+           Double_t nsigma2 = InPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
            if(nsigma1<2 || nsigma2<2){ //2 sigma band
              TLorentzVector pi0=*p1+*p2 ;
               Int_t isolation=EvalIsolation(&pi0,0) ;
@@ -1157,12 +1217,12 @@ void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
 
       //Tagging: 1,2,3 sigma
       //Emin=100,200,300 Mev
-      //IsInPi0Band(mass, Ptphoton, type Emin cut
+      //InPi0Band(mass, Ptphoton, type Emin cut
       Int_t tag1=0 ;
       for(Int_t eminType=0; eminType<3; eminType++){
         if(p2->E()>0.1*(eminType+1)){
          //Set corresponding bit
-         Double_t nsigma = IsInPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
+         Double_t nsigma = InPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
          for(Int_t isigma=0; isigma<3; isigma++){
            if(nsigma<1+isigma){
              tag1|= (1 << (3*eminType+isigma)) ;
@@ -1177,7 +1237,7 @@ void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
       for(Int_t eminType=0; eminType<3; eminType++){
         if(p1->E()>0.1*(eminType+1)){
          //Set corresponding bit
-         Double_t nsigma = IsInPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
+         Double_t nsigma = InPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
          for(Int_t isigma=0; isigma<3; isigma++){
            if(nsigma<1+isigma){
              tag2|= (1 << (3*eminType+isigma)) ;
@@ -1349,7 +1409,7 @@ void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
   if (fDebug > 1) Printf("Terminate()");
 }
 //______________________________________________________________________________
-Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
+Double_t AliAnalysisTaskTaggedPhotons::InPi0Band(Double_t m, Double_t pt)const
 {
   //Parameterization of the pi0 peak region
   //for LHC13bcdef
@@ -1357,7 +1417,7 @@ Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
 
   Double_t mpi0sigma=TMath::Sqrt(5.22245e-03*5.22245e-03 +2.86851e-03*2.86851e-03/pt) + 9.09932e-05*pt ;
  
-  return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
+  return TMath::Abs(m-mpi0mean)/mpi0sigma ;
 }
 //______________________________________________________________________________
 Int_t AliAnalysisTaskTaggedPhotons::IsSameParent(const AliCaloPhoton *p1, const AliCaloPhoton *p2)const{
@@ -1531,7 +1591,7 @@ void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const Al
   
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p1,const AliCaloPhoton * p2,Double_t x, Bool_t isRe) const{
+void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p1,const AliCaloPhoton * p2,Double_t x, Bool_t /*isRe*/) const{
 
   Double_t ptPi = (*p1 + *p2).Pt() ;
   Double_t w=TMath::Sqrt(p1->GetWeight()*p2->GetWeight()) ;

diff --git a/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.h b/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.h
index 73905a6eb8cdeae2a27d6797ff2ed7472e1ed324..67c5e3aa830553139e63b5acf3fa81fbea8498e0 100644 (file)
--- a/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.h
+++ b/PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.h
@@ -49,7 +49,7 @@ protected:
   Int_t   GetFiducialArea(const Float_t * pos)const ; //what kind of fiducial area hit the photon
   Int_t   IsSameParent(const AliCaloPhoton *p1, const AliCaloPhoton *p2) const; //Check MC genealogy; return PDG of parent
   Bool_t  IsGoodChannel(Int_t mod, Int_t ix, Int_t iz) ;
-  Bool_t  IsInPi0Band(Double_t m, Double_t pt)const; //Check if invariant mass is within pi0 peak
+  Double_t  InPi0Band(Double_t m, Double_t pt)const; //Check if invariant mass is within pi0 peak
   Bool_t  TestDisp(Double_t l0, Double_t l1, Double_t e)const  ;
   Bool_t  TestTOF(Double_t /*t*/,Double_t /*en*/)const{return kTRUE;} 
   Bool_t  TestCharged(Double_t dr,Double_t en)const ;