PWGCF/Correlations/DPhi/AliPhiCorrelationsQATask.cxx

   1 /* $Id: AliPhiCorrelationsQATask.cxx 47416 2011-02-15 13:55:09Z jgrosseo $ */
   2
   3 #include "AliPhiCorrelationsQATask.h"
   4
   5 #include <TH2F.h>
   6 #include <TParticle.h>
   7 #include <TFile.h>
   8
   9 #include <AliLog.h>
  10 #include <AliESDVertex.h>
  11 #include <AliESDEvent.h>
  12 #include <AliMCEvent.h>
  13 #include <AliStack.h>
  14 #include <AliHeader.h>
  15 #include <AliAnalysisManager.h>
  16 #include <AliCentrality.h>
  17
  18 #include "AliESDtrackCuts.h"
  19
  20 ClassImp(AliPhiCorrelationsQATask)
  21
  22 AliPhiCorrelationsQATask::AliPhiCorrelationsQATask(const char* opt) :
  23   AliAnalysisTaskSE("AliPhiCorrelationsQATask"),
  24   fOutput(0),
  25   fOption(opt),
  26   fTPCOnly(0),
  27   fEsdTrackCuts(0),
  28   fEsdTrackCuts2(0),
  29   fCheckITS(0),
  30   fGlobalTracks(0),
  31   fCentralityCorrelation(0),
  32   fDCAPrimaries(0),
  33   fDCASecondaries(0),
  34   fUseUncheckedCentrality(kFALSE)
  35 {
  36   //
  37   // Constructor. Initialization of pointers
  38   //
  39
  40   // Define input and output slots here
  41   DefineOutput(1, TList::Class());
  42
  43   AliLog::SetClassDebugLevel("AliPhiCorrelationsQATask", AliLog::kWarning);
  44 }
  45
  46 AliPhiCorrelationsQATask::~AliPhiCorrelationsQATask()
  47 {
  48   //
  49   // Destructor
  50   //
  51
  52   // histograms are in the output list and deleted when the output
  53   // list is deleted by the TSelector dtor
  54
  55   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  56     delete fOutput;
  57     fOutput = 0;
  58   }
  59 }
  60
  61 void AliPhiCorrelationsQATask::UserCreateOutputObjects()
  62 {
  63   // create result objects and add to output list
  64
  65   Printf("AliPhiCorrelationsQATask::CreateOutputObjects");
  66   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
  67
  68   fOutput = new TList;
  69   fOutput->SetOwner();
  70
  71   fCentralityCorrelation = new TH2F("fCentralityCorrelation", ";v0 centr;spd centr", 100, 0, 100.001, 100, 0, 100.001);
  72   fOutput->Add(fCentralityCorrelation);
  73
  74    if (fEsdTrackCuts)
  75      fOutput->Add(fEsdTrackCuts);
  76    if (fEsdTrackCuts2)
  77      fOutput->Add(fEsdTrackCuts2);
  78
  79 //   {
  80 //     fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  81 //     fEsdTrackCuts->SetMinNClustersTPC(70);
  82 //     // disable DCA cut
  83 //     fEsdTrackCuts->SetMaxDCAToVertexZ();
  84 //     fEsdTrackCuts->SetMaxDCAToVertexXY();
  85 //     fEsdTrackCuts->SetDCAToVertex2D(kFALSE);
  86 //     fEsdTrackCuts->SetName("cuts_quality_only");
  87 //     fEsdTrackCuts->DefineHistograms();
  88 //     fOutput->Add(fEsdTrackCuts);
  89 //   }
  90 //
  91 //   fEsdTrackCuts2 = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  92 //   fEsdTrackCuts2->SetMinNClustersTPC(70);
  93 //   fEsdTrackCuts2->SetName("cuts_quality_dca");
  94 //   fEsdTrackCuts2->DefineHistograms();
  95 //   fOutput->Add(fEsdTrackCuts2);
  96
  97   fGlobalTracks = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
  98   fGlobalTracks->SetName("global_cuts");
  99   fGlobalTracks->DefineHistograms();
 100   fOutput->Add(fGlobalTracks);
 101
 102   fCheckITS = new AliESDtrackCuts;
 103   fCheckITS->SetRequireITSRefit(kTRUE);
 104   fCheckITS->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 105                                       AliESDtrackCuts::kAny);
 106   fCheckITS->SetDCAToVertex2D(kFALSE);
 107   fCheckITS->SetRequireSigmaToVertex(kFALSE);
 108   fCheckITS->SetName("check_its");
 109   fCheckITS->DefineHistograms();
 110   fOutput->Add(fCheckITS);
 111
 112   Int_t nBins[] =   { 20, 10, 200, 100 };
 113   Double_t xMin[] = { 0, 0, -5, -5 };
 114   Double_t xMax[] = { 20, 100, 5, 5 };
 115
 116   fDCAPrimaries = new THnF("fDCAPrimaries", ";pT;centrality;dca_xy;dca_z", 4, nBins, xMin, xMax);
 117   fDCASecondaries = (THnF*) fDCAPrimaries->Clone("fDCASecondaries");
 118
 119   fOutput->Add(fDCAPrimaries);
 120   fOutput->Add(fDCASecondaries);
 121
 122   PostData(1, fOutput);
 123 }
 124
 125 void AliPhiCorrelationsQATask::UserExec(Option_t*)
 126 {
 127   // process the event
 128
 129   //Printf("AliPhiCorrelationsQATask::UserExec");
 130
 131   AliESDEvent* esd = (AliESDEvent*) fInputEvent;
 132   if (!esd)
 133     AliFatal("No input event");
 134
 135   // vertex cut
 136   const AliESDVertex* vtxESD = esd->GetPrimaryVertex();
 137   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 138   if (!vtxESD)
 139     return;
 140   if (!vtxSPD)
 141     return;
 142
 143   Double_t vtx[3];
 144   vtxESD->GetXYZ(vtx);
 145   if (TMath::Abs(vtx[2]) > 10)
 146     return;
 147
 148   AliCentrality *centralityObj = esd->GetCentrality();
 149   Float_t v0Centrality = -1;
 150   if (centralityObj)
 151   {
 152     Float_t spdCentrality = -1;
 153
 154     if (fUseUncheckedCentrality)
 155     {
 156       v0Centrality = centralityObj->GetCentralityPercentileUnchecked("V0M");
 157       spdCentrality = centralityObj->GetCentralityPercentileUnchecked("CL1");
 158     }
 159     else
 160     {
 161       v0Centrality = centralityObj->GetCentralityPercentile("V0M");
 162       spdCentrality = centralityObj->GetCentralityPercentile("CL1");
 163     }
 164
 165     fCentralityCorrelation->Fill(v0Centrality, spdCentrality);
 166   }
 167
 168   AliStack* stack = 0;
 169   if (fMCEvent)
 170   {
 171     stack = fMCEvent->Stack();
 172     if (!stack)
 173       AliFatal("Stack is 0");
 174   }
 175
 176   Int_t count = 0;
 177
 178   for (Int_t i=0; i<esd->GetNumberOfTracks(); i++)
 179   {
 180     AliESDtrack* esdTrack = esd->GetTrack(i);
 181
 182     if (0)
 183     {
 184       if (!fGlobalTracks->AcceptTrack(esdTrack))
 185         continue;
 186     }
 187
 188     Bool_t accept = kFALSE;
 189     if (fEsdTrackCuts->AcceptTrack(esdTrack))
 190       accept = kTRUE;
 191
 192     if (fEsdTrackCuts2 && fEsdTrackCuts2->AcceptTrack(esdTrack))
 193       accept = kTRUE;
 194
 195     if (!accept)
 196       continue;
 197
 198     if (fTPCOnly)
 199     {
 200       AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esd, i);
 201       if(!track) continue;
 202
 203       if(track->Pt()>0.){
 204         // only constrain tracks above threshold
 205         AliExternalTrackParam exParam;
 206         // take the B-feild from the ESD, no 3D fieldMap available at this point
 207         Bool_t relate = kFALSE;
 208         relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 209         if(!relate)
 210         {
 211 //                 Printf("relating failed");
 212           delete track;
 213           continue;
 214         }
 215         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 216       }
 217       esdTrack = track;
 218     }
 219
 220     count++;
 221
 222     //if (v0Centrality < 0)
 223     //  continue;
 224
 225     Float_t b[2];
 226     Float_t bCov[3];
 227     esdTrack->GetImpactParameters(b, bCov);
 228
 229     Bool_t primary = kTRUE;
 230     if (stack)
 231       primary = stack->IsPhysicalPrimary(TMath::Abs(esdTrack->GetLabel()));
 232
 233     Double_t fillArray[4] = { esdTrack->Pt(), v0Centrality, b[0], b[1] };
 234
 235     if (primary)
 236       fDCAPrimaries->Fill(fillArray);
 237     else
 238       fDCASecondaries->Fill(fillArray);
 239
 240     if (fTPCOnly)
 241       delete esdTrack;
 242
 243     // fill histograms of second object
 244 //     if (!fEsdTrackCuts2->AcceptTrack(esdTrack))
 245 //       continue;
 246 //
 247 //     if (fCheckITS->AcceptTrack(esdTrack))
 248 //     {
 249 //       if (!fGlobalTracks->AcceptTrack(esdTrack))
 250 //         continue;
 251 //     }
 252
 253     /*
 254     // create a tpc only track
 255     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(esd,esdTrack->GetID());
 256     if (!track)
 257       continue;
 258
 259     if (track->Pt()>0.)
 260     {
 261       // only constrain tracks above threshold
 262       AliExternalTrackParam exParam;
 263       // take the B-feild from the ESD, no 3D fieldMap available at this point
 264       Bool_t relate = false;
 265       relate = track->RelateToVertex(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 266       if(!relate){
 267         delete track;
 268         continue;
 269       }
 270       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 271     }
 272
 273     Float_t b2[2];
 274     track->GetImpactParameters(b2,bCov);
 275
 276     //Printf("%.2f %.2f %.2f %.2f", b[0], b2[0], b[1], b2[1]);
 277
 278     delete track;
 279     */
 280   }
 281
 282   //Printf("%.2f: %d %d tracks (out of %d)", v0Centrality, centralityObj->GetQuality(), count, esd->GetNumberOfTracks());
 283 }
 284
 285 void AliPhiCorrelationsQATask::Terminate(Option_t *)
 286 {
 287   // terminate
 288
 289   Printf("In AliPhiCorrelationsQATask::Terminate...");
 290
 291   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 292   if (!fOutput)
 293   {
 294     Printf("ERROR: fOutput not available");
 295     return;
 296   }
 297
 298   fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("cuts_quality_only"));
 299   if (!fEsdTrackCuts)
 300   {
 301     Printf("ERROR: fEsdTrackCuts not available");
 302     return;
 303   }
 304
 305   fEsdTrackCuts2 = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("cuts_quality_dca"));
 306   if (!fEsdTrackCuts2)
 307   {
 308     Printf("ERROR: fEsdTrackCuts2 not available");
 309     return;
 310   }
 311
 312   fGlobalTracks = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("global_cuts"));
 313   if (!fGlobalTracks)
 314   {
 315     Printf("ERROR: fGlobalTracks not available");
 316     return;
 317   }
 318
 319   TFile* file = TFile::Open("track_cuts.root", "RECREATE");
 320   fEsdTrackCuts->SaveHistograms();
 321   fEsdTrackCuts2->SaveHistograms();
 322   fGlobalTracks->SaveHistograms();
 323   file->Close();
 324 }