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