]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/AliAnalysisTaskSEImproveITS.cxx
adding histogram of multiplicity vs centrality wo vertex selection (G.Luparello)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEImproveITS.cxx
index cad0518f4eb5d0a090d266e9bce7ea65026e4e50..5521ece56d41b8a2d18dff14d5d82686d6b4e2c7 100644 (file)
@@ -30,6 +30,8 @@
 #include "AliExternalTrackParam.h"
 #include "AliAODRecoDecayHF2Prong.h"
 #include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliNeutralTrackParam.h"
 #include "AliAnalysisTaskSEImproveITS.h"
 
 //
@@ -62,7 +64,26 @@ AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
    fPt1ResPUpg  (0),
    fPt1ResKUpg  (0),
    fPt1ResPiUpg (0),
+   fD0ZResPCurSA  (0),
+   fD0ZResKCurSA  (0),
+   fD0ZResPiCurSA (0),
+   fD0RPResPCurSA (0),
+   fD0RPResKCurSA (0),
+   fD0RPResPiCurSA(0),
+   fPt1ResPCurSA  (0),
+   fPt1ResKCurSA  (0),
+   fPt1ResPiCurSA (0),
+   fD0ZResPUpgSA  (0),
+   fD0ZResKUpgSA  (0),
+   fD0ZResPiUpgSA (0),
+   fD0RPResPUpgSA (0),
+   fD0RPResKUpgSA (0),
+   fD0RPResPiUpgSA(0),
+   fPt1ResPUpgSA  (0),
+   fPt1ResKUpgSA  (0),
+   fPt1ResPiUpgSA (0),
    fRunInVertexing(kFALSE),
+   fImproveTracks(kTRUE),
    fDebugOutput (0),
    fDebugNtuple (0),
    fDebugVars   (0), 
@@ -97,7 +118,26 @@ AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
    fPt1ResPUpg  (0),
    fPt1ResKUpg  (0),
    fPt1ResPiUpg (0),
+   fD0ZResPCurSA  (0),
+   fD0ZResKCurSA  (0),
+   fD0ZResPiCurSA (0),
+   fD0RPResPCurSA (0),
+   fD0RPResKCurSA (0),
+   fD0RPResPiCurSA(0),
+   fPt1ResPCurSA  (0),
+   fPt1ResKCurSA  (0),
+   fPt1ResPiCurSA (0),
+   fD0ZResPUpgSA  (0),
+   fD0ZResKUpgSA  (0),
+   fD0ZResPiUpgSA (0),
+   fD0RPResPUpgSA (0),
+   fD0RPResKUpgSA (0),
+   fD0RPResPiUpgSA(0),
+   fPt1ResPUpgSA  (0),
+   fPt1ResKUpgSA  (0),
+   fPt1ResPiUpgSA (0),
    fRunInVertexing(isRunInVertexing),
+   fImproveTracks(kTRUE),
    fDebugOutput (0),
    fDebugNtuple (0),
    fDebugVars   (0),
@@ -130,6 +170,24 @@ AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
   fPt1ResPCur  ->SetName("Pt1ResPCur"  );
   fPt1ResKCur  ->SetName("Pt1ResKCur"  );
   fPt1ResPiCur ->SetName("Pt1ResPiCur" );
+  fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));
+  fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));
+  fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));
+  fD0ZResPCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA"  )));
+  fD0ZResKCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA"  )));
+  fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));
+  fPt1ResPCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA"  )));
+  fPt1ResKCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA"  )));
+  fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));
+  fD0RPResPCurSA ->SetName("D0RPResPCurSA" );
+  fD0RPResKCurSA ->SetName("D0RPResKCurSA" );
+  fD0RPResPiCurSA->SetName("D0RPResPiCurSA");
+  fD0ZResPCurSA  ->SetName("D0ZResPCurSA"  ); 
+  fD0ZResKCurSA  ->SetName("D0ZResKCurSA"  );
+  fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );
+  fPt1ResPCurSA  ->SetName("Pt1ResPCurSA"  );
+  fPt1ResKCurSA  ->SetName("Pt1ResKCurSA"  );
+  fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );
   delete resfileCur;
   TFile *resfileUpg=TFile::Open(resfileUpgURI );
   fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
@@ -150,6 +208,24 @@ AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
   fPt1ResPUpg  ->SetName("Pt1ResPUpg"  );
   fPt1ResKUpg  ->SetName("Pt1ResKUpg"  );
   fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
+  fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));
+  fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));
+  fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));
+  fD0ZResPUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA"  )));
+  fD0ZResKUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA"  )));
+  fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));
+  fPt1ResPUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA"  )));
+  fPt1ResKUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA"  )));
+  fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));
+  fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );
+  fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );
+  fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");
+  fD0ZResPUpgSA  ->SetName("D0ZResPUpgSA"  );
+  fD0ZResKUpgSA  ->SetName("D0ZResKUpgSA"  );
+  fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );
+  fPt1ResPUpgSA  ->SetName("Pt1ResPUpgSA"  );
+  fPt1ResKUpgSA  ->SetName("Pt1ResKUpgSA"  );
+  fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );
   delete resfileUpg;
 
   DefineOutput(1,TNtuple::Class());
@@ -214,8 +290,11 @@ void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
   // Smear all tracks
   TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
   if (!mcs) return;
-  for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack)
-    SmearTrack(ev->GetTrack(itrack),mcs);
+  if (fImproveTracks) {
+    for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
+      SmearTrack(ev->GetTrack(itrack),mcs);
+    }
+  }
 
   // TODO: recalculated primary vertex
   AliVVertex *primaryVertex=ev->GetPrimaryVertex();
@@ -283,6 +362,56 @@ void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
   }
 
 
+  // Dstar->Kpipi
+  TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
+  
+  if (arrayCascade) {
+    for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
+      AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
+      //Get D0 from D*
+      AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();
+      
+      // recalculate vertices
+      //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
+
+      //soft pion
+      AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
+       
+      //track D0
+      AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
+
+      //!!!!TODO: covariance matrix
+      
+      // update d0 
+      Double_t d0z0[2],covd0z0[3];
+      Double_t d01[2],d01err[2];
+
+      //the D*
+      et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+      d01[0]=d0z0[0];
+      d01err[0] = TMath::Sqrt(covd0z0[0]); 
+      trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+      d01[1]=d0z0[0];
+      d01err[1] = TMath::Sqrt(covd0z0[0]);  
+      decayDstar->Setd0Prongs(2,d01);
+      decayDstar->Setd0errProngs(2,d01err);
+        
+      // delete v12;
+      delete trackD0; trackD0=NULL;
+
+       // a run for D*
+      Double_t px1[2],py1[2],pz1[2];
+      for (Int_t i=0;i<2;++i) {
+       const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
+       px1[i]=t1->Px();
+       py1[i]=t1->Py();
+       pz1[i]=t1->Pz();
+      }
+      decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
+      
+     }
+  }
+
 
   // Three prong
   TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
@@ -363,6 +492,11 @@ void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArr
   if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
     return;
 
+  // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
+  if (TESTBIT(track->GetITSClusterMap(),7)) return;
+  //
+
+
   // Get reconstructed track parameters
   AliExternalTrackParam et; et.CopyFromVTrack(track);
   Double_t *param=const_cast<Double_t*>(et.GetParameter());
@@ -400,28 +534,28 @@ void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArr
   Double_t spt1o =0.;
   switch (mc->GetPdgCode()) {
   case 2212: case -2212:
-    sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
-    sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
-    spt1o =EvalGraph(fPt1ResPCur ,ptmc);
-    sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
-    sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
-    spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
+    sd0rpo=EvalGraph(ptmc,fD0RPResPCur,fD0RPResPCurSA);
+    sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
+    spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
+    sd0rpn=EvalGraph(ptmc,fD0RPResPUpg,fD0RPResPUpgSA);
+    sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
+    spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
     break;
   case 321: case -321:
-    sd0rpo=EvalGraph(fD0RPResKCur,ptmc); 
-    sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
-    spt1o =EvalGraph(fPt1ResKCur ,ptmc);
-    sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
-    sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
-    spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
+    sd0rpo=EvalGraph(ptmc,fD0RPResKCur,fD0RPResKCurSA);
+    sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
+    spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
+    sd0rpn=EvalGraph(ptmc,fD0RPResKUpg,fD0RPResKUpgSA);
+    sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
+    spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
     break;
   case 211: case -211:
-    sd0rpo=EvalGraph(fD0RPResPiCur,ptmc); 
-    sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
-    spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
-    sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
-    sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
-    spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
+    sd0rpo=EvalGraph(ptmc,fD0RPResPiCur,fD0RPResPiCurSA);
+    sd0zo =EvalGraph(ptmc,fD0ZResPiCur,fD0ZResPiCurSA);
+    spt1o =EvalGraph(ptmc,fPt1ResPiCur,fPt1ResPiCurSA);
+    sd0rpn=EvalGraph(ptmc,fD0RPResPiUpg,fD0RPResPiUpgSA);
+    sd0zn =EvalGraph(ptmc,fD0ZResPiUpg,fD0ZResPiUpgSA);
+    spt1n =EvalGraph(ptmc,fPt1ResPiUpg,fPt1ResPiUpgSA);
     break;
   default:
     return;
@@ -461,6 +595,13 @@ void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArr
   track->SetPosition(x,kFALSE);
   track->SetP(p,kTRUE);
 
+
+  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
+  UChar_t itsClusterMap = track->GetITSClusterMap();
+  SETBIT(itsClusterMap,7);
+  track->SetITSClusterMap(itsClusterMap);
+  //
+
   // write out debug infos
   if (fDebugNtuple->GetEntries()<fNDebug) {
     Int_t idbg=0;
@@ -498,7 +639,7 @@ AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *o
   return vertex;
 }
 
-Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x) const {
+Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
   //
   // Evaluates a TGraph without linear extrapolation. Instead the last
   // valid point of the graph is used when out of range.
@@ -509,7 +650,12 @@ Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x)
   Int_t    n   =graph->GetN();
   Double_t xmin=graph->GetX()[0  ];
   Double_t xmax=graph->GetX()[n-1];
-  if (x<xmin) return graph->Eval(xmin);
+  if (x<xmin) {
+    if(!graphSA) return graph->Eval(xmin);
+    Double_t xminSA=graphSA->GetX()[0];
+    if(x<xminSA) return graphSA->Eval(xminSA);
+    return graphSA->Eval(x);
+  }
   if (x>xmax) return graph->Eval(xmax);
   return graph->Eval(x);
 }