]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/GammaConv/AliV0Reader.cxx
Added new background scheeme, did some cleanup. added new bethe block parameters...
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliV0Reader.cxx
index 6a99950091581129273ed1f7f06ee5b012c7590f..a9e829f63712908a3b2aee205e825cf24c094e90 100644 (file)
@@ -33,6 +33,7 @@
 #include "AliStack.h"
 #include "AliMCEventHandler.h"
 #include "AliTPCpidESD.h"
+#include "AliGammaConversionBGHandler.h"
 
 class iostream;
 class AliESDv0;
@@ -99,8 +100,13 @@ AliV0Reader::AliV0Reader() :
   fNSigmaMass(0.),
   fUseImprovedVertex(kFALSE),
   fUseOwnXYZCalculation(kFALSE),
-  fCurrentEventGoodV0s(),
-  fPreviousEventGoodV0s()
+  fDoCF(kFALSE),
+  fUseOnFlyV0Finder(kTRUE),
+  fUpdateV0AlreadyCalled(kFALSE),
+  fCurrentEventGoodV0s(NULL),
+//  fPreviousEventGoodV0s(),
+  fBGEventHandler(NULL),
+  fBGEventInitialized(kFALSE)
 {
   fTPCpid = new AliTPCpidESD;  
 }
@@ -161,8 +167,13 @@ AliV0Reader::AliV0Reader(const AliV0Reader & original) :
   fNSigmaMass(original.fNSigmaMass),
   fUseImprovedVertex(original.fUseImprovedVertex),
   fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
+  fDoCF(original.fDoCF),
+  fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
+  fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
-  fPreviousEventGoodV0s(original.fPreviousEventGoodV0s)
+  //  fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
+  fBGEventHandler(original.fBGEventHandler),
+  fBGEventInitialized(original.fBGEventInitialized)
 {
        
 }
@@ -182,7 +193,8 @@ AliV0Reader::~AliV0Reader()
 
 void AliV0Reader::Initialize(){
   //see header file for documentation
-       
+
+  fUpdateV0AlreadyCalled = kFALSE;     
   // Get the input handler from the manager
   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
   if(fESDHandler == NULL){
@@ -199,25 +211,72 @@ void AliV0Reader::Initialize(){
   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
   if(fMCTruth == NULL){
     //print warning here
+    fDoMC = kFALSE;
   }
        
   //Get pointer to the mc stack
-  fMCStack = fMCTruth->MCEvent()->Stack();
-  if(fMCStack == NULL){
-    //print warning here
+  if(fMCTruth){
+    fMCStack = fMCTruth->MCEvent()->Stack();
+    if(fMCStack == NULL){
+      //print warning here
+    }
+    // Better parameters for MonteCarlo from A. Kalweit 2010/01/8
+    fTPCpid->SetBetheBlochParameters( 2.15898e+00/50.,
+                                     1.75295e+01,
+                                     3.40030e-09,
+                                     1.96178e+00,
+                                     3.91720e+00);
+  }
+  else{
+    // Better parameters for data from A. Kalweit 2010/01/8
+    fTPCpid->SetBetheBlochParameters(0.0283086,
+                                    2.63394e+01,
+                                    5.04114e-11,
+                                    2.12543e+00,
+                                    4.88663e+00);
   }
-       
        
   // for CF
   //Get pointer to the mc event
-  fMCEvent = fMCTruth->MCEvent();
-  if(fMCEvent == NULL){
-    //print warning here
-  }    
-       
+  if(fDoCF && fDoMC){
+    fMCEvent = fMCTruth->MCEvent();
+    if(fMCEvent == NULL){
+      //print warning here
+      fDoCF = kFALSE;
+    }  
+  }
+
+
+
        
   AliKFParticle::SetField(fESDEvent->GetMagneticField());
-       
+
+  //  fCurrentEventGoodV0s = new TClonesArray("TClonesArray", 0);
+  fCurrentEventGoodV0s = new TClonesArray("AliKFParticle", 0);
+
+  if(fBGEventInitialized == kFALSE){
+    Double_t *zBinLimitsArray = new Double_t[8];//{-7,-5,-3,-1,1,3,5,7};
+    zBinLimitsArray[0] = -7;
+    zBinLimitsArray[1] = -5;
+    zBinLimitsArray[2] = -3;
+    zBinLimitsArray[3] = -1;
+    zBinLimitsArray[4] = 1;
+    zBinLimitsArray[5] = 3;
+    zBinLimitsArray[6] = 5;
+    zBinLimitsArray[7] = 7;
+    
+    Double_t *multiplicityBinLimitsArray= new Double_t[4];//={0,10,20,500};
+    multiplicityBinLimitsArray[0] = 0;
+    multiplicityBinLimitsArray[1] = 10;
+    multiplicityBinLimitsArray[2] = 20;
+    multiplicityBinLimitsArray[3] = 500;
+    
+    
+    fBGEventHandler = new AliGammaConversionBGHandler(8,4,10);
+    
+    fBGEventHandler->Initialize(zBinLimitsArray, multiplicityBinLimitsArray);
+    fBGEventInitialized = kTRUE;
+  }
 }
 
 AliESDv0* AliV0Reader::GetV0(Int_t index){
@@ -232,6 +291,22 @@ Bool_t AliV0Reader::CheckForPrimaryVertex(){
 }
 
 
+Bool_t AliV0Reader::CheckV0FinderStatus(Int_t index){
+
+  if(fUseOnFlyV0Finder){
+    if(!GetV0(index)->GetOnFlyStatus()){
+      return kFALSE;
+    }
+  }
+  if(!fUseOnFlyV0Finder){
+    if(!GetV0(index)->GetOnFlyStatus()){
+      return kFALSE;
+    }
+  }
+  return kTRUE;
+}
+
+
 
 Bool_t AliV0Reader::NextV0(){
   //see header file for documentation
@@ -245,13 +320,15 @@ Bool_t AliV0Reader::NextV0(){
       fCurrentV0IndexNumber++;
       continue;
     }
-
     Double_t containerInput[3];
-    containerInput[0] = GetMotherCandidatePt();
-    containerInput[1] = GetMotherCandidateEta();
-    containerInput[2] = GetMotherCandidateMass();
-               
-
+    if(fDoCF){
+      containerInput[0] = GetMotherCandidatePt();
+      containerInput[1] = GetMotherCandidateEta();
+      containerInput[2] = GetMotherCandidateMass();
+    }
+    
+   
     //checks if on the fly mode is set
     if ( !fCurrentV0->GetOnFlyStatus() ){
       if(fHistograms != NULL){
@@ -260,8 +337,10 @@ Bool_t AliV0Reader::NextV0(){
       fCurrentV0IndexNumber++;
       continue;
     }
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly);            // for CF       
-               
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly);          // for CF       
+    }
+    
     //checks if we have a prim vertex
     if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) { 
       if(fHistograms != NULL){
@@ -270,8 +349,9 @@ Bool_t AliV0Reader::NextV0(){
       fCurrentV0IndexNumber++;
       continue;
     }
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors);               // for CF       
-               
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors);             // for CF       
+    }
                
     //Check the pid probability
     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
@@ -281,21 +361,8 @@ Bool_t AliV0Reader::NextV0(){
       fCurrentV0IndexNumber++;
       continue;
     }
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID);                      // for CF
-               
-               
-               
-    if(fUseOwnXYZCalculation == kFALSE){
-      fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
-    }
-    else{
-      Double_t convpos[2];
-      convpos[0]=0;
-      convpos[1]=0;
-      GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
-      fCurrentXValue = convpos[0];
-      fCurrentYValue = convpos[1];
-      fCurrentZValue = GetConvPosZ(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField());
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID);                    // for CF
     }
                
     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
@@ -304,9 +371,10 @@ Bool_t AliV0Reader::NextV0(){
       }
       fCurrentV0IndexNumber++;
       continue;
-    }          
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepR);                   // for CF
-               
+    }  
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepR);                 // for CF
+    }
                
                
     if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ){ // cuts out regions where we do not reconstruct
@@ -315,9 +383,10 @@ Bool_t AliV0Reader::NextV0(){
       }
       fCurrentV0IndexNumber++;
       continue;
-    }          
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine);                        // for CF
-               
+    }
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine);                      // for CF
+    }
                
     if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
       if(fHistograms != NULL){
@@ -325,9 +394,10 @@ Bool_t AliV0Reader::NextV0(){
       }
       fCurrentV0IndexNumber++;
       continue;
-    }          
-    fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ);           // for CF       
-               
+    }
+    if(fDoCF){
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ);         // for CF       
+    }
                
     /* Moved further up so corr framework can work
        if(UpdateV0Information() == kFALSE){
@@ -345,8 +415,9 @@ Bool_t AliV0Reader::NextV0(){
        fCurrentV0IndexNumber++;
        continue;
       }
-      fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF);               // for CF       
-                       
+      if(fDoCF){
+       fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF);              // for CF       
+      }
                        
       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
       if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
@@ -356,8 +427,9 @@ Bool_t AliV0Reader::NextV0(){
        fCurrentV0IndexNumber++;
        continue;
       }
-      fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2);                      // for CF
-                       
+      if(fDoCF){
+       fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2);                     // for CF
+      }
                        
       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){
        if(fHistograms != NULL){
@@ -366,8 +438,9 @@ Bool_t AliV0Reader::NextV0(){
        fCurrentV0IndexNumber++;
        continue;
       }
-      fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta);                       // for CF
-                       
+      if(fDoCF){
+       fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta);                      // for CF
+      }
                        
       if(fMotherCandidateLorentzVector->Pt()<fPtCut){
        if(fHistograms != NULL){
@@ -376,8 +449,9 @@ Bool_t AliV0Reader::NextV0(){
        fCurrentV0IndexNumber++;
        continue;
       }
-      fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt);                        // for CF
-                       
+      if(fDoCF){
+       fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt);                       // for CF
+      }
                        
     }
     else if(fUseESDTrack){
@@ -388,8 +462,10 @@ Bool_t AliV0Reader::NextV0(){
       fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
     }
 
-    fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
-               
+    //    fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
+
+    new((*fCurrentEventGoodV0s)[fCurrentEventGoodV0s->GetEntriesFast()])  AliKFParticle(*fCurrentMotherKFCandidate);
+
     iResult=kTRUE;//means we have a v0 who survived all the cuts applied
                
     fCurrentV0IndexNumber++;
@@ -411,7 +487,7 @@ Bool_t AliV0Reader::UpdateV0Information(){
        
   if(fCurrentNegativeESDTrack->GetSign() == fCurrentPositiveESDTrack->GetSign()){             // avoid like sign
     iResult=kFALSE;
-    if(fHistograms != NULL){
+    if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
       fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
     }
   }
@@ -426,9 +502,8 @@ Bool_t AliV0Reader::UpdateV0Information(){
       !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
     //  if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) || 
     //      !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
-               
     iResult=kFALSE;
-    if(fHistograms != NULL){
+    if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
       fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
     }
   }
@@ -437,7 +512,7 @@ Bool_t AliV0Reader::UpdateV0Information(){
       fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {                 
                
     iResult=kFALSE;
-    if(fHistograms != NULL){
+    if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
       fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
     }
   }
@@ -449,7 +524,7 @@ Bool_t AliV0Reader::UpdateV0Information(){
        fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
        fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
       iResult=kFALSE;
-      if(fHistograms != NULL){
+      if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
        fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
       }
     }
@@ -458,7 +533,7 @@ Bool_t AliV0Reader::UpdateV0Information(){
         fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
         fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
        iResult=kFALSE;
-       if(fHistograms != NULL){
+       if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
          fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
        }
       }
@@ -469,15 +544,12 @@ Bool_t AliV0Reader::UpdateV0Information(){
         fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
         fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
        iResult=kFALSE;
-       if(fHistograms != NULL){
+       if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
          fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
        }
       }
     }
   }
-
-
-
        
   if(fCurrentNegativeKFParticle != NULL){
     delete fCurrentNegativeKFParticle;
@@ -562,20 +634,35 @@ Bool_t AliV0Reader::UpdateV0Information(){
     }
   }
        
-  //  if(iResult==kTRUE){
-  //   fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate); // moved it to NextV0() after all the cuts are applied
-  //  }
-
 
   // for CF
   Double_t containerInput[3];
-  containerInput[0] = GetMotherCandidatePt();
-  containerInput[1] = GetMotherCandidateEta();
-  containerInput[2] = GetMotherCandidateMass();
+  if(fDoCF){
+    containerInput[0] = GetMotherCandidatePt();
+    containerInput[1] = GetMotherCandidateEta();
+    containerInput[2] = GetMotherCandidateMass();
+    
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);            // for CF       
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);            // for CF       
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);               // for CF       
+  }
+  
+
+  if(fUseOwnXYZCalculation == kFALSE){
+    fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
+  }
+  else{
+    Double_t convpos[2];
+    convpos[0]=0;
+    convpos[1]=0;
+    GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
+    fCurrentXValue = convpos[0];
+    fCurrentYValue = convpos[1];
+    fCurrentZValue = GetConvPosZ(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField());
+  }
 
-  fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);              // for CF       
-  fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);              // for CF       
-  fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);         // for CF       
+
+  fUpdateV0AlreadyCalled = kTRUE;
 
   return iResult;
 }
@@ -643,30 +730,10 @@ void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
 
 void AliV0Reader::UpdateEventByEventData(){
   //see header file for documentation
-       
-  if(fCurrentEventGoodV0s.size() >0 ){
-    //    fPreviousEventGoodV0s.clear();
-    //    fPreviousEventGoodV0s = fCurrentEventGoodV0s;
-    if(fPreviousEventGoodV0s.size()>19){
-      for(UInt_t nCurrent=0;nCurrent<fCurrentEventGoodV0s.size();nCurrent++){
-       fPreviousEventGoodV0s.erase(fPreviousEventGoodV0s.begin());
-       fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
-      }
-    }
-    else{
-      for(UInt_t nCurrent=0;nCurrent<fCurrentEventGoodV0s.size();nCurrent++){
-       if(fPreviousEventGoodV0s.size()<20){
-         fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
-       }
-       else{
-         fPreviousEventGoodV0s.erase(fPreviousEventGoodV0s.begin());
-         fPreviousEventGoodV0s.push_back(fCurrentEventGoodV0s.at(nCurrent));
-       }
-      }
-    }
+  if(fCurrentEventGoodV0s->GetEntriesFast() >0 ){
+    fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetZ(),fESDEvent->GetNumberOfTracks());
   }
-  fCurrentEventGoodV0s.clear();
-       
+  fCurrentEventGoodV0s->Delete();
   fCurrentV0IndexNumber=0;
 }
 
@@ -941,3 +1008,8 @@ Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Doubl
 
    return convposz;
 }
+
+AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t event){
+
+  return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fESDEvent->GetNumberOfTracks());
+}