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