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), fUseFriends(kFALSE), fHistNTPCCl(0), fHistNESDtracks(0),   fHistNESDfriendtracks(0)
  30
  31 {
  32   // Constructor
  33
  34   // Define input and output slots here
  35   // Input slot #0 works with a TChain
  36   DefineInput(0, TChain::Class());
  37   // Output slot #0 writes into a TH1 container
  38   DefineOutput(0, TList::Class());
  39 }
  40
  41 //________________________________________________________________________
  42 void AliAnalysisTaskPt::ConnectInputData(Option_t *)
  43 {
  44   // Connect ESD or AOD here
  45   // Called once
  46
  47   printf("----> AliAnalysisTaskPt::ConnectInputData\n");
  48   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  49   if (!tree) {
  50     Printf("ERROR: Could not read chain from input slot 0");
  51   } else {
  52     // Disable all branches and enable only the needed ones
  53     // The next two lines are different when data produced as AliESDEvent is read
  54     /*
  55     tree->SetBranchStatus("*", kFALSE);
  56     tree->SetBranchStatus("fTracks.*", kTRUE);
  57     */
  58
  59     AliVEventHandler *esdH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  60
  61     if (!esdH) {
  62       Printf("ERROR: Could not get ESDInputHandler");
  63     } else {
  64       Printf("----> AliAnalysisTaskPt::ConnectInputData Getting the event from handler %p", esdH);
  65       //fESD = dynamic_cast<AliESDEvent*>(esdH->GetEvent());
  66       fESD = esdH->GetVVEvent();
  67       if (fUseFriends){
  68         fESDfriend = esdH->GetVVFriendEvent();
  69       }
  70     }
  71     if (!fESD) {
  72       Printf("ERROR, no ESD event");
  73     }
  74     if (fUseFriends && !fESDfriend) {
  75       Printf("ERROR, no ESD friend");
  76     }
  77   }
  78
  79   Printf("fESD = %p, fESDfriend = %p", fESD, fESDfriend);
  80   printf("<---- AliAnalysisTaskPt::ConnectInputData\n");
  81 }
  82
  83 //________________________________________________________________________
  84 void AliAnalysisTaskPt::CreateOutputObjects()
  85 {
  86   // Create histograms
  87   // Called once
  88
  89   fListOut = new TList();
  90   fListOut->SetOwner();
  91   fListOut->SetName("listHistos");
  92
  93   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
  94   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  95   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  96   fHistPt->SetMarkerStyle(kFullCircle);
  97
  98   fHistQ = new TH1F("fHistQ", "TPC clusters: Q distribution", 1000, 0, 10000);
  99   fHistQ->GetXaxis()->SetTitle("Q");
 100   fHistQ->GetYaxis()->SetTitle("dN/dQ");
 101   fHistQ->SetMarkerStyle(kFullCircle);
 102
 103   fHistNTPCCl = new TH1F("fHistNTPCCl", "Number of TPC clusters", 160, -0.5, 159.5);
 104   fHistNTPCCl->GetXaxis()->SetTitle("n. TPC Cl.");
 105   fHistNTPCCl->GetYaxis()->SetTitle("dN/d(n. TPC Cl)");
 106   fHistNTPCCl->SetMarkerStyle(kFullCircle);
 107
 108   fHistNESDtracks = new TH1F("fHistNESDtracks", "Number of ESD friend tracks", 1000, -0.5, 999.5);
 109   fHistNESDtracks->GetXaxis()->SetTitle("n. ESD friend tracks");
 110   fHistNESDtracks->GetYaxis()->SetTitle("dN/d(n. ESD friend tracks)");
 111   fHistNESDtracks->SetMarkerStyle(kFullCircle);
 112
 113   fHistNESDfriendtracks = new TH1F("fHistNESDfriendtracks", "Number of ESD tracks", 1000, -0.5, 999.5);
 114   fHistNESDfriendtracks->GetXaxis()->SetTitle("n. ESD tracks");
 115   fHistNESDfriendtracks->GetYaxis()->SetTitle("dN/d(n. ESD tracks)");
 116   fHistNESDfriendtracks->SetMarkerStyle(kFullCircle);
 117
 118   fListOut->Add(fHistPt);
 119   fListOut->Add(fHistQ);
 120   fListOut->Add(fHistNTPCCl);
 121   fListOut->Add(fHistNESDtracks);
 122   fListOut->Add(fHistNESDfriendtracks);
 123
 124   PostData(0, fListOut);
 125
 126   fCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
 127 }
 128
 129 //________________________________________________________________________
 130 void AliAnalysisTaskPt::Exec(Option_t *)
 131 {
 132   // Main loop
 133   // Called for each event
 134
 135   if (!fESD) {
 136     Printf("ERROR: fESD not available");
 137     return;
 138   }
 139   if (!fESDfriend) {
 140     Printf("ERROR: fESDfriend not available");
 141       if (fUseFriends){
 142         return;
 143       }
 144   }
 145
 146   Int_t nESDtracks = fESD->GetNumberOfTracks();
 147   Int_t nESDfriendtracks = 0;
 148   if (fUseFriends) nESDfriendtracks = fESDfriend->GetNumberOfTracks();
 149   Printf("There are %d tracks in this event", nESDtracks);
 150   Printf("... and there are %d friends in this event", nESDfriendtracks);
 151
 152   fHistNESDtracks->Fill(nESDtracks);
 153   fHistNESDfriendtracks->Fill(nESDfriendtracks);
 154
 155   // Track loop to fill a pT spectrum
 156   for (Int_t iTracks = 0; iTracks < nESDtracks; iTracks++) {
 157     Printf("Checking track %d: Note that with Flat, the GetTrack is not yet implemented!!!", iTracks);
 158     AliVVtrack* track = fESD->GetVVTrack(iTracks);
 159     if (!track) {
 160       Printf("ERROR: Could not receive track %d", iTracks);
 161       continue;
 162     }
 163     Printf("track %d has pt = %f", iTracks, track->GetPt());
 164     fHistPt->Fill(track->GetPt());
 165     fHistNTPCCl->Fill(track->GetTPCNcls());
 166   } //track loop
 167
 168
 169   if (fUseFriends){
 170     // Friend Track loop
 171     for (Int_t iFriend = 0; iFriend < nESDfriendtracks; iFriend++) {
 172       //Printf("Getting friend %d", iFriend);
 173       const AliVVfriendTrack* friendTrack = fESDfriend->GetTrack(iFriend);
 174       if (!friendTrack) {
 175         Printf("ERROR: Could not receive track %d", iFriend);
 176         continue;
 177       }
 178
 179       AliTPCseed seed;
 180       Int_t err = friendTrack->GetTPCseed( seed );
 181       if( err==0 ){
 182         Printf("Found TPC seed" );
 183         for (Int_t irow = 0; irow < 160; irow++){
 184           AliTPCclusterMI* cluMI = seed.GetClusterPointer(irow);
 185           if (cluMI){
 186             Printf("Found cluster at row %d", irow);
 187             Float_t q = cluMI->GetQ();
 188             Printf("Q = %f", q);
 189             fHistQ->Fill(q);
 190           }
 191           else {
 192             Printf("Row %d does not contain clusters", irow);
 193           }
 194         }
 195       }
 196       else {
 197         //Printf("Schade... seed is %p", seed);
 198       }
 199     }
 200   }
 201
 202   // Post output data.
 203   PostData(0, fListOut);
 204   fEv++;
 205 }
 206
 207 //________________________________________________________________________
 208 void AliAnalysisTaskPt::Terminate(Option_t *)
 209 {
 210   // Draw result to the screen
 211   // Called once at the end of the query
 212
 213   Printf("Terminate called: fESD = %p", fESD);
 214
 215   fListOut = dynamic_cast<TList*> (GetOutputData(0));
 216
 217   if (fListOut) {
 218     fHistPt = dynamic_cast<TH1F*>(fListOut->FindObject("fHistPt"));
 219     if (!fHistPt) {
 220       Printf("ERROR: fHistPt not available");
 221       return;
 222     }
 223
 224     TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
 225     c1->cd(1)->SetLogy();
 226     fHistPt->DrawCopy("E");
 227   }
 228   else {
 229     Printf("In Terminate: no TList found");
 230   }
 231
 232 }