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 "AliFlatESDEvent.h"
  12 #include "AliESDtrackCuts.h"
  13 #include "AliVEventHandler.h"
  14 #include "AliTPCseed.h"
  15 #include "AliTPCclusterMI.h"
  16 #include "AliVfriendEvent.h"
  17 #include "AliVfriendTrack.h"
  18
  19 #include "AliAnalysisTaskPt.h"
  20
  21 // example of an analysis task creating a p_t spectrum
  22 // Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing
  23
  24 ClassImp(AliAnalysisTaskPt)
  25
  26 //________________________________________________________________________
  27 AliAnalysisTaskPt::AliAnalysisTaskPt(const char *name)
  28 : AliAnalysisTask(name, ""), fESD(0), fESDfriend(0), fHistPt(0), fCuts(0), fEv(0), fHistQ(0), fListOut(0), fUseFriends(kFALSE), fHistNTPCCl(0), fHistNESDtracks(0),   fHistNESDfriendtracks(0)
  29
  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     TString classInputHandler = esdH->ClassName();
  60
  61     Printf("----> AliAnalysisTaskPt: ClassName of handler = %s", classInputHandler.Data());
  62
  63     if (!esdH) {
  64       Printf("ERROR: Could not get ESDInputHandler");
  65     } else {
  66       Printf("----> AliAnalysisTaskPt::ConnectInputData Getting the event from handler %p", esdH);
  67       //fESD = dynamic_cast<AliFlatESDEvent*>(esdH->GetEvent());
  68       fESD = esdH->GetEvent();
  69       if (fUseFriends){
  70         fESDfriend = esdH->GetVFriendEvent();
  71       }
  72     }
  73     if (!fESD) {
  74       Printf("ERROR, no ESD event");
  75     }
  76     if (fUseFriends && !fESDfriend) {
  77       Printf("ERROR, no ESD friend");
  78     }
  79   }
  80
  81   Printf("fESD = %p, fESDfriend = %p", fESD, fESDfriend);
  82   printf("<---- AliAnalysisTaskPt::ConnectInputData\n");
  83 }
  84
  85 //________________________________________________________________________
  86 void AliAnalysisTaskPt::CreateOutputObjects()
  87 {
  88   // Create histograms
  89   // Called once
  90
  91   fListOut = new TList();
  92   fListOut->SetOwner();
  93   fListOut->SetName("listHistos");
  94
  95   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
  96   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  97   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  98   fHistPt->SetMarkerStyle(kFullCircle);
  99
 100   fHistQ = new TH1F("fHistQ", "TPC clusters: Q distribution", 1000, 0, 10000);
 101   fHistQ->GetXaxis()->SetTitle("Q");
 102   fHistQ->GetYaxis()->SetTitle("dN/dQ");
 103   fHistQ->SetMarkerStyle(kFullCircle);
 104
 105   fHistNTPCCl = new TH1F("fHistNTPCCl", "Number of TPC clusters", 160, -0.5, 159.5);
 106   fHistNTPCCl->GetXaxis()->SetTitle("n. TPC Cl.");
 107   fHistNTPCCl->GetYaxis()->SetTitle("dN/d(n. TPC Cl)");
 108   fHistNTPCCl->SetMarkerStyle(kFullCircle);
 109
 110   fHistNESDtracks = new TH1F("fHistNESDtracks", "Number of ESD friend tracks", 1000, -0.5, 999.5);
 111   fHistNESDtracks->GetXaxis()->SetTitle("n. ESD friend tracks");
 112   fHistNESDtracks->GetYaxis()->SetTitle("dN/d(n. ESD friend tracks)");
 113   fHistNESDtracks->SetMarkerStyle(kFullCircle);
 114
 115   fHistNESDfriendtracks = new TH1F("fHistNESDfriendtracks", "Number of ESD tracks", 1000, -0.5, 999.5);
 116   fHistNESDfriendtracks->GetXaxis()->SetTitle("n. ESD tracks");
 117   fHistNESDfriendtracks->GetYaxis()->SetTitle("dN/d(n. ESD tracks)");
 118   fHistNESDfriendtracks->SetMarkerStyle(kFullCircle);
 119
 120   fListOut->Add(fHistPt);
 121   fListOut->Add(fHistQ);
 122   fListOut->Add(fHistNTPCCl);
 123   fListOut->Add(fHistNESDtracks);
 124   fListOut->Add(fHistNESDfriendtracks);
 125
 126   PostData(0, fListOut);
 127
 128   fCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
 129 }
 130
 131 //________________________________________________________________________
 132 void AliAnalysisTaskPt::Exec(Option_t *)
 133 {
 134   // Main loop
 135   // Called for each event
 136
 137   if (!fESD) {
 138     Printf("ERROR: fESD not available");
 139     return;
 140   }
 141   if (!fESDfriend) {
 142     Printf("ERROR: fESDfriend not available");
 143       if (fUseFriends){
 144         return;
 145       }
 146   }
 147
 148   Int_t nESDtracks = fESD->GetNumberOfTracks();
 149   Int_t nESDfriendtracks = 0;
 150   if (fUseFriends) nESDfriendtracks = fESDfriend->GetNumberOfTracks();
 151   Printf("There are %d tracks in this event", nESDtracks);
 152   Printf("... and there are %d friends in this event", nESDfriendtracks);
 153
 154   fHistNESDtracks->Fill(nESDtracks);
 155   fHistNESDfriendtracks->Fill(nESDfriendtracks);
 156
 157   // Track loop to fill a pT spectrum
 158   for (Int_t iTracks = 0; iTracks < nESDtracks; iTracks++) {
 159     Printf("Checking track %d: Note that with Flat, the GetTrack is not yet implemented!!!", iTracks);
 160     const AliVTrack* track = dynamic_cast<AliVTrack*>(fESD->GetTrack(iTracks));
 161     if (!track) {
 162       Printf("ERROR: Could not receive track %d", iTracks);
 163       continue;
 164     }
 165     Printf("track %d has pt = %f", iTracks, track->Pt());
 166     fHistPt->Fill(track->Pt());
 167     fHistNTPCCl->Fill(track->GetTPCNcls());
 168   } //track loop
 169
 170
 171   if (fUseFriends){
 172     // Friend Track loop
 173     for (Int_t iFriend = 0; iFriend < nESDfriendtracks; iFriend++) {
 174       //Printf("Getting friend %d", iFriend);
 175       const AliVfriendTrack* friendTrack = fESDfriend->GetTrack(iFriend);
 176       if (!friendTrack) {
 177         Printf("ERROR: Could not receive track %d", iFriend);
 178         continue;
 179       }
 180
 181       AliTPCseed seed;
 182       Int_t err = friendTrack->GetTPCseed( seed );
 183       if( err==0 ){
 184         Printf("Found TPC seed" );
 185         for (Int_t irow = 0; irow < 160; irow++){
 186           AliTPCclusterMI* cluMI = seed.GetClusterPointer(irow);
 187           if (cluMI){
 188             Printf("Found cluster at row %d", irow);
 189             Float_t q = cluMI->GetQ();
 190             Printf("Q = %f", q);
 191             fHistQ->Fill(q);
 192           }
 193           else {
 194             Printf("Row %d does not contain clusters", irow);
 195           }
 196         }
 197       }
 198       else {
 199         //Printf("Schade... seed is %p", seed);
 200       }
 201     }
 202   }
 203
 204   // Post output data.
 205   PostData(0, fListOut);
 206   fEv++;
 207 }
 208
 209 //________________________________________________________________________
 210 void AliAnalysisTaskPt::Terminate(Option_t *)
 211 {
 212   // Draw result to the screen
 213   // Called once at the end of the query
 214
 215   Printf("Terminate called: fESD = %p", fESD);
 216
 217   fListOut = dynamic_cast<TList*> (GetOutputData(0));
 218
 219   if (fListOut) {
 220     fHistPt = dynamic_cast<TH1F*>(fListOut->FindObject("fHistPt"));
 221     if (!fHistPt) {
 222       Printf("ERROR: fHistPt not available");
 223       return;
 224     }
 225
 226     TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
 227     c1->cd(1)->SetLogy();
 228     fHistPt->DrawCopy("E");
 229   }
 230   else {
 231     Printf("In Terminate: no TList found");
 232   }
 233
 234 }