ANALYSIS/AliAnalysisTaskPt.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TCanvas.h"
   5 #include "TObjArray.h"
   6
   7 #include "AliAnalysisTask.h"
   8 #include "AliAnalysisManager.h"
   9
  10 #include "AliESDEvent.h"
  11 #include "AliVVevent.h"
  12 #include "AliVVtrack.h"
  13 #include "AliESDtrackCuts.h"
  14 #include "AliVEventHandler.h"
  15 #include "../TPC/Rec/AliTPCseed.h"
  16 #include "../TPC/Rec/AliTPCclusterMI.h"
  17
  18 #include "AliAnalysisTaskPt.h"
  19
  20 // example of an analysis task creating a p_t spectrum
  21 // Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing
  22
  23 ClassImp(AliAnalysisTaskPt)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskPt::AliAnalysisTaskPt(const char *name)
  27 : AliAnalysisTask(name, ""), fESD(0), fESDfriend(0), fHistPt(0), fCuts(0), fEv(0), fHistQ(0), fListOut(0)
  28 {
  29   // Constructor
  30
  31   // Define input and output slots here
  32   // Input slot #0 works with a TChain
  33   DefineInput(0, TChain::Class());
  34   // Output slot #0 writes into a TH1 container
  35   DefineOutput(0, TList::Class());
  36 }
  37
  38 //________________________________________________________________________
  39 void AliAnalysisTaskPt::ConnectInputData(Option_t *)
  40 {
  41   // Connect ESD or AOD here
  42   // Called once
  43
  44   printf("----> AliAnalysisTaskPt::ConnectInputData\n");
  45   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  46   if (!tree) {
  47     Printf("ERROR: Could not read chain from input slot 0");
  48   } else {
  49     // Disable all branches and enable only the needed ones
  50     // The next two lines are different when data produced as AliESDEvent is read
  51     /*
  52     tree->SetBranchStatus("*", kFALSE);
  53     tree->SetBranchStatus("fTracks.*", kTRUE);
  54     */
  55
  56     AliVEventHandler *esdH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  57
  58     if (!esdH) {
  59       Printf("ERROR: Could not get ESDInputHandler");
  60     } else {
  61       Printf("----> AliAnalysisTaskPt::ConnectInputData Getting the event from handler %p", esdH);
  62       //fESD = dynamic_cast<AliESDEvent*>(esdH->GetEvent());
  63       fESD = esdH->GetEvent();
  64       fESDfriend = esdH->GetFriendEvent();
  65     }
  66     if (!fESD) {
  67       Printf("ERROR, no ESD event");
  68     }
  69     if (!fESDfriend) {
  70       Printf("ERROR, no ESD friend");
  71     }
  72   }
  73
  74   Printf("fESD = %p, fESDfriend = %p", fESD, fESDfriend);
  75   printf("<---- AliAnalysisTaskPt::ConnectInputData\n");
  76 }
  77
  78 //________________________________________________________________________
  79 void AliAnalysisTaskPt::CreateOutputObjects()
  80 {
  81   // Create histograms
  82   // Called once
  83
  84   fListOut = new TList();
  85   fListOut->SetOwner();
  86   fListOut->SetName("listHistos");
  87
  88   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
  89   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  90   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  91   fHistPt->SetMarkerStyle(kFullCircle);
  92
  93   fHistQ = new TH1F("fHistQ", "TPC clusters: Q distribution", 1000, 0, 10000);
  94   fHistQ->GetXaxis()->SetTitle("Q");
  95   fHistQ->GetYaxis()->SetTitle("dN/dQ");
  96   fHistQ->SetMarkerStyle(kFullCircle);
  97
  98   fListOut->Add(fHistPt);
  99   fListOut->Add(fHistQ);
 100
 101   PostData(0, fListOut);
 102
 103   fCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
 104 }
 105
 106 //________________________________________________________________________
 107 void AliAnalysisTaskPt::Exec(Option_t *)
 108 {
 109   // Main loop
 110   // Called for each event
 111
 112   if (!fESD) {
 113     Printf("ERROR: fESD not available");
 114     return;
 115   }
 116   if (!fESDfriend) {
 117     Printf("ERROR: fESDfriend not available");
 118     return;
 119   }
 120
 121   Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 122   Printf("... and there are %d friends in this event", fESDfriend->GetNumberOfTracks());
 123
 124   // Track loop to fill a pT spectrum
 125   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 126     AliVVtrack* track = fESD->GetTrack(iTracks);
 127     if (!track) {
 128       Printf("ERROR: Could not receive track %d", iTracks);
 129       continue;
 130     }
 131
 132     fHistPt->Fill(track->Pt());
 133   } //track loop
 134
 135
 136   // Friend Track loop
 137   for (Int_t iFriend = 0; iFriend < fESDfriend->GetNumberOfTracks(); iFriend++) {
 138     AliESDfriendTrack* friendTrack = fESDfriend->GetTrack(iFriend);
 139     if (!friendTrack) {
 140       Printf("ERROR: Could not receive track %d", iFriend);
 141       continue;
 142     }
 143     TObject* calibObject;
 144     AliTPCseed* seed = NULL;
 145     for (Int_t idx = 0; (calibObject = friendTrack->GetCalibObject(idx)); ++idx) {
 146       Printf(" |Cal %d = %p", idx, calibObject);
 147       if ((seed = dynamic_cast<AliTPCseed*>(calibObject))) {
 148         Printf("Found TPC seed");
 149         for (Int_t irow = 0; irow < 160; irow++){
 150           AliTPCclusterMI* cluMI = seed->GetClusterPointer(irow);
 151           if (cluMI){
 152             Printf("Found cluster at row %d", irow);
 153             Float_t q = cluMI->GetQ();
 154             Printf("Q = %f", q);
 155             fHistQ->Fill(q);
 156           }
 157           else {
 158             Printf("Row %d does not contain clusters", irow);
 159           }
 160         }
 161       }
 162     }
 163   }
 164
 165   // Post output data.
 166   PostData(0, fListOut);
 167   fEv++;
 168 }
 169
 170 //________________________________________________________________________
 171 void AliAnalysisTaskPt::Terminate(Option_t *)
 172 {
 173   // Draw result to the screen
 174   // Called once at the end of the query
 175
 176   Printf("Terminate called: fESD = %p", fESD);
 177
 178   fListOut = dynamic_cast<TList*> (GetOutputData(0));
 179
 180   if (fListOut) {
 181     fHistPt = dynamic_cast<TH1F*>(fListOut->FindObject("fHistPt"));
 182     if (!fHistPt) {
 183       Printf("ERROR: fHistPt not available");
 184       return;
 185     }
 186
 187     TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
 188     c1->cd(1)->SetLogy();
 189     fHistPt->DrawCopy("E");
 190   }
 191   else {
 192     Printf("In Terminate: no TList found");
 193   }
 194
 195 }