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