PWG4/JetTasks/AliPWG4CosmicCandidates.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------------
  17 // Analysis task looking for cosmic candidates embedded in a p-p event
  18 // Task checks if particles are back-to-back in eta and phi
  19 //
  20 //
  21 // Author : Marta Verweij - UU - marta.verweij@cern.ch
  22 //-----------------------------------------------------------------------
  23
  24 #include "TVector3.h"
  25 #include <iostream>
  26 #include "TH1.h"
  27 #include "TH2.h"
  28 #include "TH3.h"
  29 #include "TList.h"
  30 #include "TChain.h"
  31
  32 #include "AliAnalysisTask.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37
  38 #include "AliESDtrack.h"
  39 #include "AliESDtrackCuts.h"
  40 #include "AliExternalTrackParam.h"
  41
  42 #include "AliLog.h"
  43
  44 #include "AliPWG4CosmicCandidates.h"
  45
  46 //using namespace std; //required for resolving the 'cout' symbol
  47 using namespace std;
  48
  49 ClassImp(AliPWG4CosmicCandidates)
  50
  51 //________________________________________________________________________
  52 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
  53 : AliAnalysisTaskSE(),
  54   fTrackCuts(0),
  55   fPtMin(5.),
  56   fMaxCosmicAngle(0.002),
  57   fNEventAll(0),
  58   fNEventSel(0),
  59   fPtSignedCosmicCandidates(0),
  60   fDeltaPtCosmicCandidates(0),
  61   fDeltaPhiSumEta(0),
  62   fDCAZCosmicCandidates(0),
  63   fDCARCosmicCandidates(0),
  64   fTheta(0),
  65   fThetaZoom(0),
  66   fThetaPt1Pt2(0),
  67   fThetaPt1Pt2Signed(0),
  68   fDeltaPhiSumEtaPt1(0),
  69   fDeltaPhiSumEtaPt2(0),
  70   fThetaDCAZ1DCAZ2(0),
  71   fRisol(0),
  72   fRisolTheta(0),
  73   fHistListCosmics(0)
  74 {
  75   //
  76   // Default constructor
  77   //
  78 }
  79
  80 //________________________________________________________________________
  81 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
  82   : AliAnalysisTaskSE(name),
  83     fTrackCuts(0),
  84     fPtMin(5.),
  85     fMaxCosmicAngle(0.002),
  86     fNEventAll(0),
  87     fNEventSel(0),
  88     fPtSignedCosmicCandidates(0),
  89     fDeltaPtCosmicCandidates(0),
  90     fDeltaPhiSumEta(0),
  91     fDCAZCosmicCandidates(0),
  92     fDCARCosmicCandidates(0),
  93     fTheta(0),
  94     fThetaZoom(0),
  95     fThetaPt1Pt2(0),
  96     fThetaPt1Pt2Signed(0),
  97     fDeltaPhiSumEtaPt1(0),
  98     fDeltaPhiSumEtaPt2(0),
  99     fThetaDCAZ1DCAZ2(0),
 100     fRisol(0),
 101     fRisolTheta(0),
 102     fHistListCosmics(0)
 103 {
 104   // Constructor. Initialization of Inputs and Outputs
 105
 106   // Define input and output slots here
 107   // Input slot #0 works with a TChain
 108   //  DefineInput(0, TChain::Class());
 109   // Output slot #0 id reserved by the base class for AOD
 110   // Output slot #1 writes into a TList
 111   DefineOutput(1, TList::Class());
 112   // Output slot #2 writes into a AliESDtrackCuts
 113   DefineOutput(2, AliESDtrackCuts::Class());
 114 }
 115 //________________________________________________________________________
 116 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
 117   : AliAnalysisTaskSE(res),
 118     fTrackCuts(0),
 119     fPtMin(5.),
 120     fMaxCosmicAngle(0.002),
 121     fNEventAll(0),
 122     fNEventSel(0),
 123     fPtSignedCosmicCandidates(0),
 124     fDeltaPtCosmicCandidates(0),
 125     fDeltaPhiSumEta(0),
 126     fDCAZCosmicCandidates(0),
 127     fDCARCosmicCandidates(0),
 128     fTheta(0),
 129     fThetaZoom(0),
 130     fThetaPt1Pt2(0),
 131     fThetaPt1Pt2Signed(0),
 132     fDeltaPhiSumEtaPt1(0),
 133     fDeltaPhiSumEtaPt2(0),
 134     fThetaDCAZ1DCAZ2(0),
 135     fRisol(0),
 136     fRisolTheta(0),
 137     fHistListCosmics(0)
 138 {
 139   //
 140   // Dummy copy constructor
 141   //
 142 }
 143
 144 //________________________________________________________________________
 145 AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
 146 {
 147   //
 148   // Dummy assignment operator
 149   //
 150   return *this;
 151 }
 152
 153 //________________________________________________________________________
 154 void AliPWG4CosmicCandidates::LocalInit()
 155 {
 156   //
 157   // Only called once at beginning
 158   //
 159   PostData(2,fTrackCuts);
 160 }
 161
 162 //________________________________________________________________________
 163 void AliPWG4CosmicCandidates::UserCreateOutputObjects()
 164 {
 165   //
 166   // Create output objects
 167   // Called once
 168   AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n"));
 169
 170   Bool_t oldStatus = TH1::AddDirectoryStatus();
 171   TH1::AddDirectory(kFALSE);
 172
 173   OpenFile(1);
 174   fHistListCosmics = new TList();
 175
 176   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 177   fHistListCosmics->Add(fNEventAll);
 178   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 179   fHistListCosmics->Add(fNEventSel);
 180
 181   Float_t fgkPtMin=0.;
 182   Float_t fgkPtMax=100.;
 183   Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
 184   Int_t fgkNPtBins2D= (int)((fgkPtMax-fgkPtMin)/4.);
 185
 186   Int_t fgkNPhiBins=18;
 187   Float_t kMinPhi = -0.5*TMath::Pi();
 188   Float_t kMaxPhi = 3./2.*TMath::Pi();
 189
 190   Int_t fgkNThetaBins=fgkNPhiBins*8;
 191   Float_t kMinTheta = -0.5*TMath::Pi();
 192   Float_t kMaxTheta = 3./2.*TMath::Pi();
 193
 194   Int_t fgkNDCARBins=40;
 195   Float_t fgkDCARMin = -0.2;
 196   Float_t fgkDCARMax = 0.2;
 197   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
 198   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
 199
 200   Int_t fgkNDCAZBins=40;
 201   Float_t fgkDCAZMin = -2.;
 202   Float_t fgkDCAZMax = 2.;
 203   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
 204   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
 205
 206   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*fgkNPtBins, -1.*fgkPtMax, fgkPtMax);
 207   fHistListCosmics->Add(fPtSignedCosmicCandidates);
 208
 209   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
 210   fHistListCosmics->Add(fDeltaPtCosmicCandidates);
 211
 212   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
 213   fHistListCosmics->Add(fDeltaPhiSumEta);
 214
 215   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
 216   fHistListCosmics->Add(fDCAZCosmicCandidates);
 217
 218   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
 219   fHistListCosmics->Add(fDCARCosmicCandidates);
 220
 221   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
 222   fHistListCosmics->Add(fTheta);
 223
 224   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
 225   fHistListCosmics->Add(fThetaZoom);
 226
 227   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,fgkPtMin,fgkPtMax,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 228   fHistListCosmics->Add(fThetaPt1Pt2);
 229
 230   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax);
 231   fHistListCosmics->Add(fThetaPt1Pt2Signed);
 232
 233   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 234   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
 235
 236   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 237   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
 238
 239   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
 240   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
 241
 242   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
 243   fHistListCosmics->Add(fRisol);
 244
 245   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
 246   fHistListCosmics->Add(fRisolTheta);
 247
 248   TH1::AddDirectory(oldStatus);
 249
 250   if(binsDCAR) delete [] binsDCAR;
 251   if(binsDCAZ) delete [] binsDCAZ;
 252
 253 }
 254
 255 //________________________________________________________________________
 256 void AliPWG4CosmicCandidates::UserExec(Option_t *)
 257 {
 258   // Main loop
 259   // Called for each event
 260   //
 261
 262   // All events without selection
 263   fNEventAll->Fill(0.);
 264
 265   if (!fInputEvent) {
 266     AliDebug(2,Form("ERROR: fESD not available"));
 267     cout << "ERROR: fESD not available" << endl;
 268     PostData(1, fHistListCosmics);
 269     return;
 270   }
 271
 272   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
 273   // Need vertex cut
 274   TString vtxName(vtx->GetName());
 275   if(!vtx || vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 276     // Post output data
 277     PostData(1, fHistListCosmics);
 278     return;
 279   }
 280
 281   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 282   double primVtx[3];
 283   vtx->GetXYZ(primVtx);
 284   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 285     // Post output data
 286     PostData(1, fHistListCosmics);
 287     return;
 288   }
 289   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){
 290     // Post output data
 291     PostData(1, fHistListCosmics);
 292     return;
 293   }
 294   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 295
 296   if(!fTrackCuts) {
 297    // Post output data
 298     PostData(1, fHistListCosmics);
 299     return;
 300   }
 301
 302   fNEventSel->Fill(0.);
 303
 304   Float_t dcaR[2] = {0.,0.};
 305   Float_t dcaZ[2] = {0.,0.};
 306
 307   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 308
 309     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
 310     if (!track1)  continue;
 311     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
 312     if(track1->Pt()<fPtMin) continue;
 313     //Start 2nd track loop to look for correlations
 314     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 315       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
 316       if(!track2) continue;
 317       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
 318
 319       //Check if back-to-back
 320       Double_t mom1[3],mom2[3];
 321       track1->GetPxPyPz(mom1);
 322       track2->GetPxPyPz(mom2);
 323       //     Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
 324       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 325       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 326       //Double_t theta = momv1.Angle(momv2);
 327       Double_t theta = momv1.Phi()-momv2.Phi();
 328       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 329
 330       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
 331       Float_t deltaPhi = track1->Phi()-track2->Phi();
 332       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 333       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
 334
 335       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
 336       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
 337
 338       if(track2->Pt()<0.5) continue;
 339       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
 340       fRisol->Fill(rIsol); //Fill R histogram
 341       if(track2->Pt()<fPtMin) continue;
 342
 343         fTheta->Fill(theta);
 344         fThetaZoom->Fill(theta);
 345         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
 346         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
 347         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
 348         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
 349         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
 350         fRisolTheta->Fill(rIsol,theta);
 351         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 352           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
 353           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
 354           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
 355         }
 356
 357     } // track2 loop
 358
 359   } // track1 loop
 360
 361   PostData(1,fHistListCosmics);
 362 }
 363
 364 //________________________________________________________________________
 365 void AliPWG4CosmicCandidates::Terminate(Option_t *)
 366 {
 367   //
 368   // Called once at the end of the query
 369   //
 370
 371 }