This is rather large upgrade of the analysis. The sharing correction has been improve...
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Feb 2009 11:15:12 +0000 (11:15 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Feb 2009 11:15:12 +0000 (11:15 +0000)
14 files changed:
FMD/analysis/AliFMDAnalysisTaskBackgroundCorrection.cxx
FMD/analysis/AliFMDAnalysisTaskCollector.cxx
FMD/analysis/AliFMDAnalysisTaskCollector.h
FMD/analysis/AliFMDAnalysisTaskDensity.cxx
FMD/analysis/AliFMDAnalysisTaskDensity.h
FMD/analysis/AliFMDAnalysisTaskDndeta.cxx
FMD/analysis/AliFMDAnalysisTaskDndeta.h
FMD/analysis/AliFMDAnalysisTaskSE.cxx
FMD/analysis/AliFMDAnalysisTaskSharing.cxx
FMD/analysis/AliFMDAnalysisTaskSharing.h
FMD/analysis/RunAliEnFMDAnalysis.C
FMD/analysis/RunAliEnFMDAnalysisSE.C
FMD/analysis/RunLocalFMDAnalysis.C
FMD/analysis/RunLocalFMDAnalysisSE.C

index 06cf897..d62fc00 100644 (file)
@@ -74,7 +74,7 @@ void AliFMDAnalysisTaskBackgroundCorrection::CreateOutputObjects()
   
   
   TH2F* hMult = 0;
-  
+  TH2F* hHits = 0;
   Int_t nVtxbins = pars->GetNvtxBins();
   
   for(Int_t det =1; det<=3;det++)
@@ -100,18 +100,26 @@ void AliFMDAnalysisTaskBackgroundCorrection::CreateOutputObjects()
                              nSec, 0, 2*TMath::Pi());
            hMult->Sumw2();
            fOutputList->Add(hMult);
+           hHits  = new TH2F(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),
+                             hBg->GetNbinsX(),
+                             hBg->GetXaxis()->GetXmin(),
+                             hBg->GetXaxis()->GetXmax(),
+                             nSec, 0, 2*TMath::Pi());
+           
+           hHits->Sumw2();
+           fOutputList->Add(hHits);
            vtxArray->AddAtAndExpand(hMult,i);
            
          }
        } 
     }
-  
   if(fStandalone) {
     fOutputVertexString = new TObjString();
   }
   fOutputList->Add(fOutputVertexString);
   
   
+  
 }
 //_____________________________________________________________________
 void AliFMDAnalysisTaskBackgroundCorrection::ConnectInputData(Option_t */*option*/)
@@ -161,11 +169,14 @@ void AliFMDAnalysisTaskBackgroundCorrection::Exec(Option_t */*option*/)
       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
       TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
       TH2F* hMultInput = (TH2F*)vtxInputArray->At(vtxbin);
+      TH2F* hHits      = (TH2F*)fOutputList->FindObject(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
+      
+      hHits->Add(hMultInput);
       TH2F* hBg        = pars->GetBackgroundCorrection(det, ringChar, vtxbin);
       
       TH2F* hTmp       = (TH2F*)hMultInput->Clone("hMult_from_event");
-            
-      hTmp->Divide(hTmp,hBg,1,1,"B");
+      
+      hTmp->Divide(hTmp,hBg,1,1);//,"B");
       
       hMultTotal->Add(hTmp);
       delete hTmp;
index e91d7b4..d76dbc6 100644 (file)
@@ -107,10 +107,18 @@ void AliFMDAnalysisTaskCollector::Exec(Option_t */*option*/)
   if (old) {
     fESD->CopyFromOldESD();
   }
+  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
   
   Double_t vertex[3];
-  fESD->GetVertex()->GetXYZ(vertex);
+  
+  GetVertex(vertex);
+  if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
+    return;
   fZvtxDist->Fill(vertex[2]);
+  
+  if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
+    return;
+  
   AliESDFMD* fmd = fESD->GetFMDData();
   if (!fmd) return;
   
@@ -152,7 +160,34 @@ void AliFMDAnalysisTaskCollector::Terminate(Option_t */*option*/)
   }
   */
 }
-
+//_____________________________________________________________________
+void AliFMDAnalysisTaskCollector::GetVertex(Double_t* vertexXYZ) 
+{
+  const AliESDVertex* vertex = 0;
+  vertex = fESD->GetPrimaryVertex();
+  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
+    vertex = fESD->GetPrimaryVertexSPD();
+  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
+    vertex = fESD->GetPrimaryVertexTPC();
+  
+  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))    
+    vertex = fESD->GetVertex();
+  
+  if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
+    vertex->GetXYZ(vertexXYZ);
+    //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
+    return;
+  }
+  else if (fESD->GetESDTZERO()) { 
+    vertexXYZ[0] = 0;
+    vertexXYZ[1] = 0;
+    vertexXYZ[2] = fESD->GetT0zVertex();
+    
+    return;
+  }
+  
+  return;
+}
 //____________________________________________________________________
 //
 // EOF
index a7de28f..eb0117a 100644 (file)
@@ -40,6 +40,7 @@ class AliFMDAnalysisTaskCollector : public AliAnalysisTask
     virtual void SetDebugLevel(Int_t level) {fDebug = level;}
     
  private:
+    void          GetVertex(Double_t* vertexXYZ); 
     Int_t         fDebug;        //  Debug flag
     TChain*       fChain;        //! chained files
     AliESDEvent*  fESD;          //! ESD
index 9162295..67db948 100644 (file)
@@ -32,7 +32,8 @@ AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
   fESD(0x0),
   fVertexString(),
   fVertex(0),
-  fStandalone(kTRUE)
+  fStandalone(kTRUE),
+  fStatus(kTRUE)
 {
   // Default constructor
   DefineInput (0, AliESDFMD::Class());
@@ -48,7 +49,8 @@ AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE
     fESD(0x0),
     fVertexString(),
     fVertex(0),
-    fStandalone(kTRUE)
+    fStandalone(kTRUE),
+    fStatus(kTRUE)
 {
   fStandalone = SE;
   if(fStandalone) {
@@ -68,6 +70,9 @@ void AliFMDAnalysisTaskDensity::CreateOutputObjects()
     fOutputList = new TList();
   fOutputList->SetName("density_list");
   
+  fOutputList->Add(&fArray);
+  fOutputList->Add(&fVertexString);
+  
   TH2F* hMult = 0;
   
   Int_t nVtxbins = pars->GetNvtxBins();
@@ -94,14 +99,13 @@ void AliFMDAnalysisTaskDensity::CreateOutputObjects()
                              hBg->GetXaxis()->GetXmin(),
                              hBg->GetXaxis()->GetXmax(),
                              nSec, 0, 2*TMath::Pi());
-           vtxArray->AddAtAndExpand(hMult,i);
            
+           vtxArray->AddAtAndExpand(hMult,i);
          }
        } 
     }
   
-  fOutputList->Add(&fArray);
-  fOutputList->Add(&fVertexString);
+  
   
   
 }
@@ -124,8 +128,13 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
   Double_t vertex[3];
   fVertex->GetXYZ(vertex);
   // Z Vtx cut
-  if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) 
+  if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
+    fStatus = kFALSE;
     return;
+  }
+  else
+    fStatus = kTRUE;
+  
   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
   
@@ -153,13 +162,20 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
       
       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
+      
       Char_t   ring = (ir == 0 ? 'I' : 'O');
       UShort_t nsec = (ir == 0 ? 20  : 40);
       UShort_t nstr = (ir == 0 ? 512 : 256);
+      
       for(UShort_t sec =0; sec < nsec;  sec++)  {
        for(UShort_t strip = 0; strip < nstr; strip++) {
          Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
-         if(mult < 1 || mult == AliESDFMD::kInvalidMult) continue;
+         Float_t mult_cut = 0.1;
+         if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
+         //Particle number cut goes here...
+         Float_t nParticles = 0;
+         if(mult > mult_cut)
+           nParticles = 1;
          Float_t eta = fESD->Eta(det,ring,sec,strip);
          Double_t x,y,z;
          geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
@@ -168,11 +184,14 @@ void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
            phi = phi+2*TMath::Pi();
          Float_t correction = GetAcceptanceCorrection(ring,strip);
          if(correction) mult = mult / correction;
-         hMult->Fill(eta,phi,mult);
+         hMult->Fill(eta,phi,nParticles);
+         
          
        }
       }
+      
     }
+    
        
   
   }
index 3fefd7c..33dcf7f 100644 (file)
@@ -31,7 +31,8 @@ class AliFMDAnalysisTaskDensity : public AliAnalysisTask
       fESD(o.fESD),
       fVertexString(o.fVertexString),
       fVertex(o.fVertex),
-      fStandalone(o.fStandalone) {}
+      fStandalone(o.fStandalone),
+      fStatus(o.fStatus) {}
     AliFMDAnalysisTaskDensity& operator=(const AliFMDAnalysisTaskDensity&) { return *this; }
     // Implementation of interface methods
     virtual void ConnectInputData(Option_t *option);
@@ -46,7 +47,7 @@ class AliFMDAnalysisTaskDensity : public AliAnalysisTask
     void SetOutputList(TList* outlist) {fOutputList = outlist;}
     void SetInputESDFMD(AliESDFMD* esdfmd) {fESD = esdfmd;}
     void SetInputVertex(AliESDVertex* vertex) {fVertex = vertex;}
-   
+    Bool_t GetEventStatus() { return fStatus; }
  private:
     
     Int_t         fDebug;        //  Debug flag
@@ -56,6 +57,7 @@ class AliFMDAnalysisTaskDensity : public AliAnalysisTask
     TObjString    fVertexString;
     AliESDVertex* fVertex;
     Bool_t        fStandalone;
+    Bool_t        fStatus;
     
     ClassDef(AliFMDAnalysisTaskDensity, 0); // Analysis task for FMD analysis
 };
index d2559f0..acf4104 100644 (file)
@@ -6,6 +6,7 @@
 #include <TFile.h>
 #include <TList.h>
 #include <iostream>
+#include "TH1F.h"
 #include "TH2F.h"
 #include "AliFMDAnalysisTaskDndeta.h"
 #include "AliAnalysisManager.h"
@@ -20,6 +21,7 @@
 #include "TMath.h"
 #include "AliFMDAnaParameters.h"
 #include "AliFMDGeometry.h"
+#include "AliGenEventHeader.h"
 
 ClassImp(AliFMDAnalysisTaskDndeta)
 
@@ -32,7 +34,9 @@ AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
   fInputArray(0),
   fVertexString(0x0),
   fNevents(),
-  fStandalone(kTRUE)
+  fNMCevents(),
+  fStandalone(kTRUE),
+  fMCevent(0)
 {
   // Default constructor
   DefineInput (0, TList::Class());
@@ -48,7 +52,9 @@ AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
     fInputArray(0),
     fVertexString(0x0),
     fNevents(),
-    fStandalone(kTRUE)
+    fNMCevents(),
+    fStandalone(kTRUE),
+    fMCevent(0)
 {
   fStandalone = SE;
   if(fStandalone) {
@@ -72,9 +78,19 @@ void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
   
   
   TH2F* hMult = 0;
+  TH1F* hPrimVertexBin = 0;
   
+  
+  TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', 0);
+  TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
+                           hBg->GetNbinsX(),
+                           hBg->GetXaxis()->GetXmin(),
+                           hBg->GetXaxis()->GetXmax());
+  hPrimary->Sumw2();
+  fOutputList->Add(hPrimary);
   Int_t nVtxbins = pars->GetNvtxBins();
   
+  
   for(Int_t det =1; det<=3;det++)
     {
       TObjArray* detArray = new TObjArray();
@@ -104,10 +120,24 @@ void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
        } 
     }
   
+  for(Int_t i = 0; i< nVtxbins; i++) {
+   
+    hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
+                             Form("primmult_vtxbin%d",i),
+                             hBg->GetNbinsX(),
+                             hBg->GetXaxis()->GetXmin(),
+                             hBg->GetXaxis()->GetXmax());
+    hPrimVertexBin->Sumw2();
+    fOutputList->Add(hPrimVertexBin);
+    
+  }
+  
   fNevents.SetBins(nVtxbins,0,nVtxbins);
   fNevents.SetName("nEvents");
+  fNMCevents.SetBins(nVtxbins,0,nVtxbins);
+  fNMCevents.SetName("nMCEvents");
   fOutputList->Add(&fNevents);
-   
+  fOutputList->Add(&fNMCevents);
   
 }
 //_____________________________________________________________________
@@ -142,6 +172,10 @@ void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
       
     }
   }
+  
+  if(fMCevent)
+    ProcessPrimary();
+  
   if(fStandalone) {
     PostData(0, fOutputList); 
   }
@@ -171,6 +205,44 @@ void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
   
 }
 //_____________________________________________________________________
+void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
+
+  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  
+  AliMCParticle* particle = 0;
+  AliStack* stack = fMCevent->Stack();
+  
+  TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
+  
+  Bool_t firstTrack = kTRUE;
+  Int_t nTracks = fMCevent->GetNumberOfTracks();
+  for(Int_t i = 0 ;i<nTracks;i++) {
+    particle = fMCevent->GetTrack(i);
+    if(!particle)
+      continue;
+    if(TMath::Abs(particle->Zv()) > pars->GetVtxCutZ())
+      continue;
+    if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
+      hPrimary->Fill(particle->Eta());
+      Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
+      Double_t vertexBinDouble = (particle->Zv() + pars->GetVtxCutZ()) / delta;
+      Int_t    vertexBin       = (Int_t)vertexBinDouble;
+      TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
+      hPrimVtxBin->Fill(particle->Eta());
+      if(firstTrack) {
+       fNMCevents.Fill(vertexBin);
+       firstTrack = kFALSE;
+      }
+    }
+      
+  }
+  
+  
+  
+  
+
+}
+//_____________________________________________________________________
 //
 //
 // EOF
index 7c604a2..c12b76a 100644 (file)
@@ -10,6 +10,7 @@
 #include "TObjString.h"
 #include "TArrayI.h"
 #include "TH1I.h"
+#include "AliMCEvent.h"
 
 class AliFMDAnalysisTaskDndeta : public AliAnalysisTask
 {
@@ -25,7 +26,9 @@ class AliFMDAnalysisTaskDndeta : public AliAnalysisTask
       fInputArray(o.fInputArray),
       fVertexString(o.fVertexString),
       fNevents(o.fNevents),
-      fStandalone(o.fStandalone) {}
+      fNMCevents(o.fNMCevents),
+      fStandalone(o.fStandalone),
+      fMCevent(o.fMCevent) {}
     AliFMDAnalysisTaskDndeta& operator=(const AliFMDAnalysisTaskDndeta&) { return *this; }
     // Implementation of interface methods
     virtual void ConnectInputData(Option_t *option = "");
@@ -38,7 +41,10 @@ class AliFMDAnalysisTaskDndeta : public AliAnalysisTask
     void SetInputList(TList* inputList) {fInputList = inputList;}
     void SetInputVertex(TObjString* vtxString) {fVertexString = vtxString;}
     void SetOutputList(TList* outputList) {fOutputList = outputList;}
+    void SetMCEvent(AliMCEvent* mcevent) {fMCevent = mcevent;}
+    void ProcessPrimary();
     TList* GetOutputList() {return fOutputList;}
+    
  private:
     Int_t         fDebug;        //  Debug flag
     TList*        fOutputList;
@@ -47,7 +53,9 @@ class AliFMDAnalysisTaskDndeta : public AliAnalysisTask
     TObjArray*    fInputArray;
     TObjString*   fVertexString;
     TH1I          fNevents;
+    TH1I          fNMCevents;
     Bool_t        fStandalone;
+    AliMCEvent*   fMCevent;
     ClassDef(AliFMDAnalysisTaskDndeta, 0); // Analysis task for FMD analysis
 };
  
index 09df659..3a0a2fe 100644 (file)
@@ -81,9 +81,18 @@ void AliFMDAnalysisTaskSE::UserExec(Option_t */*option*/)
   fSharing.SetInputESD(fESD);
   
   fSharing.Exec("");
-  fDensity.Exec("");
-  fBackground.Exec("");  
-  fDndeta.Exec("");
+  if(fSharing.GetEventStatus()) {
+    fDensity.Exec("");
+    if(fDensity.GetEventStatus()) {
+      fBackground.Exec("");  
+      AliMCEvent* mcevent = MCEvent();
+      fDndeta.SetMCEvent(mcevent);
+      fDndeta.Exec("");
+    }
+  }
+  else
+    return;
+  
   //fListOfHistos = fBackground.GetOutputList();
   
   PostData(1, fListOfHistos);
index acf3e40..d96d79c 100644 (file)
@@ -6,14 +6,17 @@
 #include <TFile.h>
 #include <TList.h>
 #include <iostream>
-
+#include <TMath.h>
+#include "AliFMDDebug.h"
 #include "AliFMDAnalysisTaskSharing.h"
 #include "AliAnalysisManager.h"
 #include "AliESDFMD.h"
+#include "AliFMDGeometry.h"
 #include "AliMCEventHandler.h"
 #include "AliStack.h"
 #include "AliESDVertex.h"
 #include "AliFMDAnaParameters.h"
+#include "AliFMDParameters.h"
 
 ClassImp(AliFMDAnalysisTaskSharing)
 
@@ -26,7 +29,8 @@ AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
   fSharedPrev(kFALSE),
   fDiagList(),
   fStandalone(kTRUE),
-  fEsdVertex(0)
+  fEsdVertex(0),
+  fStatus(kTRUE)
 {
   // Default constructor
   DefineInput (0, AliESDEvent::Class());
@@ -41,11 +45,14 @@ AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE
     fDebug(0),
     fESD(0x0),
     foutputESDFMD(),
+    fEnergy(0),
+    fNstrips(0),
     fSharedThis(kFALSE),
     fSharedPrev(kFALSE),
     fDiagList(),
     fStandalone(kTRUE),
-    fEsdVertex(0)
+    fEsdVertex(0),
+    fStatus(kTRUE)
 {
   fStandalone = SE;
   if(fStandalone) {
@@ -73,12 +80,18 @@ void AliFMDAnalysisTaskSharing::CreateOutputObjects()
       Char_t ringChar = (iring == 0 ? 'I' : 'O');
       TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
                                     Form("Edist_before_sharing_FMD%d%c", det, ringChar),
-                                    200,0,5);
+                                    1000,0,25);
       TH1F* hEdist_after  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
                                     Form("Edist_after_sharing_FMD%d%c", det, ringChar),
-                                    200,0,5);
+                                    1000,0,25);
+      
+      
+      TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
+                                    Form("N_strips_hit_FMD%d%c",det,ringChar),
+                                    25,0,25);
       fDiagList.Add(hEdist);
       fDiagList.Add(hEdist_after);
+      fDiagList.Add(hNstripsHit);
 
     }
   }
@@ -92,7 +105,7 @@ void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
 //_____________________________________________________________________
 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
 {
-  
+
   AliESD* old = fESD->GetAliESDOld();
   if (old) {
     fESD->CopyFromOldESD();
@@ -100,10 +113,21 @@ void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
   
   foutputESDFMD->Clear();
   
+  Double_t vertex[3];
+  GetVertex(vertex);
+  fEsdVertex->SetXYZ(vertex);
+  
+  if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
+    fStatus = kFALSE;
+    return;
+  }
+  else
+    fStatus = kTRUE;
+  
   AliESDFMD* fmd = fESD->GetFMDData();
   
   if (!fmd) return;
-  
+  Int_t nHits = 0;
   for(UShort_t det=1;det<=3;det++) {
     Int_t nRings = (det==1 ? 1 : 2);
     for (UShort_t ir = 0; ir < nRings; ir++) {
@@ -116,34 +140,48 @@ void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
       for(UShort_t sec =0; sec < nsec;  sec++) {
        fSharedThis      = kFALSE;
        fSharedPrev      = kFALSE;
+       fEnergy = 0;
+       fNstrips = 0;
        for(UShort_t strip = 0; strip < nstr; strip++) {
          foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
          Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
-
+         
+         
          if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
          
+         Double_t eta  = fmd->Eta(det,ring,sec,strip);//EtaFromStrip(det,ring,sec,strip,vertex[2]);
+         //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
          hEdist->Fill(mult);
+         if(fmd->IsAngleCorrected())
+           mult = mult/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip)));
          Float_t Eprev = 0;
          Float_t Enext = 0;
          if(strip != 0)
-           if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult)
+           if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
              Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
+             if(fmd->IsAngleCorrected())
+               Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
+           }
          if(strip != nstr - 1)
-           if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult)
-           Enext = fmd->Multiplicity(det,ring,sec,strip+1);
+           if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
+             Enext = fmd->Multiplicity(det,ring,sec,strip+1);
+             if(fmd->IsAngleCorrected())
+               Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
+           }
          
-         Float_t nParticles = GetMultiplicityOfStrip(mult,Eprev,Enext,det,ring);
-         foutputESDFMD->SetMultiplicity(det,ring,sec,strip,nParticles);
-         foutputESDFMD->SetEta(det,ring,sec,strip,fmd->Eta(det,ring,sec,strip));
+         Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);
+         if(merged_energy > 0 )
+           nHits++;
+         foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
+         foutputESDFMD->SetEta(det,ring,sec,strip,eta);
          
        }
       }
     }
   }
   
-  Double_t vertex[3];
-  GetVertex(vertex);
-  fEsdVertex->SetXYZ(vertex);
+  //std::cout<<fESD->GetEventNumberInFile()<<"    "<<nHits<<"   "<<vertex[2]<<std::endl;
+  
   if(fStandalone) {
     PostData(0, foutputESDFMD); 
     PostData(1, fEsdVertex); 
@@ -153,36 +191,94 @@ void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
 }
 //_____________________________________________________________________
 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
+                                                         Float_t eta,
                                                          Float_t Eprev,
                                                          Float_t Enext,
-                                                         Int_t   det,
-                                                         Char_t  ring) {
+                                                         UShort_t   det,
+                                                         Char_t  ring,
+                                                         UShort_t sec,
+                                                         UShort_t strip) {
   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+  AliFMDParameters* recopars = AliFMDParameters::Instance();
+  Float_t merged_energy = 0;
   Float_t nParticles = 0;
   Float_t cutLow  = 0.2;
-  Float_t cutHigh = pars->GetMPV(det,ring) - 2*pars->GetSigma(det,ring);
-  Float_t Etotal = mult;
-  /* 
-  if(mult > 3*pars->GetMPV(det,ring) && 
-     (Enext > 3*pars->GetMPV(det,ring) || (Enext > 3*pars->GetMPV(det,ring))))
-    return 0;
-  
-  if(mult > 5*pars->GetMPV(det,ring))
-    return 0;
-  */
+  //  Float_t gain = recopars->GetPulseGain(det,ring,sec,strip);
+  //Float_t pw   = recopars->GetPedestalWidth(det,ring,sec,strip);
+  //Float_t cutLow  =  (4*pw)/(gain*recopars->GetDACPerMIP());
+  //std::cout<<cutLow<<"    "<<gain<<"    "<<recopars->GetEdepMip()<<std::endl;
+  Float_t cutHigh = pars->GetMPV(det,ring) -pars->GetSigma(det,ring);
+  Float_t cutPart = pars->GetMPV(det,ring) - 5*pars->GetSigma(det,ring);
+  // Float_t cutPart = cutLow;//0.33*pars->GetMPV(det,ring);
+  
+  // if(ring == 'O')
+  //  cutPart = pars->GetMPV(det,ring) - 2*pars->GetSigma(det,ring);
+  Float_t Etotal  = mult;
+  //std::cout<<mult<<"  "<<strip<<std::endl;
+  
+  if(mult > 0 ) {
+    fEnergy = fEnergy + mult;
+    fNstrips++;
+  }
+  if((Enext <0.01 && fEnergy >0) || fNstrips >4 ) {
+    
+          
+    if((fEnergy*TMath::Cos(Eta2Theta(eta))) > cutPart || fNstrips > 1) {
+      nParticles = 1;
+      merged_energy = fEnergy*TMath::Cos(Eta2Theta(eta));
+      TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
+      hEdist->Fill(fEnergy);
+      TH1F* hNstrips = (TH1F*)fDiagList.FindObject(Form("N_strips_hit_FMD%d%c",det,ring));
+      hNstrips->Fill(fNstrips);
+      //  std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det )<<std::endl;
+      
+    }
+    // else
+    //std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
+    
+    fEnergy  = 0;
+    fNstrips = 0;
+    return merged_energy;
+  }
+  
+  return 0;
+  
+  
+  /*
   if(fSharedThis) {
     fSharedThis      = kFALSE;
     fSharedPrev      = kTRUE;
     return 0.;
   }
   
-  if(Etotal < 0.33*pars->GetMPV(det,ring)) {
-    fSharedThis      = kFALSE;
-    fSharedPrev      = kFALSE;
-    return 0.; 
-  }
+  if(mult < cutLow)
+     return 0;
+     
+  
+  //if(Etotal < 0.33*pars->GetMPV(det,ring)) {
+  //  fSharedThis      = kFALSE;
+  //  fSharedPrev      = kFALSE;
+  //  return 0.; 
+ // }
+  
+  //Experimental cut 
+  if(mult<Enext && Enext>cutHigh)
+    {
+      fSharedThis      = kFALSE;
+      fSharedPrev      = kFALSE;
+      return 0;
+    }
+   //Experimental cut 2
+  if(mult>Enext && mult > cutHigh)
+    {
+      fSharedThis      = kTRUE;
+      fSharedPrev      = kFALSE;
+      nParticles = 1;
+    }
+  
   
-  if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
+  
+     if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
     Etotal += Eprev;
   }
   
@@ -190,42 +286,89 @@ Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
     Etotal += Enext;
     fSharedThis      = kTRUE;
   }
-  TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
-  hEdist->Fill(Etotal);
-  if(Etotal > cutHigh ) {
-    
+      TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
+     hEdist->Fill(Etotal);
+  
+     Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
+     if(Etotal > cutPart) {
     nParticles = 1;
     fSharedPrev      = kTRUE;
+    // std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det )<<std::endl;
   }
   else {
+    if(Etotal > 0)
+      //std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
     fSharedThis      = kFALSE;
     fSharedPrev      = kFALSE;
   }
   
-  return nParticles;
+     return nParticles; 
+     */
 }
 //_____________________________________________________________________
 void AliFMDAnalysisTaskSharing::GetVertex(Double_t* vertexXYZ) 
 {
   const AliESDVertex* vertex = 0;
   vertex = fESD->GetPrimaryVertex();
-  if (!vertex)        vertex = fESD->GetPrimaryVertexSPD();
-  if (!vertex)        vertex = fESD->GetPrimaryVertexTPC();
-  if (!vertex)        vertex = fESD->GetVertex();
-  
-  if (vertex) {
+  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
+    vertex = fESD->GetPrimaryVertexSPD();
+  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
+    vertex = fESD->GetPrimaryVertexTPC();
+  if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))    
+    vertex = fESD->GetVertex();
+  if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
     vertex->GetXYZ(vertexXYZ);
+    //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
     return;
   }
   else if (fESD->GetESDTZERO()) { 
     vertexXYZ[0] = 0;
     vertexXYZ[1] = 0;
     vertexXYZ[2] = fESD->GetT0zVertex();
-
+    
     return;
   }
   
   return;
+  
+}
+//_____________________________________________________________________
+Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
+
+  Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
+  
+  if(eta < 0)
+    theta = theta-TMath::Pi();
+  
+  // std::cout<<"From eta2Theta: "<<theta<<std::endl;
+  return theta;
+  
+
+
+}
+//_____________________________________________________________________
+Double_t AliFMDAnalysisTaskSharing::EtaFromStrip(UShort_t det, 
+                                               Char_t ring, 
+                                               UShort_t sector, 
+                                               UShort_t strip, 
+                                               Double_t zvtx)
+{
+   
+  AliFMDGeometry* geo = AliFMDGeometry::Instance();
+  
+  Double_t x,y,z;
+  geo->Detector2XYZ(det,ring,sector,strip,x,y,z);
+  
+  Double_t r = TMath::Sqrt(x*x+y*y);
+  
+  Double_t z_real      = z-zvtx;
+  Double_t theta       = TMath::ATan2(r,z_real);
+  // std::cout<<"From EtaFromStrip "<<theta<<std::endl;
+  Double_t eta         =  -1*TMath::Log(TMath::Tan(0.5*theta));
+  
+  // std::cout<<det<<"   "<<ring<<"   "<<sector<<"   "<<strip<<"   "<<r<<"    "<<z_real<<"   "<<theta<<"    "<<eta<<std::endl;
+  
+  return eta;
 }
 //_____________________________________________________________________
 //
index 7258056..3a6e6bb 100644 (file)
@@ -26,11 +26,14 @@ class AliFMDAnalysisTaskSharing : public AliAnalysisTask
       fESD(o.fESD),
       // fOutputESD(),
       foutputESDFMD(o.foutputESDFMD),
+      fEnergy(o.fEnergy),
+      fNstrips(o.fNstrips),
       fSharedThis(o.fSharedThis),
       fSharedPrev(o.fSharedPrev),
       fDiagList(),
       fStandalone(o.fStandalone),
-      fEsdVertex(o.fEsdVertex) {}
+      fEsdVertex(o.fEsdVertex),
+      fStatus(o.fStatus)   {}
     AliFMDAnalysisTaskSharing& operator=(const AliFMDAnalysisTaskSharing&) { return *this; }
     
     // Implementation of interface methods
@@ -38,24 +41,31 @@ class AliFMDAnalysisTaskSharing : public AliAnalysisTask
     virtual void CreateOutputObjects();
     virtual void Init() {}
     virtual void LocalInit() {Init();}
-    virtual void Exec(Option_t *option);
+    virtual void Exec(Option_t */*option*/);
     virtual void Terminate(Option_t* /* option*/) {}
     virtual void SetDebugLevel(Int_t level) {fDebug = level;}
-    Float_t GetMultiplicityOfStrip(Float_t mult, Float_t Eprev, Float_t Enext, Int_t   det, Char_t  ring);
+    Float_t GetMultiplicityOfStrip(Float_t mult, Float_t eta, Float_t Eprev, Float_t Enext, UShort_t   det, Char_t  ring, UShort_t sec, UShort_t strip);
     void GetVertex(Double_t* vertexXYZ) ;
     void SetFMDData(AliESDFMD* fmd) {foutputESDFMD = fmd;}
     void SetVertex(AliESDVertex* vertex) {fEsdVertex = vertex;}
     void SetInputESD(AliESDEvent* esd) {fESD = esd;}
+    Bool_t GetEventStatus() {return fStatus;}
  private:
+    Float_t Eta2Theta(Float_t eta);
+    Double_t EtaFromStrip(UShort_t det, Char_t ring, UShort_t sector, UShort_t strip, Double_t zvtx);
     Int_t         fDebug;        //  Debug flag
     AliESDEvent*  fESD;          //! ESD
     // AliESDEvent   fOutputESD;
     AliESDFMD*    foutputESDFMD;
+    Float_t       fEnergy;
+    Int_t         fNstrips;
     Bool_t        fSharedThis;
     Bool_t        fSharedPrev;
     TList         fDiagList;
     Bool_t        fStandalone;
     AliESDVertex* fEsdVertex;
+    Bool_t        fStatus;
+
     ClassDef(AliFMDAnalysisTaskSharing, 0); // Analysis task for FMD analysis
 };
  
index b4bd963..b718406 100644 (file)
@@ -77,7 +77,7 @@ void RunAliEnFMDAnalysis(const Char_t* collectionfile = "collection.xml",
   TAlienCollection* coll =  TAlienCollection::Open(collectionfile);  
   coll->Reset();
   Int_t nFiles = 0;
-  while(coll->Next() && nFiles<10) {
+  while(coll->Next() && nFiles<2) {
     cout<<coll->GetTURL("")<<endl;
     TString test(coll->GetTURL(""));
     chain->Add(coll->GetTURL(""));
index f4c0992..e43e51b 100644 (file)
@@ -15,7 +15,10 @@ void RunAliEnFMDAnalysisSE(const Char_t* collectionName="collection.xml", const
   // ESD input handler
   AliESDInputHandler *esdHandler = new AliESDInputHandler();
   mgr->SetInputEventHandler(esdHandler);
-
+  
+  AliMCEventHandler *mcHandler = new AliMCEventHandler();
+  mgr->SetMCtruthEventHandler(mcHandler);
+  
   AliAODHandler* aodHandler   = new AliAODHandler();
   mgr->SetOutputEventHandler(aodHandler);
   aodHandler->SetOutputFileName("AliAODs.root");
@@ -61,7 +64,7 @@ void RunAliEnFMDAnalysisSE(const Char_t* collectionName="collection.xml", const
   timer.Start();
   if (mgr->InitAnalysis()) {
     mgr->PrintStatus();
-    mgr->StartAnalysis("local",chain, 1000);
+    mgr->StartAnalysis("local",chain, 100);
   }   
   timer.Stop();
   timer.Print();
index 45133cd..0af36d6 100644 (file)
@@ -9,7 +9,8 @@ void RunLocalFMDAnalysis(const Char_t* filename= "AliESDs.root",
   AliCDBManager* cdb = AliCDBManager::Instance();
   cdb->SetDefaultStorage(cdbPath);
   cdb->SetRun(0);
-  
+  AliFMDParameters* recopars = AliFMDParameters::Instance();
+  recopars->Init();
   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
   pars->Init();
   if (AliGeomManager::GetGeometry() == NULL)
@@ -42,9 +43,9 @@ void RunLocalFMDAnalysis(const Char_t* filename= "AliESDs.root",
   AliAnalysisDataContainer* cdiag1 = mgr->CreateContainer("diagSharing1",AliESDEvent::Class(),AliAnalysisManager::kExchangeContainer);
   AliAnalysisDataContainer* cdiag2 = mgr->CreateContainer("diagSharing2",TList::Class(),AliAnalysisManager::kOutputContainer,"edists.root");
   AliAnalysisDataContainer* cexchange1 = mgr->CreateContainer("exchangeESDFMD1",AliESDFMD::Class(),AliAnalysisManager::kExchangeContainer);
-  AliAnalysisDataContainer* cexchange2 = mgr->CreateContainer("listOfhists",TList::Class(),AliAnalysisManager::kExchangeContainer);
+  AliAnalysisDataContainer* cexchange2 = mgr->CreateContainer("list_of_hits",TList::Class(),AliAnalysisManager::kExchangeContainer);
   AliAnalysisDataContainer* cvertex = mgr->CreateContainer("vertex",TObjString::Class(),AliAnalysisManager::kExchangeContainer);
-  AliAnalysisDataContainer* cexchange3 = mgr->CreateContainer("BackgroundCorrectedperevent",TList::Class(),AliAnalysisManager::kOutputContainer,"testOut.root");
+  AliAnalysisDataContainer* cexchange3 = mgr->CreateContainer("list_of_hits_and_mult",TList::Class(),AliAnalysisManager::kOutputContainer,"testOut.root");
   AliAnalysisDataContainer* coutput = mgr->CreateContainer("BackgroundCorrected",TList::Class(),AliAnalysisManager::kOutputContainer,outFile);
   
   mgr->ConnectInput(FMDana0, 0 , cin_esd);   
@@ -81,7 +82,7 @@ void RunLocalFMDAnalysis(const Char_t* filename= "AliESDs.root",
   TStopwatch timer;
   timer.Start();
   cout<<"Executing analysis"<<endl;
-  mgr->StartAnalysis("local",chain);
+  mgr->StartAnalysis("local",chain,1000);
   timer.Stop();
   timer.Print();
 }
index 1120aa4..a9f985c 100644 (file)
@@ -17,7 +17,10 @@ void RunLocalFMDAnalysisSE(const Char_t* filename= "AliESDs.root", const Char_t*
   // ESD input handler
   AliESDInputHandler *esdHandler = new AliESDInputHandler();
   mgr->SetInputEventHandler(esdHandler);
-
+  
+  AliMCEventHandler *mcHandler = new AliMCEventHandler();
+  mgr->SetMCtruthEventHandler(mcHandler);
+  
   AliAODHandler* aodHandler   = new AliAODHandler();
   mgr->SetOutputEventHandler(aodHandler);
   aodHandler->SetOutputFileName("AliAODs.root");