Fixing the code for the analysis train
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 Mar 2011 21:45:23 +0000 (21:45 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 Mar 2011 21:45:23 +0000 (21:45 +0000)
PWG2/EBYE/Fluctuations/AliEbyEFluctuationAnalysisTaskTrain.cxx

index e9204ea..6ed5a92 100644 (file)
-#include "TChain.h"
-#include "TTree.h"
-#include "TH1F.h"
-#include "TH2F.h"
-#include "TCanvas.h"
-#include "TParticle.h"
-#include "TObjArray.h"
-
-#include "AliAnalysisTask.h"
-#include "AliAnalysisManager.h"
-
-#include "AliStack.h"
-#include "AliMCEvent.h"
-#include "AliESDEvent.h"
-#include "AliESDInputHandler.h"
-//#include "AliESDtrackCuts.h"
-#include "AliCentrality.h"
-
-#include "AliEbyEFluctuationAnalysisTaskTrain.h"
-
-// Event by event charge fluctuation analysis
-// Authors: Satyajit Jena and Panos Cristakoglou
-// 
-
-ClassImp(AliEbyEFluctuationAnalysisTaskTrain)
-
-//________________________________________________________________________
-AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name) 
-  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), 
-    fHistEventStats(0), fHistCentrality(0),
-    fHistNMult(0), fHistNPlusNMinus(0),
-    fHistNMultMC(0), fHistNPlusNMinusMC(0),
-    fAnalysisType("ESD"), fAnalysisMode("TPC"),
-    fCentralityEstimator("V0M"), 
-    fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),
-    fVxMax(3.), fVyMax(3.), fVzMax(10.),
-    fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), 
-    fMaxDCAxy(3.0), fMaxDCAz(3.0) {
-  // Constructor
-  
-  // Define input and output slots here
-  // Input slot #0 works with a TChain
-  DefineInput(0, TChain::Class());
-  // Output slot #0 id reserved by the base class for AOD
-  // Output slot #1 writes into a TH1 container
-  DefineOutput(1, TList::Class());
-}
-
-//________________________________________________________________________
-void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {
-  // Create histograms
-  // Called once
-
-  fOutputList = new TList();
-  fOutputList->SetOwner();
-
-  //Event stats.
-  TString gCutName[4] = {"Total","Offline trigger",
-                         "Vertex","Analyzed"};
-  fHistEventStats = new TH1F("fHistEventStats",
-                             "Event statistics;;N_{events}",
-                             4,0.5,4.5);
-  for(Int_t i = 1; i <= 4; i++)
-    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
-  fOutputList->Add(fHistEventStats);
-
-  //ESD analysis
-  if(fAnalysisType == "ESD") {
-    fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
-                              20,0.5,20.5);
-    fOutputList->Add(fHistCentrality);
-    
-    TString histName;
-    histName = "fHistNMult";
-    fHistNMult = new TH1F(histName.Data(), 
-                         ";N_{mult.}",
-                         500,0,3000);
-    fOutputList->Add(fHistNMult);
-
-    histName = "fHistNPlusNMinus";
-    fHistNPlusNMinus = new TH2F(histName.Data(), 
-                               ";N_{+};N_{-}",
-                               2000,0.5,2000.5,2000,0.5,2000.5);  
-    fOutputList->Add(fHistNPlusNMinus);
-  }//ESD analysis
-  else if(fAnalysisType == "MC") {
-    TString histName = "fHistNMultMC";
-    fHistNMultMC = new TH1F(histName.Data(), 
-                           ";N_{mult.}",
-                           600,0,6000);
-    fOutputList->Add(fHistNMultMC);
-    
-    histName = "fHistNPlusNMinusMC";
-    fHistNPlusNMinusMC = new TH2F(histName.Data(), 
-                                 ";N_{+};N_{-}",
-                                 3000,0.5,3000.5,3000,0.5,3000.5);  
-    fOutputList->Add(fHistNPlusNMinusMC);
-  }//MC analysis
-
-  PostData(1, fOutputList);
-}
-
-//________________________________________________________________________
-void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {
-  // Main loop
-  // Called for each event
-  
-  Int_t nPlus = 0, nMinus = 0;
-
-  // Post output data.
-  //ESD analysis
-  if(fAnalysisType == "ESD") {
-    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
-    if (!fESD) {
-      printf("ERROR: fESD not available\n");
-      return;
-    }
-  
-    fHistEventStats->Fill(1); //all events
-    
-    //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
-    //if(isSelected) {
-      fHistEventStats->Fill(2); //triggered + centrality
-      //Printf("Event accepted");
-
-      //Centrality stuff
-      AliCentrality *centrality = fESD->GetCentrality();
-      Int_t nCentrality = 0;
-      
-      if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
-                                             fCentralityPercentileMax,
-                                             fCentralityEstimator.Data())) {
-       if(fAnalysisMode = "TPC") {
-         const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
-         if(vertex) {
-           if(vertex->GetNContributors() > 0) {
-             if(vertex->GetZRes() != 0) {
-               fHistEventStats->Fill(3); //events with a proper vertex
-               if(TMath::Abs(vertex->GetXv()) < fVxMax) {
-                 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
-                   if(TMath::Abs(vertex->GetZv()) < fVzMax) {
-                     fHistEventStats->Fill(4); //analyzed events
-                     
-                     Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
-                            centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
-                            nCentrality,fESD->GetNumberOfTracks());
-                     
-                     Int_t nAcceptedTracks = 0;
-                     //TObjArray *gTrackArray = 0;
-                     //if(fESDtrackCuts)
-                     //gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
-                     //if(gTrackArray) {
-                     //nAcceptedTracks = gTrackArray->GetEntries();
-                     nAcceptedTracks = fESD->GetNumberOfTracks();
-                     AliESDtrack* track = 0;
-
-                     Float_t dcaXY = 0.0, dcaZ = 0.0;
-
-                     // Track loop to fill a pT spectrum
-                     for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) {
-                       //track = dynamic_cast<AliESDtrack *>(gTrackArray->At(iTracks));
-                       track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));
-                       if (!track) {
-                         printf("ERROR: Could not receive track %d\n", iTracks);
-                         continue;
-                       }
-
-                       AliESDtrack *tpcOnlyTrack = new AliESDtrack();
-                       
-                       // only true if we have a tpc track
-                       if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
-                         delete tpcOnlyTrack;
-                         return;
-                       }
-                       
-                       tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);
-                       
-                       //==============================================================//
-                       Int_t nClustersITS = track->GetITSclusters(0x0);
-                       Int_t nClustersTPC = track->GetTPCclusters(0x0);
-
-                       Float_t chi2PerClusterITS = -1;
-                       if (nClustersITS!=0)
-                         chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
-                       Float_t chi2PerClusterTPC = -1;
-                       if (nClustersTPC!=0)
-                         chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
-                       
-                       //Printf("TPC clusters: %d - %f === dca: %lf %lf",nClustersTPC, chi2PerClusterTPC,dcaXY,dcaZ);
-
-                       if(nClustersTPC < fMinNumberOfTPCClusters) continue;
-                       if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;
-                       if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;
-                       if(TMath::Abs(dcaZ) > fMaxDCAz) continue;
-                       //==============================================================//
-
-                       //if(fESDtrackCuts) {  
-                       //AliESDtrack *tpcOnlyTrack = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks);  
-                       //if(!tpcOnlyTrack) continue;  
-                           
-                       //if(!fESDtrackCuts->AcceptTrack(tpcOnlyTrack)) continue;  
-                       //}//ESD track cuts object
-
-                       Short_t gCharge = tpcOnlyTrack->Charge();
-                       if(gCharge > 0) nPlus += 1;
-                       if(gCharge < 0) nMinus += 1;
-                     }//track loop
-                     //}//TObjArray valid object
-                     //if((nCentrality >= 1)&&(nCentrality <= 20)) {
-                     cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;
-                     //Printf("===NPlus: %d - NMinus: %d===",nPlus,nMinus);
-                     fHistCentrality->Fill(nCentrality);
-                     fHistNPlusNMinus->Fill(nPlus,nMinus);
-                     fHistNMult->Fill(nPlus+nMinus);
-                     //}
-                   }//Vz cut
-                 }//Vy cut
-               }//Vx cut
-             }//Vz resolution
-           }//number of contributors
-         }//valid vertex
-       }//TPC analysis mode
-      }//centrality
-       //}//physics selection
-  }//ESD analysis level
-  
-  //MC analysis
-  if(fAnalysisType == "MC") {
-    AliMCEvent* mcEvent = MCEvent();
-    if (!mcEvent) {
-      Printf("ERROR: Could not retrieve MC event");
-      return;
-    }
-    AliStack* stack = mcEvent->Stack();
-    if (!stack) {
-      Printf("ERROR: Could not retrieve MC stack");
-      return;
-    }
-    
-    fHistEventStats->Fill(1); 
-    fHistEventStats->Fill(2); 
-    fHistEventStats->Fill(3); 
-    fHistEventStats->Fill(4); //analyzed events
-    for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
-      AliVParticle* track = mcEvent->GetTrack(iTracks);
-      if (!track) {
-       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
-       continue;
-      }
-      
-      if(!stack->IsPhysicalPrimary(iTracks)) continue;
-      if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
-      if(TMath::Abs(track->Eta()) > 0.8) continue;
-      
-      Short_t gCharge = track->Charge();
-      if(gCharge > 0) nPlus += 1;
-      if(gCharge < 0) nMinus += 1;
-    }//particle loop
-    fHistNPlusNMinusMC->Fill(nPlus,nMinus);
-    fHistNMultMC->Fill(nPlus+nMinus);
-    
-  }//MC analysis level
-
-
-}      
-
-//________________________________________________________________________
-void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {
-  // Draw result to the screen
-  // Called once at the end of the query
-
-  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
-  if (!fOutputList) {
-    printf("ERROR: Output list not available\n");
-    return;
-  }
-}
+#include "TChain.h"\r
+#include "TTree.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TCanvas.h"\r
+#include "TParticle.h"\r
+#include "TObjArray.h"\r
+\r
+#include "AliAnalysisTask.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliStack.h"\r
+#include "AliMCEvent.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDInputHandler.h"\r
+//#include "AliESDtrackCuts.h"\r
+#include "AliCentrality.h"\r
+\r
+#include "AliEbyEFluctuationAnalysisTaskTrain.h"\r
+\r
+// Event by event charge fluctuation analysis\r
+// Authors: Satyajit Jena and Panos Cristakoglou\r
+// \r
+\r
+ClassImp(AliEbyEFluctuationAnalysisTaskTrain)\r
+\r
+//________________________________________________________________________\r
+AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name) \r
+  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), \r
+    fHistEventStats(0), fHistCentrality(0),\r
+    fHistNMult(0), fHistNPlusNMinus(0),\r
+    fHistNMultMC(0), fHistNPlusNMinusMC(0),\r
+    fAnalysisType("ESD"), fAnalysisMode("TPC"),\r
+    fCentralityEstimator("V0M"), \r
+    fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),\r
+    fVxMax(3.), fVyMax(3.), fVzMax(10.),\r
+    fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), \r
+    fMaxDCAxy(3.0), fMaxDCAz(3.0) {\r
+  // Constructor\r
+  \r
+  // Define input and output slots here\r
+  // Input slot #0 works with a TChain\r
+  DefineInput(0, TChain::Class());\r
+  // Output slot #0 id reserved by the base class for AOD\r
+  // Output slot #1 writes into a TH1 container\r
+  DefineOutput(1, TList::Class());\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {\r
+  // Create histograms\r
+  // Called once\r
+\r
+  fOutputList = new TList();\r
+  fOutputList->SetOwner();\r
+\r
+  //Event stats.\r
+  TString gCutName[4] = {"Total","Offline trigger",\r
+                         "Vertex","Analyzed"};\r
+  fHistEventStats = new TH1F("fHistEventStats",\r
+                             "Event statistics;;N_{events}",\r
+                             4,0.5,4.5);\r
+  for(Int_t i = 1; i <= 4; i++)\r
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
+  fOutputList->Add(fHistEventStats);\r
+\r
+  //ESD analysis\r
+  if(fAnalysisType.CompareTo("ESD") == 0 ) {\r
+    fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",\r
+                              20,0.5,20.5);\r
+    fOutputList->Add(fHistCentrality);\r
+    \r
+    TString histName;\r
+    histName = "fHistNMult";\r
+    fHistNMult = new TH1F(histName.Data(), \r
+                         ";N_{mult.}",\r
+                         800,0,4000);\r
+    fOutputList->Add(fHistNMult);\r
+\r
+    histName = "fHistNPlusNMinus";\r
+    fHistNPlusNMinus = new TH2F(histName.Data(), \r
+                               ";N_{+};N_{-}",\r
+                               2000,0.5,2000.5,2000,0.5,2000.5);  \r
+    fOutputList->Add(fHistNPlusNMinus);\r
+  }//ESD analysis\r
+  else if(fAnalysisType.CompareTo("MC") == 0) {\r
+    TString histName = "fHistNMultMC";\r
+    fHistNMultMC = new TH1F(histName.Data(), \r
+                           ";N_{mult.}",\r
+                           800,0,8000);\r
+    fOutputList->Add(fHistNMultMC);\r
+    \r
+    histName = "fHistNPlusNMinusMC";\r
+    fHistNPlusNMinusMC = new TH2F(histName.Data(), \r
+                                 ";N_{+};N_{-}",\r
+                                 3000,0.5,3000.5,3000,0.5,3000.5);  \r
+    fOutputList->Add(fHistNPlusNMinusMC);\r
+  }//MC analysis\r
+\r
+  PostData(1, fOutputList);\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {\r
+  // Main loop\r
+  // Called for each event\r
+  \r
+  Int_t nPlus = 0, nMinus = 0;\r
+\r
+  // Post output data.\r
+  //ESD analysis\r
+  if(fAnalysisType.CompareTo("ESD") == 0) {\r
+    fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
+    if (!fESD) {\r
+      printf("ERROR: fESD not available\n");\r
+      return;\r
+    }\r
+  \r
+    fHistEventStats->Fill(1); //all events\r
+    \r
+    //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
+    //if(isSelected) {\r
+      fHistEventStats->Fill(2); //triggered + centrality\r
+      \r
+\r
+      //Centrality stuff\r
+      AliCentrality *centrality = fESD->GetCentrality();\r
+      Int_t nCentrality = 0;\r
+      \r
+      if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
+                                             fCentralityPercentileMax,\r
+                                             fCentralityEstimator.Data())) {\r
+       if(fAnalysisMode.CompareTo("TPC") == 0 ) {\r
+         const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();\r
+         if(vertex) {\r
+           if(vertex->GetNContributors() > 0) {\r
+             if(vertex->GetZRes() != 0) {\r
+               fHistEventStats->Fill(3); //events with a proper vertex\r
+               if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
+                 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
+                   if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
+                     fHistEventStats->Fill(4); //analyzed events\r
+                     \r
+                     Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",\r
+                             centrality->GetCentralityPercentile(fCentralityEstimator.Data()),\r
+                             nCentrality,fESD->GetNumberOfTracks());\r
+                     \r
+\r
+                     Int_t nAcceptedTracks = 0;\r
+                     \r
+                     nAcceptedTracks = fESD->GetNumberOfTracks();\r
+                     AliESDtrack* track = 0;\r
+                     \r
+                     Float_t dcaXY = 0.0, dcaZ = 0.0;\r
+                     \r
+                     for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) \r
+                       {\r
+                         track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));\r
+                         if (!track) {\r
+                           printf("ERROR: Could not receive track %d\n", iTracks);\r
+                           continue;\r
+                         }\r
+                         \r
+                         AliESDtrack *tpcOnlyTrack = new AliESDtrack();\r
+                         \r
+                         if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {\r
+                           delete tpcOnlyTrack;\r
+                           continue;\r
+                         }\r
+                         \r
+                         tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);\r
+                                                 \r
+                         Int_t nClustersITS = track->GetITSclusters(0x0);\r
+                         Int_t nClustersTPC = track->GetTPCclusters(0x0);\r
+                         \r
+                         Float_t chi2PerClusterITS = -1;\r
+                         if (nClustersITS!=0)\r
+                           chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);\r
+                         Float_t chi2PerClusterTPC = -1;\r
+                         if (nClustersTPC!=0)\r
+                           chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);\r
+                         \r
+                         //    \r
+                         \r
+                         if(nClustersTPC < fMinNumberOfTPCClusters) continue;\r
+                         if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;\r
+                         if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;\r
+                         if(TMath::Abs(dcaZ) > fMaxDCAz) continue;\r
+                         \r
+                         Short_t gCharge = tpcOnlyTrack->Charge();\r
+                         if(gCharge > 0) nPlus += 1;\r
+                         if(gCharge < 0) nMinus += 1;\r
+                         delete tpcOnlyTrack;\r
+                       } // track loop\r
+                     \r
+                     if( nPlus == 0 && nMinus == 0) return;\r
+\r
+                     //  cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;\r
+                     fHistCentrality->Fill(nCentrality);\r
+                     fHistNPlusNMinus->Fill(nPlus,nMinus);\r
+                     fHistNMult->Fill(nPlus+nMinus);\r
+                     \r
+                     \r
+                   }//Vz cut\r
+                 }//Vy cut\r
+               }//Vx cut\r
+             }//Vz resolution\r
+           }//number of contributors\r
+         }//valid vertex\r
+       }//TPC analysis mode\r
+      }//centrality\r
+      //}//physics selection\r
+  }//ESD analysis level\r
+  \r
+      //MC analysis\r
+  if(fAnalysisType.CompareTo("MC") == 0 ) {\r
+    AliMCEvent* mcEvent = MCEvent();\r
+    if (!mcEvent) {\r
+      Printf("ERROR: Could not retrieve MC event");\r
+      return;\r
+    }\r
+    AliStack* stack = mcEvent->Stack();\r
+    if (!stack) {\r
+      Printf("ERROR: Could not retrieve MC stack");\r
+      return;\r
+    }\r
+    \r
+    fHistEventStats->Fill(1); \r
+    fHistEventStats->Fill(2); \r
+    fHistEventStats->Fill(3); \r
+    fHistEventStats->Fill(4); //analyzed events\r
+    for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
+      AliVParticle* track = mcEvent->GetTrack(iTracks);\r
+      if (!track) {\r
+       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);\r
+       continue;\r
+      }\r
+      \r
+      if(!stack->IsPhysicalPrimary(iTracks)) continue;\r
+      if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;\r
+      if(TMath::Abs(track->Eta()) > 0.8) continue;\r
+      \r
+      Short_t gCharge = track->Charge();\r
+      if(gCharge > 0) nPlus += 1;\r
+      if(gCharge < 0) nMinus += 1;\r
+    }//particle loop\r
+    fHistNPlusNMinusMC->Fill(nPlus,nMinus);\r
+    fHistNMultMC->Fill(nPlus+nMinus);\r
+    \r
+  }//MC analysis level\r
+\r
+\r
+}      \r
+\r
+//________________________________________________________________________\r
+void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {\r
+  // Draw result to the screen\r
+  // Called once at the end of the query\r
+\r
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));\r
+  if (!fOutputList) {\r
+    printf("ERROR: Output list not available\n");\r
+    return;\r
+  }\r
+}\r