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