- update on the aod filtering (fiorella)
authorjbook <jbook@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Mar 2013 13:31:05 +0000 (13:31 +0000)
committerjbook <jbook@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Mar 2013 13:31:05 +0000 (13:31 +0000)
- implementation of tpc ep and its correction (nanoAODs) and improvements

PWGDQ/dielectron/AliAnalysisTaskDielectronFilter.cxx
PWGDQ/dielectron/AliAnalysisTaskDielectronFilter.h
PWGDQ/dielectron/AliDielectron.cxx
PWGDQ/dielectron/AliDielectron.h
PWGDQ/dielectron/AliDielectronMC.cxx
PWGDQ/dielectron/AliDielectronMC.h
PWGDQ/dielectron/AliDielectronVarManager.h

index 79165ce..370e400 100644 (file)
@@ -54,7 +54,9 @@ fTriggerLogic(kAny),
 fTriggerAnalysis(0x0),
 fStoreLikeSign(kFALSE),
 fStoreRotatedPairs(kFALSE),
+fStoreEventsWithSingleTracks(kFALSE),
 fCreateNanoAOD(kFALSE),
+fStoreHeader(kFALSE),
 fEventFilter(0x0)
 {
   //
@@ -76,7 +78,9 @@ fTriggerLogic(kAny),
 fTriggerAnalysis(0x0),
 fStoreLikeSign(kFALSE),
 fStoreRotatedPairs(kFALSE),
+fStoreEventsWithSingleTracks(kFALSE),
 fCreateNanoAOD(kFALSE),
+fStoreHeader(kFALSE),
 fEventFilter(0x0)
 {
   //
@@ -256,9 +260,13 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *)
   else hasCand = (fDielectron->HasCandidates());
 
   if(fStoreRotatedPairs) hasCand = (hasCand || fDielectron->HasCandidatesTR());
+  if(fStoreEventsWithSingleTracks) hasCand = (hasCand || fDielectron->GetTrackArray(0) || fDielectron->GetTrackArray(1));
+
   
+  AliAODHandler *aodH=(AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+  AliAODExtension *extDielectron = aodH->GetFilteredAOD("AliAOD.Dielectron.root");
   if(hasCand){
-    AliAODHandler *aodH=(AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
     AliAODEvent *aod = aodH->GetAOD();
     
     // reset bit for all tracks
@@ -285,7 +293,6 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *)
       }
     }
     
-    AliAODExtension *extDielectron = aodH->GetFilteredAOD("AliAOD.Dielectron.root");
     extDielectron->SelectEvent();
     Int_t ncandidates=fDielectron->GetPairArray(1)->GetEntriesFast();
     if (ncandidates==1) fEventStat->Fill((kNbinsEvent));
@@ -314,8 +321,14 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *)
        nanoEv->AddVertex(tmpSpd);
        nanoEv->GetVertex(0)->SetNContributors((static_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()->GetNContributors());
        nanoEv->GetVertex(1)->SetNContributors((static_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertexSPD()->GetNContributors());
+       // set event plane 
+       nanoEv->GetHeader()->SetEventplane((static_cast<AliAODEvent*>(InputEvent()))->GetHeader()->GetEventplaneP());
+       nanoEv->GetHeader()->ResetEventplanePointer(); 
+       // set multiplicity
+       nanoEv->GetHeader()->SetRefMultiplicity((Int_t)values[AliDielectronVarManager::kNTrk]);
+       nanoEv->GetHeader()->SetRefMultiplicityPos((Int_t)values[AliDielectronVarManager::kNacc]);
+       //nanoEv->GetHeader()->SetRefMultiplicityNeg(values[AliDielectronVarManager::kMatchEffITSTPC]);
 
          for(int kj=0; kj<(fDielectron->GetTrackArray(0))->GetEntries(); kj++){
          Int_t posit = nanoEv->AddTrack((AliAODTrack*)fDielectron->GetTrackArray(0)->At(kj));
          Int_t posVtx = nanoEv->AddVertex(((AliAODTrack*)fDielectron->GetTrackArray(0)->At(kj))->GetProdVertex()); 
@@ -364,7 +377,15 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *)
     }
     if(isAOD) t->Fill();
   }
-  
+  if(fCreateNanoAOD && isAOD && (!hasCand) &&  fStoreHeader)  
+   {
+   // set event plane 
+   extDielectron->GetAOD()->GetHeader()->SetEventplane((static_cast<AliAODEvent*>(InputEvent()))->GetHeader()->GetEventplaneP());
+   extDielectron->GetAOD()->GetHeader()->ResetEventplanePointer();
+   extDielectron->GetTree()->Fill(); // fill header for all events without tracks
+   }
   PostData(1, const_cast<THashList*>(fDielectron->GetHistogramList()));
   PostData(2,fEventStat);
   return;
index fa8f1fe..2f64a89 100644 (file)
@@ -61,7 +61,9 @@ public:
 
   void SetStoreLikeSignCandidates(Bool_t storeLS) { fStoreLikeSign = storeLS; }
   void SetStoreRotatedPairs(Bool_t storeTR) { fStoreRotatedPairs = storeTR; }
-  void SetStoreTrackLegs(Bool_t storeTrackRef) { fCreateNanoAOD = storeTrackRef; }
+  void SetStoreEventsWithSingleTracks(Bool_t storeSingleTrk) { fStoreEventsWithSingleTracks = storeSingleTrk; }
+  void SetCreateNanoAODs(Bool_t storeTrackRef) { fCreateNanoAOD = storeTrackRef; }
+  void SetStoreHeader(Bool_t storeHeader) { fStoreHeader = storeHeader; }
 
   void SetEventFilter(AliAnalysisCuts * const filter) {fEventFilter=filter;}
 
@@ -84,7 +86,9 @@ private:
   
   Bool_t fStoreLikeSign;        // flag to store like-sign candidates
   Bool_t fStoreRotatedPairs;    // flag to store rotation
-  Bool_t fCreateNanoAOD;       // flag to store track legs
+  Bool_t fStoreEventsWithSingleTracks;    // flag to store events with a least one reconstructed track 
+  Bool_t fCreateNanoAOD;        // flag to create nanoAODs 
+  Bool_t fStoreHeader;          // flag to store header for all events 
 
   AliAnalysisCuts *fEventFilter;     // event filter
   
index 6740cfc..6c90c05 100644 (file)
@@ -49,6 +49,9 @@ The names are available via the function PairClassName(Int_t i)
 
 #include <AliKFParticle.h>
 
+#include <AliESDInputHandler.h>
+#include <AliAnalysisManager.h>
+#include <AliEPSelectionTask.h>
 #include <AliEventplane.h>
 #include <AliVEvent.h>
 #include <AliVParticle.h>
@@ -271,9 +274,8 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
   }
 
   // TPC event plane correction
-  AliEventplane *cevplane = new AliEventplane();
-  if (ev1 && cevplane && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0)) 
-    EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
+  if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0)) 
+    EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1);
 
   if (!fNoPairing){
     // create pairs and fill pair candidate arrays
@@ -296,7 +298,7 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
   //process event mixing
   if (fMixing) {
     fMixing->Fill(ev1,this);
-//     FillHistograms(0x0,kTRUE);
+    //     FillHistograms(0x0,kTRUE);
   }
 
   //in case there is a histogram manager, fill the QA histograms
@@ -307,7 +309,6 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
   // clear arrays
   if (!fDontClearArrays) ClearArrays();
   AliDielectronVarManager::SetTPCEventPlane(0x0);
-  delete cevplane;
 }
 
 //________________________________________________________________
@@ -479,13 +480,19 @@ void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
 
   //Fill track information, separately for the track array candidates
   if (!pairInfoOnly){
+    className2.Form("Track_%s",fgkPairClassNames[1]);  // unlike sign, SE only
     for (Int_t i=0; i<4; ++i){
       className.Form("Track_%s",fgkTrackClassNames[i]);
-      if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
+      Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
+      Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
+      if (!trkClass && !mergedtrkClass) continue;
       Int_t ntracks=fTracks[i].GetEntriesFast();
       for (Int_t itrack=0; itrack<ntracks; ++itrack){
         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
-        fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+        if(trkClass)
+         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+        if(mergedtrkClass && i<2)
+         fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
       }
     }
   }
@@ -598,7 +605,7 @@ void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
 }
 
 //________________________________________________________________
-void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
+void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
 {
   //
   // Prefilter tracks and tracks from pairs
@@ -606,186 +613,287 @@ void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTra
   // remove contribution of all tracks to the Q-vector that are in invariant mass window 
   //
 
-  // TODO: how should we implement the filtering for the nanoADOs
-
   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
-  if(!evplane) return;
-  
-  // do not change these vectors qref
-  TVector2 * const qref  = evplane->GetQVector();
-  if(!qref) return;
-  // random subevents
-  TVector2 *qrsub1 = evplane->GetQsub1();
-  TVector2 *qrsub2 = evplane->GetQsub2();
-
-  // copy references
-  TVector2 *qcorr  = new TVector2(*qref);
-  TVector2 *qcsub1 = 0x0;
-  TVector2 *qcsub2 = 0x0;
-  //  printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
-
-
-  // eta gap ?
-  Bool_t etagap = kFALSE;
-  for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
-    TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
-    if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
-  }
+  if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
+    //  if(1) {
+    // get the EPselectionTask for recalculation of weighting factors
+    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+    AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
+    if(!eptask) return;
+
+    // track mapping
+    TMap mapRemovedTracks;
+
+    // eta gap ?
+    Bool_t etagap = kFALSE;
+    for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
+      TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
+      if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
+    }
 
-  // subevent separation
-  if(fLikeSignSubEvents || etagap) {
-    qcsub1 = new TVector2(*qcorr);
-    qcsub2 = new TVector2(*qcorr);
+    Double_t cQX=0., cQY=0.;
+    // apply cuts to the tracks, e.g. etagap
+    if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
+      UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
+      Int_t ntracks=ev->GetNumberOfTracks();
+      for (Int_t itrack=0; itrack<ntracks; ++itrack){
+       AliVParticle *particle=ev->GetTrack(itrack);
+       AliVTrack *track= static_cast<AliVTrack*>(particle);
+       if (!track) continue;
+       //event plane cuts
+       UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
+       //apply cut
+       if (cutMask==selectedMask) continue;
+
+       mapRemovedTracks.Add(track,track);
+       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
+      }
+    }
 
-    Int_t ntracks=ev->GetNumberOfTracks();
-    
-    // track removals
-    for (Int_t itrack=0; itrack<ntracks; ++itrack){
-      AliVParticle *particle=ev->GetTrack(itrack);
-      AliVTrack *track= static_cast<AliVTrack*>(particle);
-      if (!track) continue;
+    // POI (particle of interest) rejection
+    Int_t pairIndex=GetPairIndex(arr1,arr2);
 
-      Double_t cQX     = evplane->GetQContributionX(track);
-      Double_t cQY     = evplane->GetQContributionY(track);
-      
-      // by charge sub1+ sub2-
-      if(fLikeSignSubEvents) {
-       Short_t charge=track->Charge();
-       if (charge<0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
-       if (charge>0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
-      }
-      // by eta sub1+ sub2-
-      if(etagap) {
-       Double_t eta=track->Eta();
-       if (eta<0.0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
-       if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
+    Int_t ntrack1=arrTracks1.GetEntriesFast();
+    Int_t ntrack2=arrTracks2.GetEntriesFast();
+    AliDielectronPair candidate;
+    candidate.SetKFUsage(fUseKF);
+
+    UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
+    for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
+      Int_t end=ntrack2;
+      if (arr1==arr2) end=itrack1;
+      Bool_t accepted=kFALSE;
+      for (Int_t itrack2=0; itrack2<end; ++itrack2){
+       TObject *track1=arrTracks1.UncheckedAt(itrack1);
+       TObject *track2=arrTracks2.UncheckedAt(itrack2);
+       if (!track1 || !track2) continue;
+       //create the pair
+       candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
+                           static_cast<AliVTrack*>(track2), fPdgLeg2);
+       candidate.SetType(pairIndex);
+       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
+
+       //event plane pair cuts
+       UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
+       //apply cut
+       if (cutMask==selectedMask) continue;
+
+       accepted=kTRUE;
+       //remove the tracks from the Track arrays
+       arrTracks2.AddAt(0x0,itrack2);
       }
+      if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
     }
-  }
-  else {
-    // by a random
-    qcsub1 = new TVector2(*qrsub1);
-    qcsub2 = new TVector2(*qrsub2);
-  }
-  
-  // apply cuts, e.g. etagap 
-  if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
-    UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
-    Int_t ntracks=ev->GetNumberOfTracks();
-    for (Int_t itrack=0; itrack<ntracks; ++itrack){
-      AliVParticle *particle=ev->GetTrack(itrack);
-      AliVTrack *track= static_cast<AliVTrack*>(particle);
+    //compress the track arrays
+    arrTracks1.Compress();
+    arrTracks2.Compress();
+
+    //Modify the components: subtract the tracks
+    ntrack1=arrTracks1.GetEntriesFast();
+    ntrack2=arrTracks2.GetEntriesFast();
+    // remove leg1 contribution
+    for (Int_t itrack=0; itrack<ntrack1; ++itrack){
+      AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
       if (!track) continue;
-      
-      //event plane cuts
-      UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
-      //apply cut
-      if (cutMask==selectedMask) continue;
+      // track contribution was already removed
+      if (mapRemovedTracks.FindObject(track)) continue;
+      else mapRemovedTracks.Add(track,track);
 
-      Double_t cQX     = 0.0;
-      Double_t cQY     = 0.0;
-      if(!etagap) {
-       cQX = evplane->GetQContributionX(track);
-       cQY = evplane->GetQContributionY(track);
-      }
-      Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
-      Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
-      Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
-      Double_t cQYsub2 = evplane->GetQContributionYsub2(track);      
-
-      // update Q vectors
-      qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
-      qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
-      qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
+      cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+      cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
     }
-  }
+    // remove leg2 contribution
+    for (Int_t itrack=0; itrack<ntrack2; ++itrack){
+      AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
+      if (!track) continue;
+      // track contribution was already removed
+      if (mapRemovedTracks.FindObject(track)) continue;
+      else mapRemovedTracks.Add(track,track);
 
-  // POI (particle of interest) rejection
-  Int_t pairIndex=GetPairIndex(arr1,arr2);
-  
-  Int_t ntrack1=arrTracks1.GetEntriesFast();
-  Int_t ntrack2=arrTracks2.GetEntriesFast();
-  AliDielectronPair candidate;
-  candidate.SetKFUsage(fUseKF);
-  
-  UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
-  for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
-    Int_t end=ntrack2;
-    if (arr1==arr2) end=itrack1;
-    Bool_t accepted=kFALSE;
-    for (Int_t itrack2=0; itrack2<end; ++itrack2){
-      TObject *track1=arrTracks1.UncheckedAt(itrack1);
-      TObject *track2=arrTracks2.UncheckedAt(itrack2);
-      if (!track1 || !track2) continue;
-      //create the pair
-      candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
-                          static_cast<AliVTrack*>(track2), fPdgLeg2);
-      
-      candidate.SetType(pairIndex);
-      candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
-      
-      //event plane cuts
-      UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
-      //apply cut
-      if (cutMask==selectedMask) continue;
+      cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+      cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
+    }
 
-      accepted=kTRUE;
-      //remove the tracks from the Track arrays
-      arrTracks2.AddAt(0x0,itrack2);
+    // build a corrected alieventplane using the values from the var manager
+    // these uncorrected values are filled using the stored magnitude and angle  in the header
+    TVector2 qcorr;
+    qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
+             AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
+    // fill alieventplane
+    AliEventplane cevplane;
+    cevplane.SetQVector(&qcorr);
+    AliDielectronVarManager::SetTPCEventPlane(&cevplane);
+    cevplane.SetQVector(0);
+    return;
+  } //end: nanoAODs
+  else
+    {
+    // this is done in case of ESDs or AODs
+    Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
+    // copy event plane object
+    AliEventplane cevplane(*evplane);
+    //    Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
+
+    TVector2 *qcorr  = cevplane.GetQVector();
+    if(!qcorr) return;
+    TVector2 *qcsub1 = 0x0;
+    TVector2 *qcsub2 = 0x0;
+
+    // eta gap ?
+    Bool_t etagap = kFALSE;
+    for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
+      TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
+      if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
     }
+
+    // subevent configuration for eta gap or LS (default is rndm)
+    if(fLikeSignSubEvents && etagap) {
+      // start with the full Qvector/event in both sub events
+      qcsub1 = new TVector2(*qcorr);
+      qcsub2 = new TVector2(*qcorr);
+      cevplane.SetQsub(qcsub1,qcsub2);
+
+      Int_t ntracks=ev->GetNumberOfTracks();
+      // track removals
+      for (Int_t itrack=0; itrack<ntracks; ++itrack){
+       AliVParticle *particle=ev->GetTrack(itrack);
+       AliVTrack *track= static_cast<AliVTrack*>(particle);
+       if (!track) continue;
+       if (track->GetID()>=0 && !isESD) continue;
+       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+
+       // set contributions to zero
+       // charge sub1+ sub2-
+       if(fLikeSignSubEvents) {
+         Short_t charge=track->Charge();
+         if (charge<0) {
+           cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+           cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+         }
+         if (charge>0) {
+           cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+           cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+         }
+       }
+       // eta sub1+ sub2-
+       if(etagap) {
+         Double_t eta=track->Eta();
+         if (eta<0.0) {
+           cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+           cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+         }
+         if (eta>0.0) {
+           cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+           cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+         }
+       }
+      } // end: loop over tracks
+    } // end: sub event configuration
+
+    // apply cuts, e.g. etagap
+    if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
+      UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
+      Int_t ntracks=ev->GetNumberOfTracks();
+      for (Int_t itrack=0; itrack<ntracks; ++itrack){
+       AliVParticle *particle=ev->GetTrack(itrack);
+       AliVTrack *track= static_cast<AliVTrack*>(particle);
+       if (!track) continue;
+       if (track->GetID()>=0 && !isESD) continue;
+       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+
+       //event plane cuts
+       UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
+       //apply cut
+       if (cutMask==selectedMask) continue;
+
+       // set contributions to zero
+       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+      }
+    } // end: track cuts
+
+    // POI (particle of interest) rejection
+    Int_t pairIndex=GetPairIndex(arr1,arr2);
+    Int_t ntrack1=arrTracks1.GetEntriesFast();
+    Int_t ntrack2=arrTracks2.GetEntriesFast();
+    AliDielectronPair candidate;
+    candidate.SetKFUsage(fUseKF);
+
+    UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
+    for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
+      Int_t end=ntrack2;
+      if (arr1==arr2) end=itrack1;
+      Bool_t accepted=kFALSE;
+      for (Int_t itrack2=0; itrack2<end; ++itrack2){
+       TObject *track1=arrTracks1.UncheckedAt(itrack1);
+       TObject *track2=arrTracks2.UncheckedAt(itrack2);
+       if (!track1 || !track2) continue;
+       //create the pair
+       candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
+                           static_cast<AliVTrack*>(track2), fPdgLeg2);
+
+       candidate.SetType(pairIndex);
+       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
+
+       //event plane cuts
+       UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
+       //apply cut
+       if (cutMask==selectedMask) continue;
+
+       accepted=kTRUE;
+       //remove the tracks from the Track arrays
+       arrTracks2.AddAt(0x0,itrack2);
+      }
       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
-  }
-  //compress the track arrays
-  arrTracks1.Compress();
-  arrTracks2.Compress();
-  
-  
-  //Modify the components: subtract the tracks
-  ntrack1=arrTracks1.GetEntriesFast();
-  ntrack2=arrTracks2.GetEntriesFast();
-  
-  // remove leg1 contribution
-  for (Int_t itrack=0; itrack<ntrack1; ++itrack){
-    AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
-    if (!track) continue;
-    
-    Double_t cQX     = evplane->GetQContributionX(track);
-    Double_t cQY     = evplane->GetQContributionY(track);
-    Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
-    Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
-    Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
-    Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
-    
-    // update Q vectors
-    qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
-    qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
-    qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
-  }
-  // remove leg2 contribution
-  for (Int_t itrack=0; itrack<ntrack2; ++itrack){
-    AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
-    if (!track) continue;
-    
-    Double_t cQX     = evplane->GetQContributionX(track);
-    Double_t cQY     = evplane->GetQContributionY(track);
-    Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
-    Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
-    Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
-    Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
-    
-    // update Q vectors
-    qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
-    qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
-    qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
-  }
+    }
+    //compress the track arrays
+    arrTracks1.Compress();
+    arrTracks2.Compress();
+
+    //Modify the components: subtract the tracks
+    ntrack1=arrTracks1.GetEntriesFast();
+    ntrack2=arrTracks2.GetEntriesFast();
+    // remove leg1 contribution
+    for (Int_t itrack=0; itrack<ntrack1; ++itrack){
+      AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
+      if (!track) continue;
+      if (track->GetID()>=0 && !isESD) continue;
+      Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+      // set contributions to zero
+      cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+    }
+    // remove leg2 contribution
+    for (Int_t itrack=0; itrack<ntrack2; ++itrack){
+      AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
+      if (!track) continue;
+      if (track->GetID()>=0 && !isESD) continue;
+      Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+      // set contributions to zero
+      cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+      cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+    }
 
-  //  printf("qrsub1 %p %f \t qcsub1 %p %f \n",qrsub1,qrsub1->X(),qcsub1,qcsub1->X());
-  // set AliEventplane with corrected values
-  cevplane->SetQVector(qcorr);
-  cevplane->SetQsub(qcsub1, qcsub2);
-  AliDielectronVarManager::SetTPCEventPlane(cevplane);
-}
+    // set corrected AliEventplane and fill variables with corrected values
+    AliDielectronVarManager::SetTPCEventPlane(&cevplane);
+    delete qcsub1;
+    delete qcsub2;
+  } // end: ESD or AOD case
 
+}
 //________________________________________________________________
 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
 {
index 829a9ed..d53c188 100644 (file)
@@ -176,7 +176,7 @@ private:
   Bool_t fDontClearArrays;      //Don't clear the arrays at the end of the Process function, needed for external use of pair and tracks
   
   void FillTrackArrays(AliVEvent * const ev, Int_t eventNr=0);
-  void EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane);
+  void EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev);
   void PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2);
   void FillPairArrays(Int_t arr1, Int_t arr2);
   void FillPairArrayTR();
index e97a61f..c2efcde 100644 (file)
@@ -30,6 +30,7 @@
 #include <AliMCEvent.h>
 #include <AliMCParticle.h>
 #include <AliAODMCParticle.h>
+#include <AliAODMCHeader.h>
 #include <AliStack.h>
 #include <AliESDEvent.h>
 #include <AliAODEvent.h>
@@ -196,10 +197,10 @@ Bool_t AliDielectronMC::ConnectMCEvent()
     if (!aodHandler) return kFALSE;
     AliAODEvent *aod=aodHandler->GetEvent();
     if (!aod) return kFALSE;
+
     fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
     if (!fMcArray){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
     else fHasMC=kTRUE;
-    fMCEvent=aodHandler->MCEvent();
   }
   return kTRUE;
 }
@@ -1282,3 +1283,23 @@ Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
   }
   return 0;
 }
+
+
+Bool_t AliDielectronMC::GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ){
+
+     if(fAnaType == kESD){
+     const AliVVertex* mcVtx =  fMCEvent->GetPrimaryVertex();
+     if(!mcVtx) return kFALSE;
+     primVtxX = mcVtx->GetX();
+     primVtxY = mcVtx->GetY();
+     primVtxZ = mcVtx->GetZ();
+     }else if(fAnaType == kAOD){
+     AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
+     AliAODMCHeader *mcHead = dynamic_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
+     primVtxX = mcHead->GetVtxX();
+     primVtxY = mcHead->GetVtxY();
+     primVtxZ = mcHead->GetVtxZ();
+     }
+
+return kTRUE;
+}
index 0fa61d1..cd616a7 100644 (file)
@@ -24,6 +24,7 @@ class AliAODTrack;
 class TParticle;
 class AliMCParticle;
 class AliAODMCParticle;
+class AliAODMCHeader;
 
 #include "AliDielectronSignalMC.h"
 #include "AliDielectronPair.h"
@@ -96,6 +97,7 @@ public:
   Int_t IsJpsiPrimary(const AliDielectronPair * pair);
   Int_t IsJpsiPrimary(const AliVParticle * pair);
   Bool_t CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const;
+  Bool_t GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ);
 
   AliMCEvent* GetMCEvent() { return fMCEvent; }         // return the AliMCEvent
   
@@ -123,7 +125,6 @@ private:
   Bool_t CheckIsRadiative(Int_t label) const;
   Bool_t CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const;
 
-
   ClassDef(AliDielectronMC, 0)
 };
 
index ec88988..ff822ad 100644 (file)
@@ -1055,8 +1055,10 @@ inline void AliDielectronVarManager::FillVarMCParticle2(const AliVParticle *p1,
   values[AliDielectronVarManager::kPseudoProperTime] = -2e10;
   if(mother) {    // same mother
     FillVarVParticle(mother, values);
-    Double_t lxy = ((mother->Xv()- values[AliDielectronVarManager::kXvPrim]) * mother->Px() + 
-                   (mother->Yv()- values[AliDielectronVarManager::kYvPrim]) * mother->Py() )/mother->Pt();
+    Double_t vtxX, vtxY, vtxZ;
+    mc->GetPrimaryVertex(vtxX,vtxY,vtxZ);
+    Double_t lxy = ((mother->Xv()- vtxX) * mother->Px() + 
+                   (mother->Yv()- vtxY) * mother->Py() )/mother->Pt();
     values[AliDielectronVarManager::kPseudoProperTime] = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/mother->Pt();
   }
   // AliVParticle part
@@ -1395,15 +1397,15 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa
     // values[AliDielectronVarManager::kPseudoProperTimePull] = -1e10;
     if(samemother && fgEvent) {
       if(pair->GetFirstDaughter()->GetLabel() > 0) {
-        const AliVParticle* d1 = mc->GetMCTrackFromMCEvent(pair->GetFirstDaughter()->GetLabel());
-        const AliVParticle* motherMC = mc->GetMCTrackFromMCEvent(((AliMCParticle*)d1)->GetMother());
-        const AliMCEvent *mcevent = mc->GetMCEvent();
-       const AliVVertex* mcVtx = mcevent ? mcevent->GetPrimaryVertex() : 0x0;
-       if(motherMC && mcVtx) {
+        const AliVParticle *motherMC = 0x0;
+        if(fgEvent->IsA() == AliESDEvent::Class())  motherMC = (AliMCParticle*)mc->GetMCTrackMother((AliESDtrack*)pair->GetFirstDaughter());
+        else if(fgEvent->IsA() == AliAODEvent::Class())  motherMC = (AliAODMCParticle*)mc->GetMCTrackMother((AliAODTrack*)pair->GetFirstDaughter());
+        Double_t vtxX, vtxY, vtxZ;
+       if(motherMC && mc->GetPrimaryVertex(vtxX,vtxY,vtxZ)) {
          Int_t motherLbl = motherMC->GetLabel();
          values[AliDielectronVarManager::kHasCocktailMother]=mc->CheckParticleSource(motherLbl, AliDielectronSignalMC::kDirect);
-         const Double_t lxyMC = ( (motherMC->Xv() - mcVtx->GetX()) * motherMC->Px() +
-                                   (motherMC->Yv() - mcVtx->GetY()) * motherMC->Py()   ) / motherMC->Pt();
+         const Double_t lxyMC = ( (motherMC->Xv() - vtxX) * motherMC->Px() +
+                                   (motherMC->Yv() - vtxY) * motherMC->Py()   ) / motherMC->Pt();
          const Double_t pseudoMC = lxyMC * (TDatabasePDG::Instance()->GetParticle(443)->Mass())/motherMC->Pt();
          values[AliDielectronVarManager::kPseudoProperTimeResolution] = values[AliDielectronVarManager::kPseudoProperTime] - pseudoMC;
           if (errPseudoProperTime2 > 0)
@@ -1790,12 +1792,9 @@ inline void AliDielectronVarManager::FillVarAODEvent(const AliAODEvent *event, D
   // nanoAODs (w/o AliEventPlane branch) should have the tpc event plane angle stored in the header
   if(!header->GetEventplaneP()) {
     values[AliDielectronVarManager::kTPCrpH2uc]  = header->GetEventplane();
-    //TODO: activate after new nano production
-    /*
     values[AliDielectronVarManager::kTPCmagH2uc] = header->GetEventplaneMag();
     values[AliDielectronVarManager::kTPCxH2uc]   = TMath::Cos(header->GetEventplane())*header->GetEventplaneMag();
     values[AliDielectronVarManager::kTPCyH2uc]   = TMath::Sin(header->GetEventplane())*header->GetEventplaneMag();
-    */
   }
 
   const AliAODVertex *vtxtpc = GetVertex(event, AliAODVertex::kMainTPC);
@@ -1834,13 +1833,13 @@ inline void AliDielectronVarManager::FillVarTPCEventPlane(const AliEventplane *e
     TVector2 *qcorr  = const_cast<AliEventplane *>(evplane)->GetQVector();  // This is the "corrected" Q-Vector
     TVector2 *qcsub1 = const_cast<AliEventplane *>(evplane)->GetQsub1();
     TVector2 *qcsub2 = const_cast<AliEventplane *>(evplane)->GetQsub2();
-    if(qcorr && qcsub1 && qcsub2) {
-
+    if(qcorr) {
       values[AliDielectronVarManager::kTPCxH2]   = qcorr->X();
       values[AliDielectronVarManager::kTPCyH2]   = qcorr->Y();
       values[AliDielectronVarManager::kTPCmagH2] = qcorr->Mod();
       values[AliDielectronVarManager::kTPCrpH2]  = ((TMath::Abs(qcorr->X())>1.0e-10) ? TMath::ATan2(qcorr->Y(),qcorr->X())/2.0 : 0.0);
-
+    }
+    if(qcsub1 && qcsub2) {
       values[AliDielectronVarManager::kTPCsub1xH2]   = qcsub1->X();
       values[AliDielectronVarManager::kTPCsub1yH2]   = qcsub1->Y();
       values[AliDielectronVarManager::kTPCsub1rpH2]  = ((TMath::Abs(qcsub1->X())>1.0e-10) ? TMath::ATan2(qcsub1->Y(),qcsub1->X())/2.0 : 0.0);
@@ -1853,7 +1852,6 @@ inline void AliDielectronVarManager::FillVarTPCEventPlane(const AliEventplane *e
                                                                         values[AliDielectronVarManager::kTPCsub2rpH2]) );
       values[AliDielectronVarManager::kTPCsub12DiffH2Sin] = TMath::Sin( 2.*(values[AliDielectronVarManager::kTPCsub1rpH2] -
                                                                            values[AliDielectronVarManager::kTPCsub2rpH2]) );
-
     }
   }
 }