adding delta phi study
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 Mar 2008 10:31:39 +0000 (10:31 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 Mar 2008 10:31:39 +0000 (10:31 +0000)
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h
PWG0/dNdEta/run.C
PWG0/dNdEta/runCorrection.C

index 70b5674..8c25588 100644 (file)
@@ -40,6 +40,7 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fOption(opt),
   fAnalysisMode(AliPWG0Helper::kTPC),
   fSignMode(0),
+  fOnlyPrimaries(kFALSE),
   fEsdTrackCuts(0),
   fdNdEtaCorrection(0),
   fdNdEtaAnalysisMC(0),
@@ -49,6 +50,8 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fVertexCorrelation(0),
   fVertexProfile(0),
   fVertexShiftNorm(0),
+  fEtaCorrelation(0),
+  fEtaProfile(0),
   fSigmaVertexTracks(0),
   fSigmaVertexPrim(0),
   fMultAll(0),
@@ -66,7 +69,7 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   for (Int_t i=0; i<2; i++)
     fdNdEtaCorrectionProcessType[i] = 0;
   
-  for (Int_t i=0; i<3; i++)
+  for (Int_t i=0; i<8; i++)
     fDeltaPhi[i] = 0;
 }
 
@@ -180,12 +183,15 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
   fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
   
+  fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
+  fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
+
   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
   fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
 
-  for (Int_t i=0; i<3; i++)
-    fDeltaPhi[i] = new TH1F(Form("fDeltaPhi_%d", i), ";#Delta phi;Entries", 200, -0.5, 0.5);
+  for (Int_t i=0; i<8; i++)
+    fDeltaPhi[i] = new TH1F(Form("fDeltaPhi_%d", i), ";#Delta phi;Entries", 2000, -0.1, 0.1);
 }
 
 void AlidNdEtaCorrectionTask::Exec(Option_t*)
@@ -199,9 +205,15 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     return;
   }
 
+  if (fOnlyPrimaries)
+    Printf("WARNING: Processing only primaries. For systematical studies only!");
+
   // trigger definition
   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
 
+  if (!eventTriggered)
+    Printf("No trigger");
+
   // post the data already here
   PostData(0, fOutput);
 
@@ -261,10 +273,10 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   else
     Printf("No vertex found");
 
-
   // create list of (label, eta, pt) tuples
   Int_t inputCount = 0;
   Int_t* labelArr = 0;
+  Int_t* labelArr2 = 0; // only for case of SPD
   Float_t* etaArr = 0;
   Float_t* ptArr = 0;
   Float_t* deltaPhiArr = 0;
@@ -279,6 +291,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     }
 
     labelArr = new Int_t[mult->GetNumberOfTracklets()];
+    labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
     etaArr = new Float_t[mult->GetNumberOfTracklets()];
     ptArr = new Float_t[mult->GetNumberOfTracklets()];
     deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
@@ -298,8 +311,13 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       if (TMath::Abs(deltaPhi) > 1)
         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 
+      if (fOnlyPrimaries)
+        if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+          continue;
+
       etaArr[inputCount] = mult->GetEta(i);
-      labelArr[inputCount] = mult->GetLabel(i);
+      labelArr[inputCount] = mult->GetLabel(i, 0);
+      labelArr2[inputCount] = mult->GetLabel(i, 1);
       ptArr[inputCount] = 0; // no pt for tracklets
       deltaPhiArr[inputCount] = deltaPhi;
       ++inputCount;
@@ -320,6 +338,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     Printf("Accepted %d tracks", nGoodTracks);
 
     labelArr = new Int_t[nGoodTracks];
+    labelArr2 = new Int_t[nGoodTracks];
     etaArr = new Float_t[nGoodTracks];
     ptArr = new Float_t[nGoodTracks];
     deltaPhiArr = new Float_t[nGoodTracks];
@@ -336,6 +355,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
       etaArr[inputCount] = esdTrack->Eta();
       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+      labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
       ptArr[inputCount] = esdTrack->Pt();
       deltaPhiArr[inputCount] = 0; // no delta phi for tracks
       ++inputCount;
@@ -406,14 +426,16 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       fMultVtx->Fill(nAccepted);
   }
 
+  Int_t processed = 0;
+
   for (Int_t i=0; i<inputCount; ++i)
   {
     Int_t label = labelArr[i];
+    Int_t label2 = labelArr2[i];
 
     if (label < 0)
     {
       Printf("WARNING: cannot find corresponding mc part for track(let) %d with label %d.", i, label);
-      fDeltaPhi[2]->Fill(deltaPhiArr[i]);
       continue;
     }
 
@@ -423,14 +445,53 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
       continue;
     }
-
-    if (SignOK(particle->GetPDG()) == kFALSE)
-        continue;
+    
+    // find particle that is filled in the correction map
+    // this should be the particle that has been reconstructed
+    // for TPC this is clear and always identified by the track's label
+    // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken: 
+    //  L1 L2 (P = primary, S = secondary)
+    //   P P' : --> P
+    //   P S  : --> P
+    //   S P  : --> P
+    //   S S' : --> S
+
+    if (label != label2 && label2 >= 0)
+    {
+      if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
+      {
+        particle = stack->Particle(label2);
+        if (!particle)
+        {
+          AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
+          continue;
+        }
+      }
+    }
 
     if (eventTriggered && eventVertex)
     {
-      fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+      if (SignOK(particle->GetPDG()) == kFALSE)
+        continue;
+
+      processed++;
+
+      Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
+      // in case of primary the MC values are filled, otherwise (background) the reconstructed values
+      if (label == label2 && firstIsPrim)
+      {
+        fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+      }
+      else
+      {
+        fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], etaArr[i], ptArr[i]);
+        fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
+      }
+             
       fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+
+      fEtaCorrelation->Fill(etaArr[i], particle->Eta());
+
       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
       {
         fPIDTracks->Fill(particle->GetPdgCode());
@@ -451,15 +512,65 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
       }
 
-      if (stack->IsPhysicalPrimary(label))
+      Int_t hist = -1;
+      if (label == label2)
       {
-        fDeltaPhi[0]->Fill(deltaPhiArr[i]);
+       if (firstIsPrim)
+        {
+         hist = 0;
+        }
+        else
+          hist = 1; 
+      }
+      else if (label2 >= 0)
+      {
+        Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
+        if (firstIsPrim && secondIsPrim)
+        {
+          hist = 2;
+        }
+        else if (firstIsPrim && !secondIsPrim)
+        {
+          hist = 3;
+         
+          // check if secondary is caused by the primary or it is a fake combination
+          //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
+          Int_t mother = label2;
+          while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
+          {
+            //Printf("  %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
+            if (stack->Particle(mother)->GetMother(0) == label)
+            {
+              /*Printf("The untraceable decay was:");
+              Printf("   from:");
+              particle->Print();
+              Printf("   to (status code %d):", stack->Particle(mother)->GetStatusCode());
+              stack->Particle(mother)->Print();*/
+              hist = 4;
+            }
+            mother = stack->Particle(mother)->GetMother(0);
+          }
+        }
+        else if (!firstIsPrim && secondIsPrim)
+        {
+          hist = 5;
+        }
+        else if (!firstIsPrim && !secondIsPrim)
+        {
+          hist = 6;
+        }
+
       }
       else
-        fDeltaPhi[1]->Fill(deltaPhiArr[i]);
+        hist = 7;
+      
+      fDeltaPhi[hist]->Fill(deltaPhiArr[i]);
     }
   }
 
+  if (processed < inputCount)
+    Printf("Only %d out of %d track(let)s were used", processed, inputCount); 
+
   if (eventTriggered && eventVertex)
   {
     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
@@ -486,6 +597,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
   delete[] etaArr;
   delete[] labelArr;
+  delete[] labelArr2;
   delete[] ptArr;
   delete[] deltaPhiArr;
 
@@ -591,6 +703,12 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
     fVertexProfile->Write();
   if (fVertexShiftNorm)
     fVertexShiftNorm->Write();
+
+  if (fEtaCorrelation)
+    fEtaCorrelation->Write();
+  if (fEtaProfile)
+    fEtaProfile->Write();
+
   if (fMultAll)
     fMultAll->Write();
 
@@ -600,7 +718,7 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   if (fMultVtx)
     fMultVtx->Write();
 
-  for (Int_t i=0; i<3; ++i)
+  for (Int_t i=0; i<8; ++i)
     if (fDeltaPhi[i])
       fDeltaPhi[i]->Write();
 
index 42d067f..1fec160 100644 (file)
@@ -28,6 +28,7 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
 
     void SetTrackCuts(AliESDtrackCuts* cuts) { fEsdTrackCuts = cuts; }
     void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
+    void SetOnlyPrimaries(Bool_t flag = kTRUE) { fOnlyPrimaries = flag; }
 
  protected:
     Bool_t SignOK(TParticlePDG* particle);
@@ -38,6 +39,7 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     TString fOption;                 // option string
     AliPWG0Helper::AnalysisMode fAnalysisMode;    // detector that is used for analysis
     Int_t fSignMode;                 // if 0 process all particles, if +-1 process only particles with that sign
+    Bool_t fOnlyPrimaries;           // only process primaries (syst. studies)
 
     AliESDtrackCuts*  fEsdTrackCuts;             // Object containing the parameters of the esd track cuts
 
@@ -54,6 +56,9 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     TProfile* fVertexProfile;                    //! Profile of MC z-vtx - ESD z-vtx vs. MC z-vtx
     TH1F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed
 
+    TH2F* fEtaCorrelation;                       //! ESD eta vs MC eta
+    TProfile* fEtaProfile;                       //! Profile of MC eta - ESD eta vs. MC eta
+
     // histograms for systematic studies (must be enabled with option)
 
     TH1F* fSigmaVertexTracks;                    //! (accepted tracks) vs (n of sigma to vertex cut)
@@ -64,7 +69,7 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     TH1F* fMultTr; //! primary particles  in |eta| < 1 and pT > 0.2 in triggered events
     TH1F* fMultVtx; //! primary particles  in |eta| < 1 and pT > 0.2 in triggered events with vertex
 
-    TH1F* fDeltaPhi[3]; //! delta phi of primaries, secondaries, other (= unclear cases)
+    TH1F* fDeltaPhi[8]; //! delta phi of primaries, secondaries, other (= unclear cases)
 
     AlidNdEtaCorrection* fdNdEtaCorrectionProcessType[3]; //! correction for specific process type (ND, SD, DD)
                                                           // enable with option: process-types
index 90647ef..573ef5a 100644 (file)
@@ -45,6 +45,7 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fAnalysisMode(AliPWG0Helper::kTPC),
   fReadMC(kFALSE),
   fUseMCVertex(kFALSE),
+  fUseMCKine(kFALSE),
   fEsdTrackCuts(0),
   fdNdEtaAnalysisESD(0),
   fMult(0),
@@ -213,6 +214,57 @@ void AlidNdEtaTask::Exec(Option_t*)
   // post the data already here
   PostData(0, fOutput);
 
+  // needed for syst. studies
+  AliStack* stack = 0;
+  
+  if (fUseMCVertex || fUseMCKine) {
+    if (!fReadMC) {
+      Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
+      return;
+    }
+
+    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+    if (!eventHandler) {
+      Printf("ERROR: Could not retrieve MC event handler");
+      return;
+    }
+
+    AliMCEvent* mcEvent = eventHandler->MCEvent();
+    if (!mcEvent) {
+      Printf("ERROR: Could not retrieve MC event");
+      return;
+    }
+
+    AliHeader* header = mcEvent->Header();
+    if (!header)
+    {
+      AliDebug(AliLog::kError, "Header not available");
+      return;
+    }
+
+    if (fUseMCVertex)
+    {
+      Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
+      // get the MC vertex
+      AliGenEventHeader* genHeader = header->GenEventHeader();
+      TArrayF vtxMC(3);
+      genHeader->PrimaryVertex(vtxMC);
+
+      vtx[2] = vtxMC[2];
+    }
+
+    if (fUseMCKine)
+    {
+
+      stack = mcEvent->Stack();
+      if (!stack)
+      {
+        AliDebug(AliLog::kError, "Stack not available");
+        return;
+      }
+    }
+  }
+
   // create list of (label, eta, pt) tuples
   Int_t inputCount = 0;
   Int_t* labelArr = 0;
@@ -232,19 +284,32 @@ void AlidNdEtaTask::Exec(Option_t*)
     etaArr = new Float_t[mult->GetNumberOfTracklets()];
     ptArr = new Float_t[mult->GetNumberOfTracklets()];
 
+    if (fUseMCKine && stack)
+      Printf("Processing only primaries (MC information used). This is for systematical checks only.");
+
     // get multiplicity from ITS tracklets
     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
     {
       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 
-      // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
-      if (mult->GetDeltaPhi(i) < -1000)
-        continue;
+      if (fUseMCKine && stack)
+        if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+          continue;
+      
+      Float_t deltaPhi = mult->GetDeltaPhi(i);
+      // prevent values to be shifted by 2 Pi()
+      if (deltaPhi < -TMath::Pi())
+        deltaPhi += TMath::Pi() * 2;
+      if (deltaPhi > TMath::Pi())
+        deltaPhi -= TMath::Pi() * 2;
 
-      fDeltaPhi->Fill(mult->GetDeltaPhi(i));
+      if (TMath::Abs(deltaPhi) > 1)
+        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+
+      fDeltaPhi->Fill(deltaPhi);
 
       etaArr[inputCount] = mult->GetEta(i);
-      labelArr[inputCount] = mult->GetLabel(i);
+      labelArr[inputCount] = mult->GetLabel(i, 0);
       ptArr[inputCount] = 0; // no pt for tracklets
       ++inputCount;
     }
@@ -290,40 +355,6 @@ void AlidNdEtaTask::Exec(Option_t*)
   else
     return;
 
-  if (fUseMCVertex) {
-    Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
-    if (!fReadMC) {
-      Printf("ERROR: fUseMCVertex set without fReadMC set!");
-      return;
-    }
-
-    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-    if (!eventHandler) {
-      Printf("ERROR: Could not retrieve MC event handler");
-      return;
-    }
-
-    AliMCEvent* mcEvent = eventHandler->MCEvent();
-    if (!mcEvent) {
-      Printf("ERROR: Could not retrieve MC event");
-      return;
-    }
-
-    AliHeader* header = mcEvent->Header();
-    if (!header)
-    {
-      AliDebug(AliLog::kError, "Header not available");
-      return;
-    }
-
-    // get the MC vertex
-    AliGenEventHeader* genHeader = header->GenEventHeader();
-    TArrayF vtxMC(3);
-    genHeader->PrimaryVertex(vtxMC);
-
-    vtx[2] = vtxMC[2];
-  }
-
   // Processing of ESD information (always)
   if (eventTriggered)
   {
index 66a4073..3aaf075 100644 (file)
@@ -27,6 +27,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
     void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
     void SetReadMC(Bool_t flag = kTRUE) { fReadMC = flag; }
     void SetUseMCVertex(Bool_t flag = kTRUE) { fUseMCVertex = flag; }
+    void SetUseMCKine(Bool_t flag = kTRUE) { fUseMCKine = flag; }
 
  protected:
     AliESDEvent *fESD;    //! ESD object
@@ -37,6 +38,7 @@ class AlidNdEtaTask : public AliAnalysisTask {
 
     Bool_t  fReadMC;       // if true reads MC data (to build correlation maps)
     Bool_t  fUseMCVertex;  // the MC vtx is used instead of the ESD vertex (for syst. check)
+    Bool_t  fUseMCKine;    // Process only primaries by using the MC informatin (for syst. check)
 
     AliESDtrackCuts* fEsdTrackCuts;         // Object containing the parameters of the esd track cuts
 
index 61ff005..ce581c5 100644 (file)
@@ -45,7 +45,7 @@ void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, B
 
   task = new AlidNdEtaTask(option);
 
-  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC;
+  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD;
   task->SetAnalysisMode(analysisMode);
 
   if (analysisMode != AliPWG0Helper::kSPD)
@@ -66,6 +66,7 @@ void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, B
     task->SetReadMC();
 
   //task->SetUseMCVertex();
+  //task->SetUseMCKine();
 
   mgr->AddTask(task);
 
index 245d94f..174b22c 100644 (file)
@@ -45,9 +45,11 @@ void runCorrection(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug =
 
   task = new AlidNdEtaCorrectionTask(option);
 
-  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC;
+  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD;
   task->SetAnalysisMode(analysisMode);
 
+  //task->SetOnlyPrimaries();
+
   if (analysisMode != AliPWG0Helper::kSPD)
   {
     // selection of esd tracks
@@ -87,6 +89,8 @@ void runCorrection(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug =
     mgr->SetDebugLevel(2);
     AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kDebug+2);
   }
+  else
+    AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
 
   // Run analysis
   mgr->InitAnalysis();