[u/mrichter/AliRoot.git] / PWG1 / cosmic / AliAnalysisTaskCosmic.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

// Analysis Task for the Quality Assurance of Cosmic Data
// Two track segments in the are matched in angle and charged. T
// The quality of the matching in is checked by comparing 
// the tarnsverse momenta and starting points of the track segments
//
// Author
// Andreas Morsch
// andreas.morsch@cern.ch

#include "TChain.h"
#include "TTree.h"

#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TList.h"

#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliAnalysisTaskCosmic.h"

ClassImp(AliAnalysisTaskCosmic)

//________________________________________________________________________
AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const char *name) 
    : AliAnalysisTaskSE(name),
      fHists(0),
      fhDZvsZ(0),
      fhDZvsPhi(0),
      fhCh1Ch2(0),
      fhPh1Ph2(0),
      fhCl1Cl2G(0),
      fhCl1Cl2B(0)
{
  //
  // Constructor
  //
    
    for (Int_t i = 0; i < 6; i++) {
	fhPt[i]      = 0;
	fhTheta[i]   = 0;
	fhPhi[i]     = 0;
	fhDPhi[i]    = 0;
	fhDTheta[i]  = 0;
	fhDZ[i]      = 0;
	fhDX[i]      = 0;
	fhDY[i]      = 0;
	fhDPt[i]     = 0;
	fhD1ovPt[i]  = 0;
	fpDPt[i]     = 0;
	fhDPtovPt[i] = 0;
	fpDPtS[i]    = 0;
    }
    
  DefineOutput(1,  TList::Class());
}

AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const AliAnalysisTaskCosmic& task)
    : AliAnalysisTaskSE(task),
      fHists(0),
      fhDZvsZ(0),
      fhDZvsPhi(0),
      fhCh1Ch2(0),
      fhPh1Ph2(0),
      fhCl1Cl2G(0),
      fhCl1Cl2B(0)
{
  // Copy constructor
    fHists       = new TList();
    fhDZvsZ      = (TH2F*)    ((task.fhDZvsZ)   ->Clone()); 
    fhDZvsPhi    = (TH2F*)    ((task.fhDZvsPhi) ->Clone());
    fhCh1Ch2     = (TH2F*)    ((task.fhCh1Ch2)  ->Clone());
    fhPh1Ph2     = (TH2F*)    ((task.fhPh1Ph2)  ->Clone());
    fhCl1Cl2G    = (TH2F*)    ((task.fhCl1Cl2G) ->Clone());
    fhCl1Cl2B    = (TH2F*)    ((task.fhCl1Cl2B) ->Clone());
    for (Int_t i = 0; i < 6; i++) {
	fhPt[i]       =    (TH1F*) ((task.fhPt[i])        ->Clone());
	fhTheta[i]    =    (TH1F*) ((task.fhTheta[i])     ->Clone());
	fhPhi[i]      =    (TH1F*) ((task.fhPhi[i])       ->Clone());
	fhDPhi[i]     =    (TH1F*) ((task.fhDPhi[i])      ->Clone());
	fhDTheta[i]   =    (TH1F*) ((task.fhDTheta[i])    ->Clone());
	fhDZ[i]       =    (TH1F*) ((task.fhDZ[i])        ->Clone());
	fhDX[i]       =    (TH1F*) ((task.fhDX[i])        ->Clone());
	fhDY[i]       =    (TH1F*) ((task.fhDY[i])        ->Clone());
	fhDPt[i]      =    (TH1F*) ((task.fhDPt[i])       ->Clone());
	fhD1ovPt[i]   =    (TH1F*) ((task.fhD1ovPt[i])    ->Clone());
	fhDPtovPt[i]  =    (TH1F*) ((task.fhDPtovPt[i])   ->Clone());

	fpDPt[i]      =    (TProfile*) ((task.fpDPt[i])       ->Clone());
	fpDPtS[i]     =    (TProfile*) ((task.fpDPtS[i])      ->Clone());
    }
}

//______________________________________________________________________________
AliAnalysisTaskCosmic& AliAnalysisTaskCosmic::operator=(const AliAnalysisTaskCosmic& task)
{
  // Assignment operator
  if(this!=&task) {

    AliAnalysisTaskSE::operator=(task);

    fHists       = new TList();
    fhDZvsZ      = (TH2F*)    ((task.fhDZvsZ)   ->Clone()); 
    fhDZvsPhi    = (TH2F*)    ((task.fhDZvsPhi) ->Clone());
    fhCh1Ch2     = (TH2F*)    ((task.fhCh1Ch2)  ->Clone());
    fhPh1Ph2     = (TH2F*)    ((task.fhPh1Ph2)  ->Clone());
    fhCl1Cl2G    = (TH2F*)    ((task.fhCl1Cl2G) ->Clone());
    fhCl1Cl2B    = (TH2F*)    ((task.fhCl1Cl2B) ->Clone());
    for (Int_t i = 0; i < 6; i++) {
	fhPt[i]       =    (TH1F*) ((task.fhPt[i])        ->Clone());
	fhTheta[i]    =    (TH1F*) ((task.fhTheta[i])     ->Clone());
	fhPhi[i]      =    (TH1F*) ((task.fhPhi[i])       ->Clone());
	fhDPhi[i]     =    (TH1F*) ((task.fhDPhi[i])      ->Clone());
	fhDTheta[i]   =    (TH1F*) ((task.fhDTheta[i])    ->Clone());
	fhDZ[i]       =    (TH1F*) ((task.fhDZ[i])        ->Clone());
	fhDX[i]       =    (TH1F*) ((task.fhDX[i])        ->Clone());
	fhDY[i]       =    (TH1F*) ((task.fhDY[i])        ->Clone());
	fhDPt[i]      =    (TH1F*) ((task.fhDPt[i])       ->Clone());
	fhD1ovPt[i]   =    (TH1F*) ((task.fhD1ovPt[i])    ->Clone());
	fhDPtovPt[i]  =    (TH1F*) ((task.fhDPtovPt[i])   ->Clone());

	fpDPt[i]      =    (TProfile*) ((task.fpDPt[i])       ->Clone());
	fpDPtS[i]     =    (TProfile*) ((task.fpDPtS[i])      ->Clone());
    }
  }
  return *this;
}


//________________________________________________________________________
void AliAnalysisTaskCosmic::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
    TString ext[6] = {"PC", "NC", "PZ", "NZ", "_Good", "_Bad"};
    
    char name[12];
    
    fHists = new TList();

    for (Int_t i = 0; i < 6; i++) {
	// Pt
	snprintf(name, 12, "fhPt%2s", ext[i].Data());
	fhPt[i]   = new TH1F(name, " pT distribution",        800, 0., 200.);
	fhPt[i]->SetXTitle("p_{T} [Gev]");
	// Phi
	snprintf(name, 12, "fhPhi%2s", ext[i].Data());
	fhPhi[i]   = new TH1F(name,   "Phi distribution",        62,  0., 2. * TMath::Pi());
	fhPhi[i]->SetXTitle("#phi [rad]");
	// Theta
	snprintf(name, 12, "fhTheta%2s", ext[i].Data());
	fhTheta[i]   = new TH1F(name,   "Theta distribution",    62,  0., TMath::Pi());    
	fhTheta[i]->SetXTitle("#theta [rad]");
	// Delta Phi
	snprintf(name, 12, "fhDPhi%2s", ext[i].Data());
	fhDPhi[i]   = new TH1F(name,   "DeltaPhi distribution",        320, -0.4, 0.4);
	fhDPhi[i]->SetXTitle("#Delta#phi [rad]");
	// Delta Theta
	snprintf(name, 12, "fhDTheta%2s", ext[i].Data());
	fhDTheta[i]   = new TH1F(name,   "DeltaTheta distribution",    320, -0.4, 0.4);    
	fhDTheta[i]->SetXTitle("#Delta#theta [rad]");
	// Delta Z
	snprintf(name, 12, "fhDZ%2s", ext[i].Data());
	fhDZ[i]   = new TH1F(name,   "DeltaZ distribution",        200, -10., 10.);
	fhDZ[i]->SetXTitle("#DeltaZ [cm]");
	// Delta X
	snprintf(name, 12, "fhDX%2s", ext[i].Data());
	fhDX[i]   = new TH1F(name,   "DeltaX distribution",        200, -10., 10.);
	fhDX[i]->SetXTitle("#DeltaX [cm]");
	// Delta Y
	snprintf(name, 12, "fhDY%2s", ext[i].Data());
	fhDY[i]   = new TH1F(name,   "DeltaY distribution",        200, -10, 10.);
	fhDY[i]->SetXTitle("#DeltaY [cm]");
	// Delta Pt
	snprintf(name, 12, "fhDPt%2s", ext[i].Data());
	fhDPt[i]   = new TH1F(name,   "DeltaPt distribution",        200, -20., 20.);
	fhDPt[i]->SetXTitle("#Delta p_{T}  [GeV]");

	// Delta 1/Pt
	snprintf(name, 12, "fhD1ovPt%2s", ext[i].Data());
	fhD1ovPt[i]   = new TH1F(name,   "Delta 1/Pt distribution",        200, -1., 1.);
	fhD1ovPt[i]->SetXTitle("#Delta 1/Pt");

	// Delta Pt over Pt
	snprintf(name, 12, "fhDPtovPt%2s", ext[i].Data());
	fhDPtovPt[i]   = new TH1F(name,   "DeltaPt/Pt distribution",        200, -2., 2.);
	fhDPtovPt[i]->SetXTitle("#DeltaPt/Pt");


	// Delta Pt/ Pt vs Pt
	snprintf(name, 12, "fpDPt%2s", ext[i].Data());
	fpDPt[i] = new TProfile(name, "#Delta Pt / Pt", 20, 0., 20., -1, 1., "S");
	fpDPt[i]->SetXTitle("p_{T} [GeV]");
	fpDPt[i]->SetYTitle("#Delta 1/p_{T} [GeV^{-1}]");
	// Delta Pt error
	snprintf(name, 12, "fpDPtS%2s", ext[i].Data());
	fpDPtS[i] = new TProfile(name, "#Delta Pt / <sigma>", 20, 0., 20., 0., 10.);
	fpDPtS[i]->SetXTitle("p_{T}");
	fpDPtS[i]->SetYTitle("#Delta p_{T} / <#sigma_{p_{T}}>");


	fHists->Add(fhPt[i]);
	fHists->Add(fhPhi[i]);
	fHists->Add(fhTheta[i]);
	fHists->Add(fhDPhi[i]);
	fHists->Add(fhDPt[i]);
	fHists->Add(fhD1ovPt[i]);
	fHists->Add(fhDPtovPt[i]);
	fHists->Add(fhDTheta[i]);
	fHists->Add(fhDZ[i]);
	fHists->Add(fhDX[i]);
	fHists->Add(fhDY[i]);
	fHists->Add(fpDPt[i]);
	fHists->Add(fpDPtS[i]);
    }

    fhDZvsZ      = new TH2F("fhDZvsZ",    "dz vs z", 500, -250., 250., 100, -10., 10.);
    fhDZvsZ->SetXTitle("z_{in} * sign(z_{in}) * sign(z_{out}) [cm]");
    fhDZvsZ->SetYTitle("#Delta z [cm]");
    
    
    fhDZvsPhi    = new TH2F("fhDZvsPhi", "dz vs phi", 36, 0., TMath::Pi(),  50, -2., 2.);

    fhCh1Ch2   = new TH2F("fCh1Ch2",   "ch1 vs ch2", 8, -2., 2., 8, -2., 2.);
    fhPh1Ph2   = new TH2F("fPh1Ph2",   "ph1 vs ph2", 128, 0., 2. * TMath::Pi(), 128, 0., 2. * TMath::Pi());
    fhCl1Cl2G  = new TH2F("fCl1Cl2G",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);
    fhCl1Cl2B  = new TH2F("fCl1Cl2B",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);    

    fHists->Add(fhDZvsZ);
    fHists->Add(fhDZvsPhi);
    fHists->Add(fhCh1Ch2);
    fHists->Add(fhPh1Ph2);
    fHists->Add(fhCl1Cl2G);
    fHists->Add(fhCl1Cl2B);
    fHists->SetOwner();

}

//________________________________________________________________________
void AliAnalysisTaskCosmic::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event

  if (!fInputEvent) {
    Printf("ERROR: fESD not available");
    return;
  }
  
  AliESDEvent* esdE = (AliESDEvent*)fInputEvent;
  
  if (esdE->GetNumberOfTracks() != 2) {
      PostData(1, fHists);
      return;
  }
  
  
  for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
      AliESDtrack* track = esdE->GetTrack(iTracks);


//      const AliExternalTrackParam * track = trackG->GetTPCInnerParam();
//      if (!track) continue;
      // Pt cut
      Float_t pt     = track->Pt();
      Short_t charge = track->Charge();
      Float_t phi    = track->Phi();
      if (phi > 0. && phi  < TMath::Pi()) charge*=(-1);


      // TPC track
      UInt_t status = track->GetStatus();
      if ((status&AliESDtrack::kTPCrefit) ==0) continue;

      Int_t nClustersTPC = track->GetTPCclusters(0);
      if (nClustersTPC < 50) continue;
      
      // 

      Float_t z      = track->GetZ();
      Float_t x      = track->Xv();
      Float_t y      = track->Yv();
      Float_t theta  = track->Theta();


      const AliExternalTrackParam * trackOut = track->GetOuterParam();
      Float_t zOut = 0.;
      if (trackOut)zOut = trackOut->Zv();
      

      if (charge > 0) {
	  fhPt[kPosC]   ->Fill(pt);
      } else {
	  fhPt[kNegC]   ->Fill(pt);
      }

      //      if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (pt < 1. || pt > 50.)) continue;

      
      //
      if (charge > 0) {
	  fhPhi[kPosC]  ->Fill(phi);
	  fhTheta[kPosC]->Fill(theta);
      } else {
	  fhPhi[kNegC]  ->Fill(phi);
	  fhTheta[kNegC]->Fill(theta);
      }

      
      if (z > 0) {
	  fhPt[kPosZ]   ->Fill(pt);
	  fhPhi[kPosZ]  ->Fill(phi);
	  fhTheta[kPosZ]->Fill(theta);
      } else {
	  fhPt[kNegZ]   ->Fill(pt);
	  fhPhi[kNegZ]  ->Fill(phi);
	  fhTheta[kNegZ]->Fill(theta);
      }


      // Tracks coming from above
      if (phi > 0. &&  phi < TMath::Pi()) {
	  
//	  printf("Track#1 %5d %5d %13.3f %13.3f \n", Int_t(Entry()), iTracks, phi, zOut);
	  
	  Float_t dphiMin = 999.;
	  Float_t rMin    = 999.;
	  Int_t   jMin    = -1;
	  // Search for a matching track (in dphi-dtheta)
	  for (Int_t jTracks = 0; jTracks < esdE->GetNumberOfTracks(); jTracks++) {
	      if (jTracks == iTracks) continue;

	      AliESDtrack* track2 = esdE->GetTrack(jTracks);

	      UInt_t status2 = track2->GetStatus();
	      if ((status2&AliESDtrack::kTPCrefit) ==0) continue;

	      Int_t nClustersTPC2 = track2->GetTPCclusters(0);
	      if (nClustersTPC2 < 50) continue;
	      //	      if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (track2->Pt() < 1. || track2->Pt() > 50.)) continue;
//	      const AliExternalTrackParam * track2 = trackG2->GetTPCInnerParam();
//	      if (!track2) continue;
	      
	      Float_t phi2   = track2->Phi() - TMath::Pi();
	      Float_t theta2 = TMath::Pi()   - track2->Theta();
	      Float_t dphi   = phi2   - phi;
	      Float_t dtheta = theta2 - theta;
	      
	      if (dphi >  TMath::Pi()) dphi -= 2. * TMath::Pi();
	      if (dphi < -TMath::Pi()) dphi += 2. * TMath::Pi();

	      Float_t dR = TMath::Sqrt(dphi * dphi + dtheta * dtheta);

	      if (dR < rMin) {
		  rMin    = dR;
		  dphiMin = dphi;
		  jMin = jTracks;
	      }

	  } // tracks 2
	  if (jMin != -1) {
	      // we found a matching track candidate ...
	      AliESDtrack* track2 = esdE->GetTrack(jMin);
	      const AliExternalTrackParam * trackOut2 = track2->GetOuterParam();
	      Float_t zOut2 = 0.;
	      if (trackOut2) zOut2 = trackOut2->Zv();


	      Float_t theta2 = - track2->Theta() + TMath::Pi();
	      Float_t z2 = track2->GetZ();
	      Float_t x2 = track2->Xv();
	      Float_t y2 = track2->Yv();
	      Short_t charge2 = track2->Charge();
	      Float_t dz  = z2 - z;
	      Float_t dx  = x2 - x;
	      Float_t dy  = y2 - y;
	      Float_t pt1 = track->Pt();
	      Float_t pt2 = track2->Pt();
	      Float_t dpt = pt2 - pt1;
	      Float_t d1pt = 1./pt2 - 1./pt1;
	      
	      Float_t ptm = 0.5 * (pt1 + pt2);

	      Int_t nClustersTPC2 = track2->GetTPCclusters(0);

	      if (charge > 0.) {
		  fhDPhi[kPosC]  ->Fill(dphiMin);
		  fhDTheta[kPosC]->Fill(theta2 - theta);
	      } else {
		  fhDPhi[kNegC]  ->Fill(dphiMin);
		  fhDTheta[kNegC]->Fill(theta2 - theta);
	      }

	      if (z > 0.) {
		  fhDPhi[kPosZ]  ->Fill(dphiMin);
		  fhDTheta[kPosZ]->Fill(theta2 - theta);
	      } else {
		  fhDPhi[kNegZ]  ->Fill(dphiMin);
		  fhDTheta[kNegZ]->Fill(theta2 - theta);
	      }

	      if (TMath::Abs(dpt)/ptm > 0.5) {
		  fhDPhi[kBad]  ->Fill(dphiMin);
		  fhDTheta[kBad]->Fill(theta2 - theta);
	      } else {
		  fhDPhi[kGood]  ->Fill(dphiMin);
		  fhDTheta[kGood]->Fill(theta2 - theta);
	      }
	      // Good matches ...	      
//	      if (TMath::Abs(rMin < 0.04) && (charge == charge2) && TMath::Abs(dz) < 0.5 &&  TMath::Abs(dy) && TMath::Abs(dx) < 0.5 ) 
//	      if (TMath::Abs(rMin < 0.04) && (charge == charge2) && (zOut * zOut2) < 0.) 
	      if (TMath::Abs(rMin < 0.04) && (charge == charge2)) 
	      {

		  if (TMath::Abs(dpt)/ptm  < .1) {
		      fhCl1Cl2G->Fill(nClustersTPC, nClustersTPC2);
		  }
		  
		  if (TMath::Abs(dpt)/ptm  > .5) fhCl1Cl2B->Fill(nClustersTPC, nClustersTPC2);
		  fhPh1Ph2->Fill(track->Phi(), track2->Phi());


		  Double_t sigmaPt1 = TMath::Sqrt(track->GetSigma1Pt2());
		  Double_t sigmaPt2 = TMath::Sqrt(track->GetSigma1Pt2());
		  Double_t sigmaPt = 0.5 * TMath::Sqrt(sigmaPt1 * sigmaPt1 + sigmaPt2 * sigmaPt2);
		  if (TMath::Abs(dpt)/ptm > 0.2 && pt > 10. &&  charge > 0.)
//		      printf("Track#2 %5d %5d %13.3f %13.3f %13.3f %13.3f\n", Entry(), jMin, track2->Phi(), dz, dpt, zOut2);	      		  
		  if (TMath::Abs(zOut) > 1.e-10 && TMath::Abs(zOut2) > 1.e-10) fhDZvsZ->Fill(TMath::Abs(zOut) * zOut/zOut2, dz);
		  if (zOut * zOut2 > 0. && zOut < 0) fhDZvsPhi->Fill(phi, dz);
		  
		  fhCh1Ch2->Fill(charge, track2->Charge());
		  
		  
		  if (charge > 0.) {
		      fhDPt[kPosC]->Fill(dpt);
		      fhD1ovPt[kPosC]->Fill(d1pt);
		      fpDPt[kPosC]->Fill(ptm, d1pt);
		      if (ptm > 5. && ptm < 15.) fhDPtovPt[kPosC]->Fill(dpt/ptm);
		      fhDZ[kPosC]->Fill(dz);
		      fhDX[kPosC]->Fill(dx);
		      fhDY[kPosC]->Fill(dy);
		      fpDPtS[kPosC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  } else {
		      fhDPt[kNegC]->Fill(dpt);
		      fhD1ovPt[kNegC]->Fill(d1pt);
		      fpDPt[kNegC]->Fill(ptm, d1pt);
		      if (ptm > 5. && ptm < 15.) fhDPtovPt[kNegC]->Fill(dpt/ptm);
		      fhDZ[kNegC]->Fill(dz);
		      fhDX[kNegC]->Fill(dx);
		      fhDY[kNegC]->Fill(dy);
		      fpDPtS[kNegC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  }

		  if (z > 0.) {
		      fhDPt[kPosZ]->Fill(dpt);
		      fhD1ovPt[kPosZ]->Fill(d1pt);
		      fpDPt[kPosZ]->Fill(ptm, d1pt);
		      fhDZ[kPosZ]->Fill(dz);
		      fhDX[kPosZ]->Fill(dx);
		      fhDY[kPosZ]->Fill(dy);
		      fpDPtS[kPosZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  } else {
		      fhDPt[kNegZ]->Fill(dpt);
		      fhD1ovPt[kNegZ]->Fill(d1pt);
		      fpDPt[kNegZ]->Fill(ptm, d1pt);
		      fhDZ[kNegZ]->Fill(dz);
		      fhDX[kNegZ]->Fill(dx);
		      fhDY[kNegZ]->Fill(dy);
		      fpDPtS[kNegZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  } 


		  if (TMath::Abs(dpt)/ptm > 0.5) {
		      fhDPt[kBad]->Fill(dpt);
		      fpDPt[kBad]->Fill(ptm, TMath::Abs(dpt)/ptm);
		      fhDZ[kBad]->Fill(dz);
		      fhDX[kBad]->Fill(dx);
		      fhDY[kBad]->Fill(dy);
		      fpDPtS[kBad]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  } else {
		      fhDPt[kGood]->Fill(dpt);
		      fpDPt[kGood]->Fill(ptm, TMath::Abs(dpt)/ptm);
		      fhDZ[kGood]->Fill(dz);
		      fhDX[kGood]->Fill(dx);
		      fhDY[kGood]->Fill(dy);
		      fpDPtS[kGood]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
		  } 
		  
	      } // good matches
	  } // found possible match
      } // upper
  } // tracks 1
  PostData(1, fHists);
}      

//________________________________________________________________________
void AliAnalysisTaskCosmic::Terminate(Option_t *) 
{
}
Commit	Line	Data
13698385	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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	/* $Id$ */
	17
	18	// Analysis Task for the Quality Assurance of Cosmic Data
	19	// Two track segments in the are matched in angle and charged. T
	20	// The quality of the matching in is checked by comparing
	21	// the tarnsverse momenta and starting points of the track segments
	22	//
	23	// Author
	24	// Andreas Morsch
	25	// andreas.morsch@cern.ch
	26
	27	#include "TChain.h"
	28	#include "TTree.h"
0ee094cc	29
c2174181	30	#include "TH1F.h"
	31	#include "TH2F.h"
	32	#include "TProfile.h"
	33	#include "TList.h"
c2174181	34
	35	#include "AliESDEvent.h"
	36	#include "AliESDtrack.h"
	37	#include "AliAnalysisTaskCosmic.h"
c2174181	38
	39	ClassImp(AliAnalysisTaskCosmic)
	40
	41	//________________________________________________________________________
	42	AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const char *name)
	43	: AliAnalysisTaskSE(name),
	44	fHists(0),
	45	fhDZvsZ(0),
	46	fhDZvsPhi(0),
	47	fhCh1Ch2(0),
	48	fhPh1Ph2(0),
	49	fhCl1Cl2G(0),
	50	fhCl1Cl2B(0)
	51	{
	52	//
	53	// Constructor
	54	//
13698385	55
c2174181	56	for (Int_t i = 0; i < 6; i++) {
	57	fhPt[i] = 0;
	58	fhTheta[i] = 0;
	59	fhPhi[i] = 0;
	60	fhDPhi[i] = 0;
	61	fhDTheta[i] = 0;
	62	fhDZ[i] = 0;
	63	fhDX[i] = 0;
	64	fhDY[i] = 0;
	65	fhDPt[i] = 0;
	66	fhD1ovPt[i] = 0;
	67	fpDPt[i] = 0;
	68	fhDPtovPt[i] = 0;
	69	fpDPtS[i] = 0;
	70	}
	71
	72	DefineOutput(1, TList::Class());
	73	}
	74
13698385	75	AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const AliAnalysisTaskCosmic& task)
	76	: AliAnalysisTaskSE(task),
	77	fHists(0),
	78	fhDZvsZ(0),
	79	fhDZvsPhi(0),
	80	fhCh1Ch2(0),
	81	fhPh1Ph2(0),
	82	fhCl1Cl2G(0),
	83	fhCl1Cl2B(0)
	84	{
	85	// Copy constructor
	86	fHists = new TList();
	87	fhDZvsZ = (TH2F*) ((task.fhDZvsZ) ->Clone());
	88	fhDZvsPhi = (TH2F*) ((task.fhDZvsPhi) ->Clone());
	89	fhCh1Ch2 = (TH2F*) ((task.fhCh1Ch2) ->Clone());
	90	fhPh1Ph2 = (TH2F*) ((task.fhPh1Ph2) ->Clone());
	91	fhCl1Cl2G = (TH2F*) ((task.fhCl1Cl2G) ->Clone());
	92	fhCl1Cl2B = (TH2F*) ((task.fhCl1Cl2B) ->Clone());
	93	for (Int_t i = 0; i < 6; i++) {
	94	fhPt[i] = (TH1F*) ((task.fhPt[i]) ->Clone());
	95	fhTheta[i] = (TH1F*) ((task.fhTheta[i]) ->Clone());
	96	fhPhi[i] = (TH1F*) ((task.fhPhi[i]) ->Clone());
	97	fhDPhi[i] = (TH1F*) ((task.fhDPhi[i]) ->Clone());
	98	fhDTheta[i] = (TH1F*) ((task.fhDTheta[i]) ->Clone());
	99	fhDZ[i] = (TH1F*) ((task.fhDZ[i]) ->Clone());
	100	fhDX[i] = (TH1F*) ((task.fhDX[i]) ->Clone());
	101	fhDY[i] = (TH1F*) ((task.fhDY[i]) ->Clone());
	102	fhDPt[i] = (TH1F*) ((task.fhDPt[i]) ->Clone());
	103	fhD1ovPt[i] = (TH1F*) ((task.fhD1ovPt[i]) ->Clone());
	104	fhDPtovPt[i] = (TH1F*) ((task.fhDPtovPt[i]) ->Clone());
	105
	106	fpDPt[i] = (TProfile*) ((task.fpDPt[i]) ->Clone());
	107	fpDPtS[i] = (TProfile*) ((task.fpDPtS[i]) ->Clone());
	108	}
	109	}
	110
	111	//______________________________________________________________________________
	112	AliAnalysisTaskCosmic& AliAnalysisTaskCosmic::operator=(const AliAnalysisTaskCosmic& task)
	113	{
	114	// Assignment operator
	115	if(this!=&task) {
	116
	117	AliAnalysisTaskSE::operator=(task);
	118
	119	fHists = new TList();
	120	fhDZvsZ = (TH2F*) ((task.fhDZvsZ) ->Clone());
	121	fhDZvsPhi = (TH2F*) ((task.fhDZvsPhi) ->Clone());
	122	fhCh1Ch2 = (TH2F*) ((task.fhCh1Ch2) ->Clone());
	123	fhPh1Ph2 = (TH2F*) ((task.fhPh1Ph2) ->Clone());
	124	fhCl1Cl2G = (TH2F*) ((task.fhCl1Cl2G) ->Clone());
	125	fhCl1Cl2B = (TH2F*) ((task.fhCl1Cl2B) ->Clone());
	126	for (Int_t i = 0; i < 6; i++) {
	127	fhPt[i] = (TH1F*) ((task.fhPt[i]) ->Clone());
	128	fhTheta[i] = (TH1F*) ((task.fhTheta[i]) ->Clone());
	129	fhPhi[i] = (TH1F*) ((task.fhPhi[i]) ->Clone());
	130	fhDPhi[i] = (TH1F*) ((task.fhDPhi[i]) ->Clone());
	131	fhDTheta[i] = (TH1F*) ((task.fhDTheta[i]) ->Clone());
	132	fhDZ[i] = (TH1F*) ((task.fhDZ[i]) ->Clone());
	133	fhDX[i] = (TH1F*) ((task.fhDX[i]) ->Clone());
	134	fhDY[i] = (TH1F*) ((task.fhDY[i]) ->Clone());
	135	fhDPt[i] = (TH1F*) ((task.fhDPt[i]) ->Clone());
	136	fhD1ovPt[i] = (TH1F*) ((task.fhD1ovPt[i]) ->Clone());
	137	fhDPtovPt[i] = (TH1F*) ((task.fhDPtovPt[i]) ->Clone());
	138
139	fpDPt[i] = (TProfile*) ((task.fpDPt[i]) ->Clone());
140	fpDPtS[i] = (TProfile*) ((task.fpDPtS[i]) ->Clone());
141	}
142	}
143	return *this;
144	}
145
c2174181	146
	147	//________________________________________________________________________
	148	void AliAnalysisTaskCosmic::UserCreateOutputObjects()
	149	{
	150	// Create histograms
	151	// Called once
	152	TString ext[6] = {"PC", "NC", "PZ", "NZ", "_Good", "_Bad"};
	153
	154	char name[12];
	155
	156	fHists = new TList();
	157
	158	for (Int_t i = 0; i < 6; i++) {
	159	// Pt
8df6b7b4	160	snprintf(name, 12, "fhPt%2s", ext[i].Data());
8df6b7b4	161	fhPt[i] = new TH1F(name, " pT distribution", 800, 0., 200.);
c2174181	162	fhPt[i]->SetXTitle("p_{T} [Gev]");
c2174181	163	// Phi
8df6b7b4	164	snprintf(name, 12, "fhPhi%2s", ext[i].Data());
c2174181	165	fhPhi[i] = new TH1F(name, "Phi distribution", 62, 0., 2. * TMath::Pi());
	166	fhPhi[i]->SetXTitle("#phi [rad]");
	167	// Theta
8df6b7b4	168	snprintf(name, 12, "fhTheta%2s", ext[i].Data());
c2174181	169	fhTheta[i] = new TH1F(name, "Theta distribution", 62, 0., TMath::Pi());
	170	fhTheta[i]->SetXTitle("#theta [rad]");
	171	// Delta Phi
8df6b7b4	172	snprintf(name, 12, "fhDPhi%2s", ext[i].Data());
c2174181	173	fhDPhi[i] = new TH1F(name, "DeltaPhi distribution", 320, -0.4, 0.4);
	174	fhDPhi[i]->SetXTitle("#Delta#phi [rad]");
	175	// Delta Theta
8df6b7b4	176	snprintf(name, 12, "fhDTheta%2s", ext[i].Data());
c2174181	177	fhDTheta[i] = new TH1F(name, "DeltaTheta distribution", 320, -0.4, 0.4);
	178	fhDTheta[i]->SetXTitle("#Delta#theta [rad]");
	179	// Delta Z
8df6b7b4	180	snprintf(name, 12, "fhDZ%2s", ext[i].Data());
c2174181	181	fhDZ[i] = new TH1F(name, "DeltaZ distribution", 200, -10., 10.);
	182	fhDZ[i]->SetXTitle("#DeltaZ [cm]");
	183	// Delta X
8df6b7b4	184	snprintf(name, 12, "fhDX%2s", ext[i].Data());
c2174181	185	fhDX[i] = new TH1F(name, "DeltaX distribution", 200, -10., 10.);
	186	fhDX[i]->SetXTitle("#DeltaX [cm]");
	187	// Delta Y
8df6b7b4	188	snprintf(name, 12, "fhDY%2s", ext[i].Data());
c2174181	189	fhDY[i] = new TH1F(name, "DeltaY distribution", 200, -10, 10.);
	190	fhDY[i]->SetXTitle("#DeltaY [cm]");
	191	// Delta Pt
8df6b7b4	192	snprintf(name, 12, "fhDPt%2s", ext[i].Data());
c2174181	193	fhDPt[i] = new TH1F(name, "DeltaPt distribution", 200, -20., 20.);
	194	fhDPt[i]->SetXTitle("#Delta p_{T} [GeV]");
	195
	196	// Delta 1/Pt
8df6b7b4	197	snprintf(name, 12, "fhD1ovPt%2s", ext[i].Data());
c2174181	198	fhD1ovPt[i] = new TH1F(name, "Delta 1/Pt distribution", 200, -1., 1.);
	199	fhD1ovPt[i]->SetXTitle("#Delta 1/Pt");
	200
	201	// Delta Pt over Pt
8df6b7b4	202	snprintf(name, 12, "fhDPtovPt%2s", ext[i].Data());
c2174181	203	fhDPtovPt[i] = new TH1F(name, "DeltaPt/Pt distribution", 200, -2., 2.);
	204	fhDPtovPt[i]->SetXTitle("#DeltaPt/Pt");
	205
	206
	207
	208	// Delta Pt/ Pt vs Pt
8df6b7b4	209	snprintf(name, 12, "fpDPt%2s", ext[i].Data());
c2174181	210	fpDPt[i] = new TProfile(name, "#Delta Pt / Pt", 20, 0., 20., -1, 1., "S");
	211	fpDPt[i]->SetXTitle("p_{T} [GeV]");
	212	fpDPt[i]->SetYTitle("#Delta 1/p_{T} [GeV^{-1}]");
	213	// Delta Pt error
8df6b7b4	214	snprintf(name, 12, "fpDPtS%2s", ext[i].Data());
c2174181	215	fpDPtS[i] = new TProfile(name, "#Delta Pt / <sigma>", 20, 0., 20., 0., 10.);
	216	fpDPtS[i]->SetXTitle("p_{T}");
	217	fpDPtS[i]->SetYTitle("#Delta p_{T} / <#sigma_{p_{T}}>");
	218
	219
	220	fHists->Add(fhPt[i]);
	221	fHists->Add(fhPhi[i]);
	222	fHists->Add(fhTheta[i]);
	223	fHists->Add(fhDPhi[i]);
	224	fHists->Add(fhDPt[i]);
	225	fHists->Add(fhD1ovPt[i]);
	226	fHists->Add(fhDPtovPt[i]);
	227	fHists->Add(fhDTheta[i]);
	228	fHists->Add(fhDZ[i]);
	229	fHists->Add(fhDX[i]);
	230	fHists->Add(fhDY[i]);
	231	fHists->Add(fpDPt[i]);
	232	fHists->Add(fpDPtS[i]);
	233	}
	234
	235	fhDZvsZ = new TH2F("fhDZvsZ", "dz vs z", 500, -250., 250., 100, -10., 10.);
	236	fhDZvsZ->SetXTitle("z_{in} * sign(z_{in}) * sign(z_{out}) [cm]");
	237	fhDZvsZ->SetYTitle("#Delta z [cm]");
	238
	239
	240	fhDZvsPhi = new TH2F("fhDZvsPhi", "dz vs phi", 36, 0., TMath::Pi(), 50, -2., 2.);
	241
	242	fhCh1Ch2 = new TH2F("fCh1Ch2", "ch1 vs ch2", 8, -2., 2., 8, -2., 2.);
	243	fhPh1Ph2 = new TH2F("fPh1Ph2", "ph1 vs ph2", 128, 0., 2. * TMath::Pi(), 128, 0., 2. * TMath::Pi());
	244	fhCl1Cl2G = new TH2F("fCl1Cl2G", "#cl vs #cl", 200, 0., 200., 200, 0., 200);
	245	fhCl1Cl2B = new TH2F("fCl1Cl2B", "#cl vs #cl", 200, 0., 200., 200, 0., 200);
	246
	247	fHists->Add(fhDZvsZ);
	248	fHists->Add(fhDZvsPhi);
	249	fHists->Add(fhCh1Ch2);
	250	fHists->Add(fhPh1Ph2);
	251	fHists->Add(fhCl1Cl2G);
	252	fHists->Add(fhCl1Cl2B);
	253	fHists->SetOwner();
	254
	255	}
	256
	257	//________________________________________________________________________
	258	void AliAnalysisTaskCosmic::UserExec(Option_t *)
	259	{
	260	// Main loop
	261	// Called for each event
	262
	263	if (!fInputEvent) {
	264	Printf("ERROR: fESD not available");
	265	return;
	266	}
	267
	268	AliESDEvent* esdE = (AliESDEvent*)fInputEvent;
86e75636	269
	270	if (esdE->GetNumberOfTracks() != 2) {
	271	PostData(1, fHists);
	272	return;
	273	}
	274
c2174181	275
	276	for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
	277	AliESDtrack* track = esdE->GetTrack(iTracks);
	278
	279
	280
	281	// const AliExternalTrackParam * track = trackG->GetTPCInnerParam();
	282	// if (!track) continue;
	283	// Pt cut
	284	Float_t pt = track->Pt();
0ee094cc	285	Short_t charge = track->Charge();
c2174181	286	Float_t phi = track->Phi();
13698385	287	if (phi > 0. && phi < TMath::Pi()) charge*=(-1);
c2174181	288
	289
	290
	291	// TPC track
	292	UInt_t status = track->GetStatus();
0ee094cc	293	if ((status&AliESDtrack::kTPCrefit) ==0) continue;
c2174181	294
	295	Int_t nClustersTPC = track->GetTPCclusters(0);
	296	if (nClustersTPC < 50) continue;
	297
	298	//
	299
	300	Float_t z = track->GetZ();
	301	Float_t x = track->Xv();
	302	Float_t y = track->Yv();
	303	Float_t theta = track->Theta();
	304
	305
	306	const AliExternalTrackParam * trackOut = track->GetOuterParam();
	307	Float_t zOut = 0.;
	308	if (trackOut)zOut = trackOut->Zv();
	309
	310
	311	if (charge > 0) {
	312	fhPt[kPosC] ->Fill(pt);
	313	} else {
	314	fhPt[kNegC] ->Fill(pt);
	315	}
	316
	317	// if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (pt < 1. \|\| pt > 50.)) continue;
	318
	319
	320	//
	321	if (charge > 0) {
	322	fhPhi[kPosC] ->Fill(phi);
	323	fhTheta[kPosC]->Fill(theta);
	324	} else {
	325	fhPhi[kNegC] ->Fill(phi);
	326	fhTheta[kNegC]->Fill(theta);
	327	}
	328
	329
	330	if (z > 0) {
	331	fhPt[kPosZ] ->Fill(pt);
	332	fhPhi[kPosZ] ->Fill(phi);
	333	fhTheta[kPosZ]->Fill(theta);
	334	} else {
	335	fhPt[kNegZ] ->Fill(pt);
	336	fhPhi[kNegZ] ->Fill(phi);
	337	fhTheta[kNegZ]->Fill(theta);
	338	}
	339
	340
	341
	342	// Tracks coming from above
	343	if (phi > 0. && phi < TMath::Pi()) {
	344
	345	// printf("Track#1 %5d %5d %13.3f %13.3f \n", Int_t(Entry()), iTracks, phi, zOut);
	346
	347	Float_t dphiMin = 999.;
	348	Float_t rMin = 999.;
	349	Int_t jMin = -1;
	350	// Search for a matching track (in dphi-dtheta)
	351	for (Int_t jTracks = 0; jTracks < esdE->GetNumberOfTracks(); jTracks++) {
	352	if (jTracks == iTracks) continue;
	353
	354	AliESDtrack* track2 = esdE->GetTrack(jTracks);
	355
113bb16c	356	UInt_t status2 = track2->GetStatus();
113bb16c	357	if ((status2&AliESDtrack::kTPCrefit) ==0) continue;
c2174181	358
	359	Int_t nClustersTPC2 = track2->GetTPCclusters(0);
	360	if (nClustersTPC2 < 50) continue;
	361	// if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (track2->Pt() < 1. \|\| track2->Pt() > 50.)) continue;
	362	// const AliExternalTrackParam * track2 = trackG2->GetTPCInnerParam();
	363	// if (!track2) continue;
	364
	365	Float_t phi2 = track2->Phi() - TMath::Pi();
	366	Float_t theta2 = TMath::Pi() - track2->Theta();
	367	Float_t dphi = phi2 - phi;
	368	Float_t dtheta = theta2 - theta;
	369
	370	if (dphi > TMath::Pi()) dphi -= 2. * TMath::Pi();
	371	if (dphi < -TMath::Pi()) dphi += 2. * TMath::Pi();
	372
	373	Float_t dR = TMath::Sqrt(dphi * dphi + dtheta * dtheta);
	374
	375	if (dR < rMin) {
	376	rMin = dR;
	377	dphiMin = dphi;
	378	jMin = jTracks;
	379	}
	380
	381	} // tracks 2
	382	if (jMin != -1) {
	383	// we found a matching track candidate ...
	384	AliESDtrack* track2 = esdE->GetTrack(jMin);
	385	const AliExternalTrackParam * trackOut2 = track2->GetOuterParam();
	386	Float_t zOut2 = 0.;
	387	if (trackOut2) zOut2 = trackOut2->Zv();
	388
	389
	390	Float_t theta2 = - track2->Theta() + TMath::Pi();
	391	Float_t z2 = track2->GetZ();
	392	Float_t x2 = track2->Xv();
	393	Float_t y2 = track2->Yv();
0ee094cc	394	Short_t charge2 = track2->Charge();
c2174181	395	Float_t dz = z2 - z;
	396	Float_t dx = x2 - x;
	397	Float_t dy = y2 - y;
	398	Float_t pt1 = track->Pt();
	399	Float_t pt2 = track2->Pt();
	400	Float_t dpt = pt2 - pt1;
	401	Float_t d1pt = 1./pt2 - 1./pt1;
	402
	403	Float_t ptm = 0.5 * (pt1 + pt2);
	404
	405	Int_t nClustersTPC2 = track2->GetTPCclusters(0);
	406
	407	if (charge > 0.) {
	408	fhDPhi[kPosC] ->Fill(dphiMin);
	409	fhDTheta[kPosC]->Fill(theta2 - theta);
	410	} else {
	411	fhDPhi[kNegC] ->Fill(dphiMin);
	412	fhDTheta[kNegC]->Fill(theta2 - theta);
	413	}
	414
	415	if (z > 0.) {
	416	fhDPhi[kPosZ] ->Fill(dphiMin);
	417	fhDTheta[kPosZ]->Fill(theta2 - theta);
	418	} else {
	419	fhDPhi[kNegZ] ->Fill(dphiMin);
	420	fhDTheta[kNegZ]->Fill(theta2 - theta);
	421	}
	422
	423	if (TMath::Abs(dpt)/ptm > 0.5) {
	424	fhDPhi[kBad] ->Fill(dphiMin);
	425	fhDTheta[kBad]->Fill(theta2 - theta);
	426	} else {
	427	fhDPhi[kGood] ->Fill(dphiMin);
	428	fhDTheta[kGood]->Fill(theta2 - theta);
	429	}
	430	// Good matches ...
	431	// if (TMath::Abs(rMin < 0.04) && (charge == charge2) && TMath::Abs(dz) < 0.5 && TMath::Abs(dy) && TMath::Abs(dx) < 0.5 )
	432	// if (TMath::Abs(rMin < 0.04) && (charge == charge2) && (zOut * zOut2) < 0.)
	433	if (TMath::Abs(rMin < 0.04) && (charge == charge2))
	434	{
	435
	436	if (TMath::Abs(dpt)/ptm < .1) {
	437	fhCl1Cl2G->Fill(nClustersTPC, nClustersTPC2);
	438	}
	439
	440	if (TMath::Abs(dpt)/ptm > .5) fhCl1Cl2B->Fill(nClustersTPC, nClustersTPC2);
	441	fhPh1Ph2->Fill(track->Phi(), track2->Phi());
	442
	443
	444	Double_t sigmaPt1 = TMath::Sqrt(track->GetSigma1Pt2());
	445	Double_t sigmaPt2 = TMath::Sqrt(track->GetSigma1Pt2());
	446	Double_t sigmaPt = 0.5 * TMath::Sqrt(sigmaPt1 * sigmaPt1 + sigmaPt2 * sigmaPt2);
	447	if (TMath::Abs(dpt)/ptm > 0.2 && pt > 10. && charge > 0.)
	448	// printf("Track#2 %5d %5d %13.3f %13.3f %13.3f %13.3f\n", Entry(), jMin, track2->Phi(), dz, dpt, zOut2);
13698385	449	if (TMath::Abs(zOut) > 1.e-10 && TMath::Abs(zOut2) > 1.e-10) fhDZvsZ->Fill(TMath::Abs(zOut) * zOut/zOut2, dz);
c2174181	450	if (zOut * zOut2 > 0. && zOut < 0) fhDZvsPhi->Fill(phi, dz);
	451
	452	fhCh1Ch2->Fill(charge, track2->Charge());
	453
	454
	455	if (charge > 0.) {
	456	fhDPt[kPosC]->Fill(dpt);
	457	fhD1ovPt[kPosC]->Fill(d1pt);
	458	fpDPt[kPosC]->Fill(ptm, d1pt);
	459	if (ptm > 5. && ptm < 15.) fhDPtovPt[kPosC]->Fill(dpt/ptm);
	460	fhDZ[kPosC]->Fill(dz);
	461	fhDX[kPosC]->Fill(dx);
	462	fhDY[kPosC]->Fill(dy);
	463	fpDPtS[kPosC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	464	} else {
	465	fhDPt[kNegC]->Fill(dpt);
	466	fhD1ovPt[kNegC]->Fill(d1pt);
	467	fpDPt[kNegC]->Fill(ptm, d1pt);
	468	if (ptm > 5. && ptm < 15.) fhDPtovPt[kNegC]->Fill(dpt/ptm);
	469	fhDZ[kNegC]->Fill(dz);
	470	fhDX[kNegC]->Fill(dx);
	471	fhDY[kNegC]->Fill(dy);
	472	fpDPtS[kNegC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	473	}
	474
	475	if (z > 0.) {
	476	fhDPt[kPosZ]->Fill(dpt);
	477	fhD1ovPt[kPosZ]->Fill(d1pt);
	478	fpDPt[kPosZ]->Fill(ptm, d1pt);
	479	fhDZ[kPosZ]->Fill(dz);
	480	fhDX[kPosZ]->Fill(dx);
	481	fhDY[kPosZ]->Fill(dy);
	482	fpDPtS[kPosZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	483	} else {
	484	fhDPt[kNegZ]->Fill(dpt);
	485	fhD1ovPt[kNegZ]->Fill(d1pt);
	486	fpDPt[kNegZ]->Fill(ptm, d1pt);
	487	fhDZ[kNegZ]->Fill(dz);
	488	fhDX[kNegZ]->Fill(dx);
	489	fhDY[kNegZ]->Fill(dy);
	490	fpDPtS[kNegZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	491	}
	492
	493
	494	if (TMath::Abs(dpt)/ptm > 0.5) {
	495	fhDPt[kBad]->Fill(dpt);
	496	fpDPt[kBad]->Fill(ptm, TMath::Abs(dpt)/ptm);
	497	fhDZ[kBad]->Fill(dz);
	498	fhDX[kBad]->Fill(dx);
	499	fhDY[kBad]->Fill(dy);
	500	fpDPtS[kBad]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	501	} else {
	502	fhDPt[kGood]->Fill(dpt);
	503	fpDPt[kGood]->Fill(ptm, TMath::Abs(dpt)/ptm);
	504	fhDZ[kGood]->Fill(dz);
	505	fhDX[kGood]->Fill(dx);
	506	fhDY[kGood]->Fill(dy);
	507	fpDPtS[kGood]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
	508	}
	509
	510	} // good matches
	511	} // found possible match
	512	} // upper
	513	} // tracks 1
514	PostData(1, fHists);
515	}
516
517	//________________________________________________________________________
518	void AliAnalysisTaskCosmic::Terminate(Option_t *)
519	{
520	}