protection for ProdVtxRange2D added

[u/mrichter/AliRoot.git] / CORRFW / AliCFParticleGenCuts.cxx

diff --git a/CORRFW/AliCFParticleGenCuts.cxx b/CORRFW/AliCFParticleGenCuts.cxx
index 8c5b69933fb7c06ded76dc0d32333740baa70910..b66f4e7407292c2520174001ffdaaff914acd76c 100644 (file)
--- a/CORRFW/AliCFParticleGenCuts.cxx
+++ b/CORRFW/AliCFParticleGenCuts.cxx
@@ -50,6 +50,8 @@ AliCFParticleGenCuts::AliCFParticleGenCuts() :
   fRequireIsPrimary(0),
   fRequireIsSecondary(0),
   fRequirePdgCode(0),
+  fRequireAbsolutePdg(0),
+  fProdVtxRange2D(0),
   fPdgCode(0),
   fProdVtxXMin (-1.e+09),
   fProdVtxYMin (-1.e+09),
@@ -79,6 +81,9 @@ AliCFParticleGenCuts::AliCFParticleGenCuts() :
   for (int i=0; i<kNCuts; i++) 
     for (int j=0; j<kNStepQA; j++) 
       fhQA[i][j]=0x0;
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=0x0;
 }
 
 //______________________________
@@ -91,6 +96,8 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* tit
   fRequireIsPrimary(0),
   fRequireIsSecondary(0),
   fRequirePdgCode(0),
+  fRequireAbsolutePdg(0),
+  fProdVtxRange2D(0),
   fPdgCode(0),
   fProdVtxXMin (-1.e+09),
   fProdVtxYMin (-1.e+09),
@@ -120,6 +127,9 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* tit
   for (int i=0; i<kNCuts; i++) 
     for (int j=0; j<kNStepQA; j++) 
       fhQA[i][j]=0x0;
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=0x0;
 }
 
 //______________________________
@@ -132,6 +142,8 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
   fRequireIsPrimary(c.fRequireIsPrimary),
   fRequireIsSecondary(c.fRequireIsSecondary),
   fRequirePdgCode(c.fRequirePdgCode),
+  fRequireAbsolutePdg(c.fRequireAbsolutePdg),
+  fProdVtxRange2D(c.fProdVtxRange2D),
   fPdgCode(c.fPdgCode),
   fProdVtxXMin (c.fProdVtxXMin),
   fProdVtxYMin (c.fProdVtxYMin),
@@ -161,6 +173,9 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
   for (int i=0; i<kNCuts; i++)
     for (int j=0; j<kNStepQA; j++)
       fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
 }
 
 //______________________________
@@ -178,6 +193,8 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts
     fRequireIsPrimary=c.fRequireIsPrimary;
     fRequireIsSecondary=c.fRequireIsSecondary;
     fRequirePdgCode=c.fRequirePdgCode;
+    fRequireAbsolutePdg=c.fRequireAbsolutePdg;
+    fProdVtxRange2D=c.fProdVtxRange2D;
     fPdgCode=c.fPdgCode;
     fProdVtxXMin=c.fProdVtxXMin;
     fProdVtxYMin=c.fProdVtxYMin;
@@ -205,6 +222,9 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts
     for (int i=0; i<kNCuts; i++)
       for (int j=0; j<kNStepQA; j++)
        fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
+
+    for (int j=0; j<kNStepQA; j++)
+       fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
   }
   return *this ;
 }
@@ -248,6 +268,14 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
   Double32_t partVy=(Double32_t)mcPart->Yv();
   Double32_t partVz=(Double32_t)mcPart->Zv();
 
+  // calculate the production vertex ellipse
+  Double32_t prodVtxXYmin = 0.;
+  if (fProdVtxXMin>0 && fProdVtxYMin>0)
+       prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
+  Double32_t prodVtxXYmax = 0.;
+  if(fProdVtxXMax>0 && fProdVtxYMax>0)
+       prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
+
   Double32_t decayVx=0.;
   Double32_t decayVy=0.;
   Double32_t decayVz=0.;
@@ -301,7 +329,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
   
   // cut on PDG code
   if ( fRequirePdgCode ) {
-    if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+    if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
   }
   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
   
@@ -310,6 +338,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
     Bool_t goodDecay = kTRUE ;
     Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
+    //now number of daughters is OK
     if (goodDecay) {
       // now check if decay channel is respected
       // first try
@@ -325,9 +354,26 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
          if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
        }
       }
-    }
+      if (!goodDecay && fRequireAbsolutePdg) {
+       //now tries inverting the sign of the daughters in case the anti-particle is also looked at
+       // third try
+       goodDecay = kTRUE ;
+       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
+         if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+       }
+       if (!goodDecay) {
+         //fourth try inverting the order of the daughters
+         goodDecay = kTRUE ;
+         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
+           if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+         }
+       }
+      } //end check anti-particle
+    } //end # daughters OK
     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
-  }
+  } //end require decay channel
   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
   
   
@@ -336,10 +382,31 @@ void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
+    || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
+   fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
+    || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
+    || ( fProdVtxRange2D &&  (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
+    || ( fProdVtxRange2D &&  (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
+    || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
+
   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
@@ -373,6 +440,14 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
   Double32_t partVy=(Double32_t)mcPart->Yv();
   Double32_t partVz=(Double32_t)mcPart->Zv();
 
+  // calculate the production vertex ellipse
+  Double32_t prodVtxXYmin = 0.;
+  if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
+       prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
+  Double32_t prodVtxXYmax = 0.;
+  if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
+       prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
+
   Double32_t decayVx=0.;
   Double32_t decayVy=0.;
   Double32_t decayVz=0.;
@@ -429,7 +504,7 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
   
   // cut on PDG code
   if ( fRequirePdgCode ) {
-    if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+    if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
   }
   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
   
@@ -466,10 +541,23 @@ void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
-  if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
+    || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
+    || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
+    || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
+    || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
@@ -496,6 +584,8 @@ void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) 
     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
 
+  fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
+
   // fill cut statistics and cut correlation histograms with information from the bitmap
   if (afterCuts) return;
 
@@ -524,6 +614,7 @@ void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
   qaList->Add(fhCutCorrelation);
 
   for (Int_t j=0; j<kNStepQA; j++) {
+    qaList->Add(fhProdVtxXY[j]);
     for(Int_t i=0; i<kNCuts; i++)
       qaList->Add(fhQA[i][j]);
   }
@@ -576,12 +667,12 @@ void AliCFParticleGenCuts::DefineHistograms() {
     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
-    fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10);
-    fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10);
-    fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10);
-    fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10);
-    fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10);
-    fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10);
+    fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
@@ -593,6 +684,11 @@ void AliCFParticleGenCuts::DefineHistograms() {
     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
+    fhProdVtxXY                  [i] = new TH2F(Form("%s_prodVtxXY%s"   ,GetName(),str),"",100,0,10,100,0,10);
+    fhProdVtxXY                  [i] ->GetXaxis()->SetTitle("x_{production vertex}");
+    fhProdVtxXY                  [i] ->GetYaxis()->SetTitle("y_{production vertex}");
+    fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
+    fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
   }
   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
 }
@@ -644,7 +740,7 @@ Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
   //
   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
-  //absolute value. By default is set to kFALSE.
+  //absolute value. 
   //
   TParticle* part = mcPart->Particle();
   Int_t pdgCode = part->GetPdgCode();
@@ -660,7 +756,7 @@ Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
   //
   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
-  //absolute value. By default is set to kFALSE.
+  //absolute value. 
   //
   Int_t pdgCode = mcPart->GetPdgCode();
   
@@ -672,7 +768,7 @@ Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs
   return kTRUE;
 }
 //______________________________
-void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
+void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
   //
   // Sets pointer to MC event information (AliMCEvent)
   //
@@ -687,7 +783,6 @@ void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
     AliError("argument must point to an AliMCEvent or an AliAODEvent !");
     return ;
   }
-  
-  if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ;
-  else          fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ;
+
+  fMCInfo = (AliVEvent*)mcEvent ;
 }