class for multiplicity analysis and macro to run it on Proof
authoralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jun 2009 13:49:27 +0000 (13:49 +0000)
committeralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jun 2009 13:49:27 +0000 (13:49 +0000)
T0/AliT0MultiplicityTask.cxx [new file with mode: 0644]
T0/AliT0MultiplicityTask.h [new file with mode: 0644]
T0/runProofT0analysis.C [new file with mode: 0644]

diff --git a/T0/AliT0MultiplicityTask.cxx b/T0/AliT0MultiplicityTask.cxx
new file mode 100644 (file)
index 0000000..cd6ee40
--- /dev/null
@@ -0,0 +1,304 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "TClonesArray.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TLine.h"
+#include "TText.h"
+#include "TFile.h"
+#include "TBenchmark.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliLog.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliESDtrack.h"
+#include "AliStack.h"
+#include "AliTrackReference.h"
+
+#include "AliT0MultiplicityTask.h"
+
+ClassImp(AliT0MultiplicityTask)
+
+
+AliT0MultiplicityTask::AliT0MultiplicityTask()
+  : AliAnalysisTaskSE("AliT0MultiplicityTask"),
+    fListOfHistos(0),
+    fOrA(0),
+    fOrC(0),
+    fMean(0),
+    fVertex(0),
+    fTime(0),
+    fAmp(0),
+    fTotalMult(0),
+    fMultRecTotal(0),
+    fMultRecRealA(0),
+    fMultRecRealC(0),
+    fPrim(0),
+    fRef(0),
+    fRec(0)
+{
+  // Default constructor
+  AliInfo("Default constructor AliT0MultiplicityTask");
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #1 TList
+  DefineOutput(1, TList::Class());
+}
+
+
+AliT0MultiplicityTask::AliT0MultiplicityTask(const char* name)
+  : AliAnalysisTaskSE(name),
+    fListOfHistos(0),
+    fOrA(0),
+    fOrC(0),
+    fMean(0),
+    fVertex(0),
+    fTime(0),
+    fAmp(0),   
+    fTotalMult(0),
+    fMultRecTotal(0),
+    fMultRecRealA(0),
+    fMultRecRealC(0),
+    fPrim(0),
+    fRef(0),
+    fRec(0)
+{
+  // Constructor
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #1 TList
+  DefineOutput(1, TList::Class());
+}
+
+
+void AliT0MultiplicityTask::UserCreateOutputObjects() 
+{
+  // Create histograms
+  // Called once
+  // Create output container
+  fListOfHistos = new TList();
+  
+  fOrA = new TH1F("hOrA"," T0 A", 100,2450,2700);
+  fOrC = new TH1F("hOrC"," T0 C", 100,2450,2700);
+  fMean= new TH1F("hMean"," T0 ",100, 12000, 13000);
+  fVertex = new TH1F("hVertex","Z position of vertex",   100,-30, 30);
+   
+  fTime = new TH1F("Time", " Amp vs Time",100, 12000, 13000);
+  fAmp = new TH1F("fAmp", " Amplitude", 100, 0, 200);
+
+  //  fHighMult  = new TH1F("fHighMult", " events with high amp ",100, 0, 10);
+  
+  fTotalMult = new TH1F("fTotalMult","total multiplicity",500,0,5000);
+  
+  fMultRecRealA = new TH2F("fMultRecRealA"," ",100,0.,200,100,0,200); 
+  fMultRecRealC = new TH2F("fMultRecRealC"," ",100,0,200,100,0,200); 
+  fMultRecTotal = new TH2F("fMultRecTotal"," ",100,0,200,100,0,200); 
+  
+  fPrim = new TH1F("fPrim", " primary",100, 0, 200);
+  fRef  = new TH1F("fRef", " from TR ",100, 0, 200);
+  fRec  = new TH1F("fRec", " in ESD ",100, 0, 200);
+   
+  fListOfHistos->Add(fOrA);
+  fListOfHistos->Add(fOrC);
+  fListOfHistos->Add(fMean);
+  fListOfHistos->Add(fVertex);
+  fListOfHistos->Add(fAmp);
+  fListOfHistos->Add(fTime);
+  fListOfHistos->Add(fTotalMult);
+  fListOfHistos->Add(fMultRecRealA);
+  fListOfHistos->Add(fMultRecRealC);
+  fListOfHistos->Add(fMultRecTotal);
+  fListOfHistos->Add(fPrim);
+  fListOfHistos->Add(fRef);
+  fListOfHistos->Add(fRec);
+}  
+
+
+void AliT0MultiplicityTask::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+    
+  // MC information
+  AliMCEvent* mcEvent = MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+  
+  
+  Int_t primaryA=0;
+  Int_t primaryC=0;
+  Int_t numPrim = 0;
+  Int_t refT0=0;
+
+  AliStack* stack = mcEvent->Stack();
+  // printf("AliT0MultiplicityAnalysis: Number of tracks on stack %5d\n", stack->GetNtrack());
+  Int_t nTracks  = mcEvent->GetNumberOfTracks();
+  for (Int_t ipart=0; ipart<nTracks; ipart++)
+    {
+      //    TParticle* particle = stack->Particle(ipart);
+      AliMCParticle* track = mcEvent->GetTrack(ipart);
+      if (!track) {
+       Printf("ERROR: Could not receive track %d (mc loop)", ipart);
+       continue;
+      }
+      Int_t label = track->GetLabel();
+      if(stack->IsPhysicalPrimary(label) == kFALSE)
+       continue;
+      
+      numPrim++;
+      Double_t eta=track->Eta();
+      
+      if(eta<-2.97 && eta>-3.28)  primaryC++;
+      if (eta >4.61 && eta<4.92)  primaryA++;
+    
+       // Loop over Track References
+       AliTrackReference* trackRef = 0;
+
+      for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
+       trackRef = track->GetTrackReference(iTrackRef);
+       if(trackRef) {
+         Int_t detectorId = trackRef->DetectorId(); 
+         if (detectorId == 7)  refT0++;
+       }      
+      }
+    }
+  fTotalMult->Fill(numPrim);
+  fPrim->Fill(primaryC+primaryA);
+  fRef->Fill(refT0);
+  //     printf (" tracks %i primaries %d \n",nTracks, numPrim);
+
+  // ESD information  
+  AliVEvent* event = InputEvent();
+  if (!event) {
+    Printf("ERROR: Could not retrieve event");
+    return;
+  }
+    
+  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
+
+  Bool_t eventTriggered = false;
+  ULong64_t triggerMask = esd->GetTriggerMask();
+  // definitions from p-p.cfg
+  ULong64_t spdFO = (1 << 14);
+  ULong64_t v0left = (1 << 11);
+  ULong64_t v0right = (1 << 12);
+
+  if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
+    eventTriggered == true;
+    //= AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
+  //  printf("!!!!! eventTriggered %i",eventTriggered);
+  // if(!eventTriggered) return;
+   Double_t besttimeC=99999.;
+  Double_t besttimeA=99999.;
+  
+//  Float_t coefA = 0.891466;
+//  Float_t coefC = 0.922772;
+  Float_t coefA = 1;
+  Float_t coefC = 1;
+  Float_t sumampA=0;
+  Float_t sumampC=0;
+  Int_t highA=0, highC=0;
+  const Double_t *amp = esd->GetT0amplitude();
+  const Double_t *time = esd->GetT0time();
+  Float_t vertex = esd->GetT0zVertex();
+  if(vertex<999) fVertex->Fill(vertex);
+  Float_t   mean = esd->GetT0();
+  fMean->Fill(mean);
+  
+  for (Int_t i=0; i<12; i++) {
+    sumampC += amp[i];  
+    fTime->Fill(amp[i],time[i]); 
+    fAmp->Fill(amp[i]); 
+    if(time[i]<besttimeC && time[i]>0) besttimeC=time[i]; //timeC
+  }
+  fOrC->Fill(besttimeC);
+  
+  for (Int_t i=12; i<24; i++){
+    sumampA += amp[i];  
+    fAmp->Fill(amp[i]); 
+    fTime->Fill(time[i]);
+    if(time[i]<besttimeA && time[i]>0) besttimeA=time[i]; //timeC
+  }
+
+  fOrA->Fill(besttimeA);
+  
+  fMultRecRealA ->Fill(primaryA, sumampA*coefA);
+  fMultRecRealC ->Fill(primaryC, sumampC*coefC);
+  fMultRecTotal->Fill(numPrim, sumampA*coefA +  sumampC*coefC);
+  fRec->Fill(Int_t(sumampA+sumampC));
+  // Post output data.
+  PostData(1, fListOfHistos);
+}      
+
+
+void AliT0MultiplicityTask::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query 
+  //  printf(" AliT0MultiplicityTask::Terminate ");
+  
+  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
+  if (!fListOfHistos) {
+    Printf("ERROR: fListOfHistos not available");
+    return;
+  }
+  //  printf(" before cast ");
+  fOrA = dynamic_cast<TH1F*>(fListOfHistos->At(0));
+  fOrC = dynamic_cast<TH1F*>(fListOfHistos->At(1));
+  fMean  = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
+  fVertex  = dynamic_cast<TH1F*>(fListOfHistos->At(3));
+  fAmp = dynamic_cast<TH1F*>(fListOfHistos->At(4));
+  fTime = dynamic_cast<TH1F*>(fListOfHistos->At(5));
+  fTotalMult = dynamic_cast<TH1F*>(fListOfHistos->At(6));
+  fMultRecTotal = dynamic_cast<TH2F*>(fListOfHistos->At(7));
+  fMultRecRealA = dynamic_cast<TH2F*>(fListOfHistos->At(8));
+  fMultRecRealC= dynamic_cast<TH2F*>(fListOfHistos->At(9));
+  fPrim= dynamic_cast<TH1F*>(fListOfHistos->At(10));
+  fRef= dynamic_cast<TH1F*>(fListOfHistos->At(11));
+  fRec= dynamic_cast<TH1F*>(fListOfHistos->At(12));
+
+  TFile fc("MultHist.root\n", "RECREATE");
+  //  printf(" File MultHist.root recreated\n");
+  fOrC->Write();
+  //  printf("1\n");
+  fOrA->Write();
+  //  printf("2\n");
+  fMean->Write();
+  //  printf("3\n");
+  fVertex->Write();
+  //  printf("4\n");
+  fAmp->Write();
+  //  printf("5\n");
+  fTime->Write();
+ fTotalMult->Write();
+ //  printf("6\n");
+ fMultRecTotal->Write();
+ //  printf("7\n");
+ fMultRecRealA->Write();
+ // printf("8\n");
+  fMultRecRealC->Write();
+  // printf("9\n");
+  fPrim->Write();
+  fRef->Write();
+  fRec->Write();
+  
+  fc.Close();
+ printf(" fc.Close()\n");
+
+}
diff --git a/T0/AliT0MultiplicityTask.h b/T0/AliT0MultiplicityTask.h
new file mode 100644 (file)
index 0000000..f95bd89
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef AliT0MultiplicityTask_cxx
+#define AliT0MultiplicityTask_cxx
+
+class TList;
+class TH1F;
+class TH2F;
+class TClonesArray;
+
+
+#include "AliAnalysisTaskSE.h"
+
+class AliT0MultiplicityTask: public AliAnalysisTaskSE
+{
+  public:
+   AliT0MultiplicityTask ();
+   AliT0MultiplicityTask (const char* name);
+    virtual ~AliT0MultiplicityTask() {}
+    
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *);
+  
+  private:
+    TList* fListOfHistos;
+    TH1F* fOrA; // Or A
+    TH1F* fOrC; // Or C
+    TH1F* fMean; // multiplicity histogram
+    TH1F* fVertex; // multiplicity histogram
+    TH1F* fTime;
+    TH1F* fAmp;
+    TH1F* fTotalMult;
+    TH2F* fMultRecTotal;
+    TH2F* fMultRecRealA;
+    TH2F* fMultRecRealC;
+    TH1F*  fPrim;
+    TH1F* fRef;
+    TH1F* fRec;
+    
+  AliT0MultiplicityTask (const AliT0MultiplicityTask&); // not implemented
+  AliT0MultiplicityTask& operator=(const AliT0MultiplicityTask&); // not implemented
+
+    ClassDef(AliT0MultiplicityTask, 1); // example of analysis 
+};
+
+#endif
diff --git a/T0/runProofT0analysis.C b/T0/runProofT0analysis.C
new file mode 100644 (file)
index 0000000..57e9eb2
--- /dev/null
@@ -0,0 +1,52 @@
+void runProofT0analysis(const char * dataset = "/COMMON/COMMON/LHC09a4_10TeV_200k#esdTree",Long64_t nentries=20000, Long64_t firstentry=0)
+{
+// Connect to Proof
+  TProof::Open("proof://alla@alicecaf.cern.ch"); 
+  //TProof::Open("lxb6046");
+
+  // Upload and enable packages: please use the correct version!
+  gProof->UploadPackage("AF-v4-16");
+  gProof->EnablePackage("AF-v4-16");
+  gProof->ShowDataSets();      
+  // Create the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("AliT0MultiplicityTask");
+
+  AliVEventHandler* esdH = new AliESDInputHandler();
+  mgr->SetInputEventHandler(esdH);
+
+  // Enable MC event handler
+  AliVEventHandler* handler = new AliMCEventHandler;
+  mgr->SetMCtruthEventHandler(handler);
+
+  // Create task
+  //  gProof->Load("AliMCComparisonTrack.cxx++g");
+  gProof->Load("AliT0MultiplicityTask.cxx++g");
+  AliAnalysisTask *task = new AliT0MultiplicityTask("AliT0MultiplicityTask");
+
+  // Add task
+  mgr->AddTask(task);
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput = 
+    mgr->CreateContainer("cchain", TChain::Class(), AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput = 
+    mgr->CreateContainer("coutput", TList::Class(), 
+    AliAnalysisManager::kOutputContainer, "MultHist.root");
+
+  // Connect input/output
+  mgr->ConnectInput(task, 0, cinput);
+  mgr->ConnectOutput(task, 1, coutput);
+
+
+  // Enable debug printouts
+  mgr->SetDebugLevel(3);
+
+  if (!mgr->InitAnalysis())
+    return;
+
+  mgr->PrintStatus();
+
+   mgr->StartAnalysis("proof",dataset,nentries,firstentry);
+}
+