* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
-/* $Id$ */
-
//____________________________________________________________________
-//
+//
// AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets.
-//
-// This class has been developed from AliITSMultReconstructor (see
-// it for more details). It is the class for the Trackleter used to estimate
-// SPD plane efficiency.
+//
+// This class has been developed from AliITSMultReconstructor (see
+// it for more details). It is the class for the Trackleter used to estimate
+// SPD plane efficiency.
// The trackleter prediction is built using the vertex and 1 cluster.
//
-//
+//
// Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili
// email: giuseppe.bruno@ba.infn.it
-//
+//
//____________________________________________________________________
+/* $Id$ */
+
#include <TFile.h>
+#include <TTree.h>
#include <TParticle.h>
+#include <TParticlePDG.h>
#include <TSystem.h>
#include <Riostream.h>
+#include <TClonesArray.h>
-#include "AliITSMultReconstructor.h"
+#include "AliTracker.h"
#include "AliITSTrackleterSPDEff.h"
#include "AliITSgeomTGeo.h"
#include "AliLog.h"
#include "AliITSPlaneEffSPD.h"
#include "AliStack.h"
-
+#include "AliTrackReference.h"
+#include "AliRunLoader.h"
+#include "AliITSReconstructor.h"
+#include "AliITSRecPoint.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
//____________________________________________________________________
ClassImp(AliITSTrackleterSPDEff)
//____________________________________________________________________
AliITSTrackleterSPDEff::AliITSTrackleterSPDEff():
-AliITSMultReconstructor(),
+AliTracker(),
+//
+fClustersLay1(0),
+fClustersLay2(0),
+fTracklets(0),
+fAssociationFlag(0),
+fNClustersLay1(0),
+fNClustersLay2(0),
+fNTracklets(0),
+fOnlyOneTrackletPerC2(0),
+fPhiWindowL2(0),
+fZetaWindowL2(0),
+fPhiOverlapCut(0),
+fZetaOverlapCut(0),
+fHistOn(0),
+fhClustersDPhiAcc(0),
+fhClustersDThetaAcc(0),
+fhClustersDZetaAcc(0),
+fhClustersDPhiAll(0),
+fhClustersDThetaAll(0),
+fhClustersDZetaAll(0),
+fhDPhiVsDThetaAll(0),
+fhDPhiVsDThetaAcc(0),
+fhDPhiVsDZetaAll(0),
+fhDPhiVsDZetaAcc(0),
+fhetaTracklets(0),
+fhphiTracklets(0),
+fhetaClustersLay1(0),
+fhphiClustersLay1(0),
+//
fAssociationFlag1(0),
fChipPredOnLay2(0),
fChipPredOnLay1(0),
fPhiWindowL1(0),
fZetaWindowL1(0),
fOnlyOneTrackletPerC1(0),
+fUpdateOncePerEventPlaneEff(0),
+fMinContVtx(0),
+fChipUpdatedInEvent(0),
fPlaneEffSPD(0),
+fPlaneEffBkg(0),
+fReflectClusterAroundZAxisForLayer0(kFALSE),
+fReflectClusterAroundZAxisForLayer1(kFALSE),
+fLightBkgStudyInParallel(kFALSE),
fMC(0),
fUseOnlyPrimaryForPred(0),
fUseOnlySecondaryForPred(0),
fUseOnlySameParticle(0),
fUseOnlyDifferentParticle(0),
-fUseOnlyStableParticle(0),
+fUseOnlyStableParticle(1),
fPredictionPrimary(0),
fPredictionSecondary(0),
fClusterPrimary(0),
fClusterSecondary(0),
+fSuccessPP(0),
+fSuccessTT(0),
+fSuccessS(0),
+fSuccessP(0),
+fFailureS(0),
+fFailureP(0),
+fRecons(0),
+fNonRecons(0),
fhClustersDPhiInterpAcc(0),
fhClustersDThetaInterpAcc(0),
fhClustersDZetaInterpAcc(0),
fhDPhiVsDZetaInterpAll(0),
fhDPhiVsDZetaInterpAcc(0),
fhetaClustersLay2(0),
-fhphiClustersLay2(0)
+fhphiClustersLay2(0),
+fhClustersInChip(0),
+fhClustersInModuleLay1(0),
+fhClustersInModuleLay2(0)
{
-
- // Method to check the SPD chips efficiencies by using tracklets
-
-
+ // default constructor
+// from AliITSMultReconstructor
+ Init();
+}
+//______________________________________________________________________
+void AliITSTrackleterSPDEff::Init() {
+ SetPhiWindowL2();
+ SetZetaWindowL2();
+ SetOnlyOneTrackletPerC2();
+ fClustersLay1 = new Float_t*[300000];
+ fClustersLay2 = new Float_t*[300000];
+ fTracklets = new Float_t*[300000];
+ fAssociationFlag = new Bool_t[300000];
+//
SetPhiWindowL1();
SetZetaWindowL1();
SetOnlyOneTrackletPerC1();
fAssociationFlag1 = new Bool_t[300000];
fChipPredOnLay2 = new UInt_t[300000];
fChipPredOnLay1 = new UInt_t[300000];
+ fChipUpdatedInEvent = new Bool_t[1200];
for(Int_t i=0; i<300000; i++) {
+ // from AliITSMultReconstructor
+ fClustersLay1[i] = new Float_t[6];
+ fClustersLay2[i] = new Float_t[6];
+ fTracklets[i] = new Float_t[5];
+ fAssociationFlag[i] = kFALSE;
+ //
fAssociationFlag1[i] = kFALSE;
}
+ for(Int_t i=0;i<1200; i++) fChipUpdatedInEvent[i] = kFALSE;
if (GetHistOn()) BookHistos();
fPlaneEffSPD = new AliITSPlaneEffSPD();
+ SetLightBkgStudyInParallel();
}
//______________________________________________________________________
-AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr),
+AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) :
+AliTracker(mr),
+// from AliITSMultReconstructor
+fClustersLay1(mr.fClustersLay1),
+fClustersLay2(mr.fClustersLay2),
+fTracklets(mr.fTracklets),
+fAssociationFlag(mr.fAssociationFlag),
+fNClustersLay1(mr.fNClustersLay1),
+fNClustersLay2(mr.fNClustersLay2),
+fNTracklets(mr.fNTracklets),
+fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
+fPhiWindowL2(mr.fPhiWindowL2),
+fZetaWindowL2(mr.fZetaWindowL2),
+fPhiOverlapCut(mr.fPhiOverlapCut),
+fZetaOverlapCut(mr.fZetaOverlapCut),
+fHistOn(mr.fHistOn),
+fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
+fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
+fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
+fhClustersDPhiAll(mr.fhClustersDPhiAll),
+fhClustersDThetaAll(mr.fhClustersDThetaAll),
+fhClustersDZetaAll(mr.fhClustersDZetaAll),
+fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
+fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
+fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
+fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
+fhetaTracklets(mr.fhetaTracklets),
+fhphiTracklets(mr.fhphiTracklets),
+fhetaClustersLay1(mr.fhetaClustersLay1),
+fhphiClustersLay1(mr.fhphiClustersLay1),
+//
fAssociationFlag1(mr.fAssociationFlag1),
fChipPredOnLay2(mr.fChipPredOnLay2),
fChipPredOnLay1(mr.fChipPredOnLay1),
fPhiWindowL1(mr.fPhiWindowL1),
fZetaWindowL1(mr.fZetaWindowL1),
fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1),
+fUpdateOncePerEventPlaneEff(mr.fUpdateOncePerEventPlaneEff),
+fMinContVtx(mr.fMinContVtx),
+fChipUpdatedInEvent(mr.fChipUpdatedInEvent),
fPlaneEffSPD(mr.fPlaneEffSPD),
+fPlaneEffBkg(mr.fPlaneEffBkg),
+fReflectClusterAroundZAxisForLayer0(mr.fReflectClusterAroundZAxisForLayer0),
+fReflectClusterAroundZAxisForLayer1(mr.fReflectClusterAroundZAxisForLayer1),
+fLightBkgStudyInParallel(mr.fLightBkgStudyInParallel),
fMC(mr.fMC),
fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred),
fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred),
fPredictionSecondary(mr.fPredictionSecondary),
fClusterPrimary(mr.fClusterPrimary),
fClusterSecondary(mr.fClusterSecondary),
+fSuccessPP(mr.fSuccessPP),
+fSuccessTT(mr.fSuccessTT),
+fSuccessS(mr.fSuccessS),
+fSuccessP(mr.fSuccessP),
+fFailureS(mr.fFailureS),
+fFailureP(mr.fFailureP),
+fRecons(mr.fRecons),
+fNonRecons(mr.fNonRecons),
fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc),
fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc),
fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc),
fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll),
fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc),
fhetaClustersLay2(mr.fhetaClustersLay2),
-fhphiClustersLay2(mr.fhphiClustersLay2)
+fhphiClustersLay2(mr.fhphiClustersLay2),
+fhClustersInChip(mr.fhClustersInChip),
+fhClustersInModuleLay1(mr.fhClustersInModuleLay1),
+fhClustersInModuleLay2(mr.fhClustersInModuleLay2)
{
// Copy constructor
}
//______________________________________________________________________
AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){
// Destructor
-
+// from AliITSMultReconstructor
+ // delete arrays
+ for(Int_t i=0; i<300000; i++) {
+ delete [] fClustersLay1[i];
+ delete [] fClustersLay2[i];
+ delete [] fTracklets[i];
+ }
+ delete [] fClustersLay1;
+ delete [] fClustersLay2;
+ delete [] fTracklets;
+ delete [] fAssociationFlag;
+//
// delete histograms
DeleteHistos();
delete [] fChipPredOnLay2;
delete [] fChipPredOnLay1;
+ delete [] fChipUpdatedInEvent;
+
delete [] fPredictionPrimary;
delete [] fPredictionSecondary;
delete [] fClusterPrimary;
delete [] fClusterSecondary;
+ delete [] fSuccessPP;
+ delete [] fSuccessTT;
+ delete [] fSuccessS;
+ delete [] fSuccessP;
+ delete [] fFailureS;
+ delete [] fFailureP;
+ delete [] fRecons;
+ delete [] fNonRecons;
// delete PlaneEff
delete fPlaneEffSPD;
+ fPlaneEffSPD=0;
+ if(fPlaneEffBkg) {
+ delete fPlaneEffBkg;
+ fPlaneEffBkg=0;
+
+ }
}
//____________________________________________________________________
void
-AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/,AliStack *pStack) {
+AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef, Bool_t lbkg) {
//
- // - calls LoadClusterArray that finds the position of the clusters
- // (in global coord)
- // - convert the cluster coordinates to theta, phi (seen from the
- // interaction vertex). Find the extrapolation/interpolation point.
+ // - you have to take care of the following, before of using Reconstruct
+ // 1) call LoadClusters(TTree* cl) that finds the position of the clusters (in global coord)
+ // and convert the cluster coordinates to theta, phi (seen from the
+ // interaction vertex).
+ // 2) call SetVertex(vtxPos, vtxErr) which set the position of the vertex
+ // - Find the extrapolation/interpolation point.
// - Find the chip corresponding to that
// - Check if there is a cluster near that point
//
-
// reset counters
- fNClustersLay1 = 0;
- fNClustersLay2 = 0;
+ if(lbkg && !GetLightBkgStudyInParallel()) {
+ AliError("You asked for lightBackground in the Reconstruction without proper call to SetLightBkgStudyInParallel(1)");
+ return;
+ }
+ AliITSPlaneEffSPD *pe;
+ if(lbkg) {
+ pe=fPlaneEffBkg;
+ } else {
+ pe=fPlaneEffSPD;
+ }
fNTracklets = 0;
- fNSingleCluster = 0;
- // loading the clusters
- LoadClusterArrays(clusterTree);
- if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
+ // retrieve the vertex position
+ Float_t vtx[3];
+ vtx[0]=(Float_t)GetX();
+ vtx[1]=(Float_t)GetY();
+ vtx[2]=(Float_t)GetZ();
+ // to study residual background (i.e. contribution from TT' to measured efficiency)
+ if(fReflectClusterAroundZAxisForLayer0 && !lbkg) ReflectClusterAroundZAxisForLayer(0);
+ if(fReflectClusterAroundZAxisForLayer1 && !lbkg) ReflectClusterAroundZAxisForLayer(1);
+ //
+ if(fMC && !pStack && !lbkg) {AliError("You asked for MC infos but AliStack not properly loaded"); return;}
+ if(fMC && !tRef && !lbkg) {AliError("You asked for MC infos but TrackRef Tree not properly loaded"); return;}
Bool_t found;
Int_t nfTraPred1=0; Int_t ntTraPred1=0;
Int_t nfTraPred2=0; Int_t ntTraPred2=0;
Int_t nfClu1=0; Int_t ntClu1=0;
Int_t nfClu2=0; Int_t ntClu2=0;
+ // Set fChipUpdatedInEvent=kFALSE for all the chips (none of the chip efficiency already updated
+ // for this new event)
+ for(Int_t i=0;i<1200;i++) fChipUpdatedInEvent[i] = kFALSE;
// find the tracklets
AliDebug(1,"Looking for tracklets... ");
fChipPredOnLay2[iC1] = key;
fAssociationFlag1[iC1] = kFALSE;
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
Float_t eta=fClustersLay1[iC1][0];
eta= TMath::Tan(eta/2.);
eta=-TMath::Log(eta);
- fhetaClustersLay1->Fill(eta);
+ fhetaClustersLay1->Fill(eta);
fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
+ fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
}
}
-
// Loop on layer 2 : finding theta, phi and r
for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
Float_t x = fClustersLay2[iC2][0] - vtx[0];
fChipPredOnLay1[iC2] = key;
fAssociationFlag[iC2] = kFALSE;
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
Float_t eta=fClustersLay2[iC2][0];
eta= TMath::Tan(eta/2.);
eta=-TMath::Log(eta);
fhetaClustersLay2->Fill(eta);
fhphiClustersLay2->Fill(fClustersLay2[iC2][1]);
+ fhClustersInChip->Fill(fhClustersInChip->GetBinCenter(key+1)); // if found=kFALSE -> overflow
}
}
// Loop on layer 1
for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
+ // here the control to check whether the efficiency of the chip traversed by this tracklet
+ // prediction has already been updated in this event using another tracklet prediction
+ if(fUpdateOncePerEventPlaneEff && fChipPredOnLay2[iC1]<1200 && fChipUpdatedInEvent[fChipPredOnLay2[iC1]]) continue;
+
// reset of variables for multiple candidates
Int_t iC2WithBestDist = 0; // reset
Float_t distmin = 100.; // just to put a huge number!
Float_t dZetamin = 0.; // Used for histograms only!
// in any case, if MC has been required, store statistics of primaries and secondaries
- if (fMC) {
+ Bool_t primary=kFALSE; Bool_t secondary=kFALSE; // it is better to have both since chip might not be found
+ if (fMC && !lbkg) {
Int_t lab1=(Int_t)fClustersLay1[iC1][3];
Int_t lab2=(Int_t)fClustersLay1[iC1][4];
Int_t lab3=(Int_t)fClustersLay1[iC1][5];
(lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
{ // this cluster is from a primary particle
fClusterPrimary[key]++;
+ primary=kTRUE;
if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
} else { // this cluster is from a secondary particle
fClusterSecondary[key]++;
+ secondary=kTRUE;
if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
}
}
(lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
(lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay2[iC1]]++;
else fPredictionSecondary[fChipPredOnLay2[iC1]]++;
+ if((lab1 != -2 && IsReconstructableAt(1,iC1,lab1,vtx,pStack,tRef)) ||
+ (lab2 != -2 && IsReconstructableAt(1,iC1,lab2,vtx,pStack,tRef)) ||
+ (lab3 != -2 && IsReconstructableAt(1,iC1,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay2[iC1]]++;
+ else fNonRecons[fChipPredOnLay2[iC1]]++;
}
}
Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2];
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
fhClustersDPhiAll->Fill(dPhi);
fhClustersDThetaAll->Fill(dTheta);
fhClustersDZetaAll->Fill(dZeta);
}
// make "elliptical" cut in Phi and Zeta!
- Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow +
- dZeta*dZeta/fZetaWindow/fZetaWindow);
+ Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
+ dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
if (d>1) continue;
if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
fhClustersDPhiAcc->Fill(dPhimin);
fhClustersDThetaAcc->Fill(dThetamin);
fhClustersDZetaAcc->Fill(dZetamin);
fTracklets[fNTracklets][3] = -2;
}
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
Float_t eta=fTracklets[fNTracklets][0];
eta= TMath::Tan(eta/2.);
eta=-TMath::Log(eta);
nfClu2+=(Int_t)found; // this for debugging purpose
ntClu2++; // to check efficiency of the method FindChip
if(key<1200) { // the Chip has been found
- if(fMC) { // this part only for MC
+ if(fMC && !lbkg) { // this part only for MC
// Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3];
// Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4];
// Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5];
+ if (label2 < 3) {
+ fSuccessTT[key]++;
+ if(primary) fSuccessPP[key]++;
+ }
if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
}
if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked !
- // OK, success
- fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
+ // OK, success
+ pe->UpDatePlaneEff(kTRUE,key); // success
+ fChipUpdatedInEvent[key]=kTRUE;
+ if(fMC && !lbkg) {
+ if(primary) fSuccessP[key]++;
+ if(secondary) fSuccessS[key]++;
+ }
}
else {
- fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
- // (might be in the tracking tollerance)
+ pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
+ fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
+ if(fMC && !lbkg) {
+ if(primary) fSuccessP[key]++;
+ if(secondary) fSuccessS[key]++;
+ }
}
}
fNTracklets++;
} // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction)
- else if (fChipPredOnLay2[iC1]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
-
+ else if (fChipPredOnLay2[iC1]<1200) {
+ pe->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]);
+ fChipUpdatedInEvent[fChipPredOnLay2[iC1]]=kTRUE;
+ if(fMC && !lbkg) {
+ if(primary) fFailureP[fChipPredOnLay2[iC1]]++;
+ if(secondary) fFailureS[fChipPredOnLay2[iC1]]++;
+ }
+ }
} // end of loop over clusters in layer 1
fNTracklets1=fNTracklets;
// Loop on layer 2
for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
+ // here the control to check whether the efficiency of the chip traversed by this tracklet
+ // prediction has already been updated in this event using another tracklet prediction
+ if(fUpdateOncePerEventPlaneEff && fChipPredOnLay1[iC2]<1200 && fChipUpdatedInEvent[fChipPredOnLay1[iC2]]) continue;
+
// reset of variables for multiple candidates
Int_t iC1WithBestDist = 0; // reset
Float_t distmin = 100.; // just to put a huge number!
Float_t dZetamin = 0.; // Used for histograms only!
// in any case, if MC has been required, store statistics of primaries and secondaries
- if (fMC) {
+ Bool_t primary=kFALSE; Bool_t secondary=kFALSE;
+ if (fMC && !lbkg) {
Int_t lab1=(Int_t)fClustersLay2[iC2][3];
Int_t lab2=(Int_t)fClustersLay2[iC2][4];
Int_t lab3=(Int_t)fClustersLay2[iC2][5];
(lab3 != -2 && PrimaryTrackChecker(lab3,pStack)))
{ // this cluster is from a primary particle
fClusterPrimary[key]++;
+ primary=kTRUE;
if(fUseOnlySecondaryForPred) continue; // skip this tracklet built with a primary track
} else { // this cluster is from a secondary particle
fClusterSecondary[key]++;
+ secondary=kTRUE;
if(fUseOnlyPrimaryForPred) continue; // skip this tracklet built with a secondary track
}
}
(lab2 != -2 && PrimaryTrackChecker(lab2,pStack) ) ||
(lab3 != -2 && PrimaryTrackChecker(lab3,pStack))) fPredictionPrimary[fChipPredOnLay1[iC2]]++;
else fPredictionSecondary[fChipPredOnLay1[iC2]]++;
+ if((lab1 != -2 && IsReconstructableAt(0,iC2,lab1,vtx,pStack,tRef)) ||
+ (lab2 != -2 && IsReconstructableAt(0,iC2,lab2,vtx,pStack,tRef)) ||
+ (lab3 != -2 && IsReconstructableAt(0,iC2,lab3,vtx,pStack,tRef))) fRecons[fChipPredOnLay1[iC2]]++;
+ else fNonRecons[fChipPredOnLay1[iC2]]++;
}
}
Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2];
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
fhClustersDPhiInterpAll->Fill(dPhi);
fhClustersDThetaInterpAll->Fill(dTheta);
fhClustersDZetaInterpAll->Fill(dZeta);
}
} // end of loop over clusters in layer 1
- if (distmin<100) { // This means that a cluster in layer 1 was found that mathes with iC2
+ if (distmin<100) { // This means that a cluster in layer 1 was found that matches with iC2
- if (fHistOn) {
+ if (fHistOn && !lbkg) {
fhClustersDPhiInterpAcc->Fill(dPhimin);
fhClustersDThetaInterpAcc->Fill(dThetamin);
fhClustersDZetaInterpAcc->Fill(dZetamin);
nfClu1+=(Int_t)found; // this for debugging purpose
ntClu1++; // to check efficiency of the method FindChip
if(key<1200) {
- if(fMC) { // this part only for MC
+ if(fMC && !lbkg) { // this part only for MC
// Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3];
// Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4];
// Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5];
+ if (label2 < 3) { // same label
+ fSuccessTT[key]++;
+ if(primary) fSuccessPP[key]++;
+ }
if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it)
if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it)
}
if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked !
- // OK, success
- fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success
+ // OK, success
+ pe->UpDatePlaneEff(kTRUE,key); // success
+ fChipUpdatedInEvent[key]=kTRUE;
+ if(fMC && !lbkg) {
+ if(primary) fSuccessP[key]++;
+ if(secondary) fSuccessS[key]++;
+ }
} else {
- fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure
- // (might be in the tracking tollerance)
+ pe->UpDatePlaneEff(kTRUE,key); // this should not be a failure
+ fChipUpdatedInEvent[key]=kTRUE; // (might be in the tracking tollerance)
+ if(fMC && !lbkg) {
+ if(primary) fSuccessP[key]++;
+ if(secondary) fSuccessS[key]++;
+ }
}
}
fNTracklets++;
} // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction)
- else if (fChipPredOnLay1[iC2]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
-
+ else if (fChipPredOnLay1[iC2]<1200) {
+ pe->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]);
+ fChipUpdatedInEvent[fChipPredOnLay1[iC2]]=kTRUE;
+ if(fMC && !lbkg) {
+ if(primary) fFailureP[fChipPredOnLay1[iC2]]++;
+ if(secondary) fFailureS[fChipPredOnLay1[iC2]]++;
+ }
+ }
} // end of loop over clusters in layer 2
AliDebug(1,Form("%d tracklets found", fNTracklets));
AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2));
}
//____________________________________________________________________
-Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx,
+Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer,const Float_t* vtx,
Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) {
//
// Input: a) layer number in the range [0,1]
Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference
// of the intersection of the tracklet with the pixel layer.
if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction
- else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
+ else {
+ if(TMath::Abs(thetaVtx)<1E-6) return kFALSE;
+ zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction
+ }
AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2]));
Double_t vtxy[2]={vtx[0],vtx[1]};
if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices
}
//______________________________________________________________________________
Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) {
+//
+// Return the average radius of a layer from Geometry
+//
if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;}
Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
//______________________________________________________________________________
Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) {
//--------------------------------------------------------------------
- //This function finds the detector crossed by the track
+ // This function finds the detector crossed by the track
+ // Input: layer in range [0,1]
+ // phi in ALICE absolute reference system
+ // z " " " " "
//--------------------------------------------------------------------
if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;}
Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 !
r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle
In the same system, the equation of a semi-line is: phi=phiVtx;
Hence you get one interception only: P=(r,phiVtx)
-Finally you want P in the ABSOLUTE ALICE system.
+Finally you want P in the ABSOLUTE ALICE reference system.
*/
Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin
Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin
return kTRUE;
}
//___________________________________________________________
-Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) {
+Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) const {
+//
+// simple method to reduce all angles (in rad)
+// in range [0,2pi[
+//
+//
while(angle >=2*TMath::Pi() || angle<0) {
if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi();
if(angle < 0) angle+=2*TMath::Pi();
}
//___________________________________________________________
Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) {
+//
+// This method check if a particle is primary; i.e.
+// it comes from the main vertex and it is a "stable" particle, according to
+// AliStack::IsPhysicalPrimary() (note that there also Sigma0 are considered as
+// a stable particle: it has no effect on this analysis).
+// This method can be called only for MC events, where Kinematics is available.
+// if fUseOnlyStableParticle is kTRUE (via SetUseOnlyStableParticle) then it
+// returns kTRUE if also AliITSTrackleterSPDEff::DecayingTrackChecker() return 0.
+// The latter (see below) try to verify if a primary particle is also "detectable".
+//
if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;}
if(!stack) {AliError("null pointer to MC stack"); return kFALSE;}
if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;}
part0->Vx(),part0->Vy(),part->Vx(),part->Vy()));
return kFALSE; }// primary if within 500 microns from true Vertex
- if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)<2) return kFALSE;
+ if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)>0) return kFALSE;
return kTRUE;
}
//_____________________________________________________________________________________________
Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) {
+//
+// This private method can be applied on MC particles (if stack is available),
+// provided they have been identified as "primary" from PrimaryTrackChecker() (see above).
+//
+// It define "detectable" a primary particle according to the following criteria:
+//
+// - if no decay products can be found in the stack (note that this does not
+// means it is stable, since a particle is stored in stack if it has at least 1 hit in a
+// sensitive detector)
+// - if it has at least one decay daughter produced outside or just on the outer pixel layer
+// - if the last decay particle is an electron (or a muon) which is not produced in-between
+// the two pixel layers (this is likely to be a kink).
if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;}
if(!stack) {AliError("null pointer to MC stack"); return 0;}
if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;}
Int_t nret=0;
TParticle* dau = 0;
Int_t nDau = 0;
- Int_t firstDau = part->GetFirstDaughter();
- if (firstDau > 0) {
+ Int_t pdgDau;
+ Int_t firstDau = part->GetFirstDaughter(); // if no daugther stored then no way to understand i
+ // its real fate ! But you have to take it !
+ if (firstDau > 0) { // if it has daugther(s) try to infer if it is "detectable" as a tracklet
Int_t lastDau = part->GetLastDaughter();
nDau = lastDau - firstDau + 1;
- //printf("number of daugthers %d \n",nDau);
- if (nDau > 0) {
- //for(Int_t j=firstDau; j<=lastDau; j++)
- for(Int_t j=firstDau; j<=firstDau; j++)
- { // only first one
- dau = stack->Particle(j);
- Double_t distx = dau->Vx()-part->Vx();
- Double_t disty = dau->Vy()-part->Vy();
- Double_t distz = dau->Vz()-part->Vz();
- Double_t distR = TMath::Sqrt(distx*distx+disty*disty+distz*distz);
- if (distR > GetRLayer(0)+0.5) nret=1; // decay after first pixel layer
- if (distR > GetRLayer(1)+0.5) nret=2; // decay after second pixel layer
- }
+ Double_t distMax=0.;
+ Int_t jmax=0;
+ for(Int_t j=firstDau; j<=lastDau; j++) {
+ dau = stack->Particle(j);
+ Double_t distx = dau->Vx();
+ Double_t disty = dau->Vy();
+ //Double_t distz = dau->Vz();
+ Double_t distR = TMath::Sqrt(distx*distx+disty*disty);
+ if(distR<distMax) continue; // considere only the daughter produced at largest radius
+ distMax=distR;
+ jmax=j;
}
- } else nret = 3; // stable particle
-return nret;
+ dau = stack->Particle(jmax);
+ pdgDau=dau->GetPdgCode();
+ if (pdgDau == 11 || pdgDau == 13 ) {
+ if(distMax < GetRLayer(1)-0.25 && distMax > GetRLayer(0)+0.27) nret=1; // can be a kink (reject it)
+ else nret =0; // delta-ray emission in material (keep it)
+ }
+ else {// not ele or muon
+ if (distMax < GetRLayer(1)-0.25 ) nret= 1;} // decay before the second pixel layer (reject it)
+ }
+return nret;
}
//_________________________________________________________________
void AliITSTrackleterSPDEff::InitPredictionMC() {
+//
+// this method allocate memory for the MC related informations
+// all the counters are set to 0
+//
+//
if(!fMC) {AliError("This method works only if SetMC() has been called"); return;}
fPredictionPrimary = new Int_t[1200];
fPredictionSecondary = new Int_t[1200];
fClusterPrimary = new Int_t[1200];
fClusterSecondary = new Int_t[1200];
+fSuccessPP = new Int_t[1200];
+fSuccessTT = new Int_t[1200];
+fSuccessS = new Int_t[1200];
+fSuccessP = new Int_t[1200];
+fFailureS = new Int_t[1200];
+fFailureP = new Int_t[1200];
+fRecons = new Int_t[1200];
+fNonRecons = new Int_t[1200];
for(Int_t i=0; i<1200; i++) {
fPredictionPrimary[i]=0;
fPredictionSecondary[i]=0;
fPredictionSecondary[i]=0;
fClusterSecondary[i]=0;
+ fSuccessPP[i]=0;
+ fSuccessTT[i]=0;
+ fSuccessS[i]=0;
+ fSuccessP[i]=0;
+ fFailureS[i]=0;
+ fFailureP[i]=0;
+ fRecons[i]=0;
+ fNonRecons[i]=0;
+}
+return;
+}
+//_________________________________________________________________
+void AliITSTrackleterSPDEff::DeletePredictionMC() {
+//
+// this method deallocate memory for the MC related informations
+// all the counters are set to 0
+//
+//
+if(fMC) {AliInfo("This method works only if fMC=kTRUE"); return;}
+if(fPredictionPrimary) {
+ delete fPredictionPrimary; fPredictionPrimary=0;
+}
+if(fPredictionSecondary) {
+ delete fPredictionSecondary; fPredictionSecondary=0;
+}
+if(fClusterPrimary) {
+ delete fClusterPrimary; fClusterPrimary=0;
+}
+if(fClusterSecondary) {
+ delete fClusterSecondary; fClusterSecondary=0;
+}
+if(fSuccessPP) {
+ delete fSuccessPP; fSuccessPP=0;
+}
+if(fSuccessTT) {
+ delete fSuccessTT; fSuccessTT=0;
+}
+if(fSuccessS) {
+ delete fSuccessS; fSuccessS=0;
+}
+if(fSuccessP) {
+ delete fSuccessP; fSuccessP=0;
+}
+if(fFailureS) {
+ delete fFailureS; fFailureS=0;
+}
+if(fFailureP) {
+ delete fFailureP; fFailureP=0;
+}
+if(fRecons) {
+ delete fRecons; fRecons=0;
+}
+if(fNonRecons) {
+ delete fNonRecons; fNonRecons=0;
}
return;
}
//______________________________________________________________________
Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const {
+//
+// This method return the Data menmber fPredictionPrimary [1200].
+// You can call it only for MC events.
+// fPredictionPrimary[key] contains the number of tracklet predictions on the
+// given chip key built using a cluster on the other layer produced (at least)
+// from a primary particle.
+// Key refers to the chip crossed by the prediction
+//
+//
if (!fMC) {CallWarningMC(); return 0;}
if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
return fPredictionPrimary[(Int_t)key];
}
//______________________________________________________________________
Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const {
+//
+// This method return the Data menmber fPredictionSecondary [1200].
+// You can call it only for MC events.
+// fPredictionSecondary[key] contains the number of tracklet predictions on the
+// given chip key built using a cluster on the other layer produced (only)
+// from a secondary particle
+// Key refers to the chip crossed by the prediction
+//
+//
if (!fMC) {CallWarningMC(); return 0;}
if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
return fPredictionSecondary[(Int_t)key];
}
//______________________________________________________________________
Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const {
+//
+// This method return the Data menmber fClusterPrimary [1200].
+// You can call it only for MC events.
+// fClusterPrimary[key] contains the number of tracklet predictions
+// built using a cluster on that layer produced (only)
+// from a primary particle
+// Key refers to the chip used to build the prediction
+//
+//
if (!fMC) {CallWarningMC(); return 0;}
if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
return fClusterPrimary[(Int_t)key];
}
//______________________________________________________________________
Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const {
+//
+// This method return the Data menmber fClusterSecondary [1200].
+// You can call it only for MC events.
+// fClusterSecondary[key] contains the number of tracklet predictions
+// built using a cluster on that layer produced (only)
+// from a secondary particle
+// Key refers to the chip used to build the prediction
+//
if (!fMC) {CallWarningMC(); return 0;}
if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
return fClusterSecondary[(Int_t)key];
}
//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetSuccessPP(const UInt_t key) const {
+//
+// This method return the Data menmber fSuccessPP [1200].
+// You can call it only for MC events.
+// fSuccessPP[key] contains the number of successes (i.e. a tracklet prediction matching
+// with a cluster on the other layer) built by using the same primary particle
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fSuccessPP[(Int_t)key];
+}
+//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetSuccessTT(const UInt_t key) const {
+//
+// This method return the Data menmber fSuccessTT [1200].
+// You can call it only for MC events.
+// fSuccessTT[key] contains the number of successes (i.e. a tracklet prediction matching
+// with a cluster on the other layer) built by using the same particle (whatever)
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fSuccessTT[(Int_t)key];
+}
+//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetSuccessS(const UInt_t key) const {
+//
+// This method return the Data menmber fSuccessS [1200].
+// You can call it only for MC events.
+// fSuccessS[key] contains the number of successes (i.e. a tracklet prediction matching
+// with a cluster on the other layer) built by using a secondary particle
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fSuccessS[(Int_t)key];
+}
+//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetSuccessP(const UInt_t key) const {
+//
+// This method return the Data menmber fSuccessP [1200].
+// You can call it only for MC events.
+// fSuccessP[key] contains the number of successes (i.e. a tracklet prediction matching
+// with a cluster on the other layer) built by using a primary particle
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fSuccessP[(Int_t)key];
+}
+//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetFailureS(const UInt_t key) const {
+//
+// This method return the Data menmber fFailureS [1200].
+// You can call it only for MC events.
+// fFailureS[key] contains the number of failures (i.e. a tracklet prediction not matching
+// with a cluster on the other layer) built by using a secondary particle
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fFailureS[(Int_t)key];
+}
+//______________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetFailureP(const UInt_t key) const {
+//
+// This method return the Data menmber fFailureP [1200].
+// You can call it only for MC events.
+// fFailureP[key] contains the number of failures (i.e. a tracklet prediction not matching
+// with a cluster on the other layer) built by using a primary particle
+// the unique chip key refers to the chip which get updated its efficiency
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fFailureP[(Int_t)key];
+}
+//_____________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetRecons(const UInt_t key) const {
+//
+// This method return the Data menmber fRecons [1200].
+// You can call it only for MC events.
+// fRecons[key] contains the number of reconstractable tracklets (i.e. a tracklet prediction which
+// has an hit in the detector)
+// the unique chip key refers to the chip where fall the prediction
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fRecons[(Int_t)key];
+}
+//_____________________________________________________________________
+Int_t AliITSTrackleterSPDEff::GetNonRecons(const UInt_t key) const {
+//
+// This method return the Data menmber fNonRecons [1200].
+// You can call it only for MC events.
+// fRecons[key] contains the number of unreconstractable tracklets (i.e. a tracklet prediction which
+// has not any hit in the detector)
+// the unique chip key refers to the chip where fall the prediction
+//
+if (!fMC) {CallWarningMC(); return 0;}
+if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;}
+return fNonRecons[(Int_t)key];
+}
+//______________________________________________________________________
void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{
// Print out some class data values in Ascii Form to output stream
// Inputs:
// none.
// Return:
// none.
- *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow ;
+ *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2
+ << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2
+ << " " << fUpdateOncePerEventPlaneEff << " " << fMinContVtx
+ << " " << fReflectClusterAroundZAxisForLayer0
+ << " " << fReflectClusterAroundZAxisForLayer1;
*os << " " << fMC;
if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;}
*os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred
for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ;
for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ;
for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessPP(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessTT(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessS(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetSuccessP(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureS(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetFailureP(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetRecons(i) ;
+ for(Int_t i=0;i<1200;i++) *os <<" "<< GetNonRecons(i) ;
return;
}
//______________________________________________________________________
// Return:
// none.
- *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow;
+ Bool_t tmp= fMC;
+ *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2
+ >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2
+ >> fUpdateOncePerEventPlaneEff >> fMinContVtx
+ >> fReflectClusterAroundZAxisForLayer0
+ >> fReflectClusterAroundZAxisForLayer1;
+ //if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
*is >> fMC;
- if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;}
- *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
- >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
- >> fUseOnlyStableParticle;
- for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
- for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
- for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
- for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
+ if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
+ else {
+ if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
+ *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred
+ >> fUseOnlySameParticle >> fUseOnlyDifferentParticle
+ >> fUseOnlyStableParticle;
+ for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fSuccessPP[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fSuccessTT[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fSuccessS[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fSuccessP[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fFailureS[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fFailureP[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fRecons[i] ;
+ for(Int_t i=0;i<1200;i++) *is >> fNonRecons[i] ;
+ }
return;
}
//______________________________________________________________________
}
//______________________________________________________________________
void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const {
- if(!fMC) {CallWarningMC(); return;}
- ofstream out(filename.Data(),ios::out | ios::binary);
- out << *this;
- out.close();
+//
+// This Method write into an either asci or root file
+// the used cuts and the statistics of the MC related quantities
+// The method SetMC() has to be called before
+// Input TString filename: name of file for output (it deletes already existing
+// file)
+// Output: none
+//
+//
+ //if(!fMC) {CallWarningMC(); return;}
+ if (!filename.Contains(".root")) {
+ ofstream out(filename.Data(),ios::out | ios::binary);
+ out << *this;
+ out.close();
+ return;
+ }
+ else {
+ TFile* mcfile = TFile::Open(filename, "RECREATE");
+ TH1F* cuts = new TH1F("cuts", "list of cuts", 11, 0, 11); // TH1I containing cuts
+ cuts->SetBinContent(1,fPhiWindowL1);
+ cuts->SetBinContent(2,fZetaWindowL1);
+ cuts->SetBinContent(3,fPhiWindowL2);
+ cuts->SetBinContent(4,fZetaWindowL2);
+ cuts->SetBinContent(5,fOnlyOneTrackletPerC1);
+ cuts->SetBinContent(6,fOnlyOneTrackletPerC2);
+ cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff);
+ cuts->SetBinContent(8,fMinContVtx);
+ cuts->SetBinContent(9,fReflectClusterAroundZAxisForLayer0);
+ cuts->SetBinContent(10,fReflectClusterAroundZAxisForLayer1);
+ cuts->SetBinContent(11,fMC);
+ cuts->Write();
+ delete cuts;
+ if(!fMC) {AliInfo("Writing only cuts, no MC info");}
+ else {
+ TH1C* mc0 = new TH1C("mc0", "mc cuts", 5, 0, 5);
+ mc0->SetBinContent(1,fUseOnlyPrimaryForPred);
+ mc0->SetBinContent(2,fUseOnlySecondaryForPred);
+ mc0->SetBinContent(3,fUseOnlySameParticle);
+ mc0->SetBinContent(4,fUseOnlyDifferentParticle);
+ mc0->SetBinContent(5,fUseOnlyStableParticle);
+ mc0->Write();
+ delete mc0;
+ TH1I *mc1;
+ mc1 = new TH1I("mc1", "mc info PredictionPrimary", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionPrimary(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc2", "mc info PredictionSecondary", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetPredictionSecondary(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc3", "mc info ClusterPrimary", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterPrimary(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc4", "mc info ClusterSecondary", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetClusterSecondary(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc5", "mc info SuccessPP", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessPP(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc6", "mc info SuccessTT", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessTT(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc7", "mc info SuccessS", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessS(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc8", "mc info SuccessP", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetSuccessP(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc9", "mc info FailureS", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureS(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc10", "mc info FailureP", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetFailureP(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc11", "mc info Recons", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetRecons(i)) ;
+ mc1->Write();
+ mc1 = new TH1I("mc12", "mc info NonRecons", 1200, 0, 1200);
+ for(Int_t i=0;i<1200;i++) mc1->SetBinContent(i+1,GetNonRecons(i)) ;
+ mc1->Write();
+ delete mc1;
+ }
+ mcfile->Close();
+ }
return;
}
//____________________________________________________________________
void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) {
- if(!fMC) {CallWarningMC(); return;}
+//
+// This Method read from an asci file (do not know why binary does not work)
+// the cuts to be used and the statistics of the MC related quantities
+// Input TString filename: name of input file for output
+// The method SetMC() has to be called before
+// Output: none
+//
+//
+ //if(!fMC) {CallWarningMC(); return;}
if( gSystem->AccessPathName( filename.Data() ) ) {
AliError( Form( "file (%s) not found", filename.Data() ) );
return;
}
- ifstream in(filename.Data(),ios::in | ios::binary);
- in >> *this;
- in.close();
+ if (!filename.Contains(".root")) {
+ ifstream in(filename.Data(),ios::in | ios::binary);
+ in >> *this;
+ in.close();
+ return;
+ }
+ else {
+ Bool_t tmp= fMC;
+ TFile *mcfile = TFile::Open(filename);
+ TH1F *cuts = (TH1F*)mcfile->Get("cuts");
+ fPhiWindowL1=(Float_t)cuts->GetBinContent(1);
+ fZetaWindowL1=(Float_t)cuts->GetBinContent(2);
+ fPhiWindowL2=(Float_t)cuts->GetBinContent(3);
+ fZetaWindowL2=(Float_t)cuts->GetBinContent(4);
+ fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5);
+ fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6);
+ fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7);
+ fMinContVtx=(Int_t)cuts->GetBinContent(8);
+ fReflectClusterAroundZAxisForLayer0=(Bool_t)cuts->GetBinContent(9);
+ fReflectClusterAroundZAxisForLayer1=(Bool_t)cuts->GetBinContent(10);
+ fMC=(Bool_t)cuts->GetBinContent(11);
+ if(!fMC) {AliInfo("Reading only cuts, no MC info"); if(tmp) SetMC(kFALSE); }
+ else { // only if file with MC predictions
+ if(!tmp) {AliInfo("Calling SetMC() to read this file wtih MC info"); SetMC();}
+ TH1C *mc0 = (TH1C*)mcfile->Get("mc0");
+ fUseOnlyPrimaryForPred=(Bool_t)mc0->GetBinContent(1);
+ fUseOnlySecondaryForPred=(Bool_t)mc0->GetBinContent(2);
+ fUseOnlySameParticle=(Bool_t)mc0->GetBinContent(3);
+ fUseOnlyDifferentParticle=(Bool_t)mc0->GetBinContent(4);
+ fUseOnlyStableParticle=(Bool_t)mc0->GetBinContent(5);
+ TH1I *mc1;
+ mc1 =(TH1I*)mcfile->Get("mc1");
+ for(Int_t i=0;i<1200;i++) fPredictionPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc2");
+ for(Int_t i=0;i<1200;i++) fPredictionSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc3");
+ for(Int_t i=0;i<1200;i++) fClusterPrimary[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc4");
+ for(Int_t i=0;i<1200;i++) fClusterSecondary[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc5");
+ for(Int_t i=0;i<1200;i++) fSuccessPP[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc6");
+ for(Int_t i=0;i<1200;i++) fSuccessTT[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc7");
+ for(Int_t i=0;i<1200;i++) fSuccessS[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc8");
+ for(Int_t i=0;i<1200;i++) fSuccessP[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc9");
+ for(Int_t i=0;i<1200;i++) fFailureS[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc10");
+ for(Int_t i=0;i<1200;i++) fFailureP[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc11");
+ for(Int_t i=0;i<1200;i++) fRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ mc1 =(TH1I*)mcfile->Get("mc12");
+ for(Int_t i=0;i<1200;i++) fNonRecons[i]=(Int_t)mc1->GetBinContent(i+1) ;
+ }
+ mcfile->Close();
+ }
return;
}
//____________________________________________________________________
Bool_t AliITSTrackleterSPDEff::SaveHists() {
- // This method save the histograms on the output file
+ // This (private) method save the histograms on the output file
// (only if fHistOn is TRUE).
+ // Also the histograms from the base class are saved through the
+ // AliITSMultReconstructor::SaveHists() call
if (!GetHistOn()) return kFALSE;
- AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
+// AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class
+ fhClustersDPhiAll->Write();
+ fhClustersDThetaAll->Write();
+ fhClustersDZetaAll->Write();
+ fhDPhiVsDThetaAll->Write();
+ fhDPhiVsDZetaAll->Write();
+
+ fhClustersDPhiAcc->Write();
+ fhClustersDThetaAcc->Write();
+ fhClustersDZetaAcc->Write();
+ fhDPhiVsDThetaAcc->Write();
+ fhDPhiVsDZetaAcc->Write();
+
+ fhetaTracklets->Write();
+ fhphiTracklets->Write();
+ fhetaClustersLay1->Write();
+ fhphiClustersLay1->Write();
fhClustersDPhiInterpAll->Write();
fhClustersDThetaInterpAll->Write();
fhetaClustersLay2->Write();
fhphiClustersLay2->Write();
+ fhClustersInChip->Write();
+ for (Int_t nhist=0;nhist<80;nhist++){
+ fhClustersInModuleLay1[nhist]->Write();
+ }
+ for (Int_t nhist=0;nhist<160;nhist++){
+ fhClustersInModuleLay2[nhist]->Write();
+ }
return kTRUE;
}
//__________________________________________________________
Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) {
//
// Saves the histograms into a tree and saves the trees into a file
+ // Also the histograms from the base class are saved
//
if (!GetHistOn()) return kFALSE;
- if (filename.Data()=="") {
+ if (!strcmp(filename.Data(),"")) {
AliWarning("WriteHistosToFile: null output filename!");
return kFALSE;
}
}
//____________________________________________________________
void AliITSTrackleterSPDEff::BookHistos() {
+//
+// This method books addtitional histograms
+// w.r.t. those of the base class.
+// In particular, the differences of cluster coordinate between the two SPD
+// layers are computed in the interpolation phase
+//
if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;}
+//
+ fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
+ fhClustersDPhiAcc->SetDirectory(0);
+ fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
+ fhClustersDThetaAcc->SetDirectory(0);
+ fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
+ fhClustersDZetaAcc->SetDirectory(0);
+
+ fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
+ fhDPhiVsDZetaAcc->SetDirectory(0);
+ fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
+ fhDPhiVsDThetaAcc->SetDirectory(0);
+
+ fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
+ fhClustersDPhiAll->SetDirectory(0);
+ fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
+ fhClustersDThetaAll->SetDirectory(0);
+ fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
+ fhClustersDZetaAll->SetDirectory(0);
+
+ fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
+ fhDPhiVsDZetaAll->SetDirectory(0);
+ fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
+ fhDPhiVsDThetaAll->SetDirectory(0);
+
+ fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
+ fhetaTracklets->SetDirectory(0);
+ fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
+ fhphiTracklets->SetDirectory(0);
+ fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
+ fhetaClustersLay1->SetDirectory(0);
+ fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
+ fhphiClustersLay1->SetDirectory(0);
+//
fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1);
fhClustersDPhiInterpAcc->SetDirectory(0);
fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1);
fhetaClustersLay2->SetDirectory(0);
fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi());
fhphiClustersLay2->SetDirectory(0);
+ fhClustersInChip = new TH1F("fhClustersInChip", "ClustersPerChip", 1200, -0.5, 1199.5);
+ fhClustersInChip->SetDirectory(0);
+// each chip is divided 8(z) x 4(y), i.e. in 32 squares, each containing 4 columns and 64 rows.
+ Float_t bz[160]; const Float_t kconv = 1.0E-04; // converts microns to cm.
+ for(Int_t i=0;i<160;i++) bz[i] = 425.0; // most are 425 microns except below
+ bz[ 31] = bz[ 32] = 625.0; // first chip boundry
+ bz[ 63] = bz[ 64] = 625.0; // first chip boundry
+ bz[ 95] = bz[ 96] = 625.0; // first chip boundry
+ bz[127] = bz[128] = 625.0; // first chip boundry
+ Double_t xbins[41]; // each bin in x (Z loc coordinate) includes 4 columns
+ //xbins[0]=0;
+ Float_t xmn,xmx,zmn,zmx;
+ if(!fPlaneEffSPD->GetBlockBoundaries(0,xmn,xmx,zmn,zmx)) AliWarning("Could not book histo properly");
+ xbins[0]=(Double_t)zmn;
+ for(Int_t i=0;i<40;i++) {
+ xbins[i+1]=xbins[i] + (bz[4*i]+bz[4*i+1]+bz[4*i+2]+bz[4*i+3])*kconv;
+ }
+ TString histname="ClustersLay1_mod_",aux;
+ fhClustersInModuleLay1 =new TH2F*[80];
+ for (Int_t nhist=0;nhist<80;nhist++){
+ aux=histname;
+ aux+=nhist;
+ //
+ fhClustersInModuleLay1[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
+ fhClustersInModuleLay1[nhist]->SetName(aux.Data());
+ fhClustersInModuleLay1[nhist]->SetTitle(aux.Data());
+ fhClustersInModuleLay1[nhist]->SetDirectory(0);
+ }
+ histname="ClustersLay2_mod_";
+ fhClustersInModuleLay2 =new TH2F*[160];
+ for (Int_t nhist=0;nhist<160;nhist++){
+ aux=histname;
+ aux+=nhist;
+ fhClustersInModuleLay2[nhist]=new TH2F("histname","histname",40,xbins,4,(Double_t)xmn,(Double_t)xmx);
+ fhClustersInModuleLay2[nhist]->SetName(aux.Data());
+ fhClustersInModuleLay2[nhist]->SetTitle(aux.Data());
+ fhClustersInModuleLay2[nhist]->SetDirectory(0);
+ }
+//
return;
}
//____________________________________________________________
void AliITSTrackleterSPDEff::DeleteHistos() {
+//
+// Private method to delete Histograms from memory
+// it is called. e.g., by the destructor.
+//
+// form AliITSMultReconstructor
+ if(fhClustersDPhiAcc) {delete fhClustersDPhiAcc; fhClustersDPhiAcc=0;}
+ if(fhClustersDThetaAcc) {delete fhClustersDThetaAcc; fhClustersDThetaAcc=0;}
+ if(fhClustersDZetaAcc) {delete fhClustersDZetaAcc; fhClustersDZetaAcc=0;}
+ if(fhClustersDPhiAll) {delete fhClustersDPhiAll; fhClustersDPhiAll=0;}
+ if(fhClustersDThetaAll) {delete fhClustersDThetaAll; fhClustersDThetaAll=0;}
+ if(fhClustersDZetaAll) {delete fhClustersDZetaAll; fhClustersDZetaAll=0;}
+ if(fhDPhiVsDThetaAll) {delete fhDPhiVsDThetaAll; fhDPhiVsDThetaAll=0;}
+ if(fhDPhiVsDThetaAcc) {delete fhDPhiVsDThetaAcc; fhDPhiVsDThetaAcc=0;}
+ if(fhDPhiVsDZetaAll) {delete fhDPhiVsDZetaAll; fhDPhiVsDZetaAll=0;}
+ if(fhDPhiVsDZetaAcc) {delete fhDPhiVsDZetaAcc; fhDPhiVsDZetaAcc=0;}
+ if(fhetaTracklets) {delete fhetaTracklets; fhetaTracklets=0;}
+ if(fhphiTracklets) {delete fhphiTracklets; fhphiTracklets=0;}
+ if(fhetaClustersLay1) {delete fhetaClustersLay1; fhetaClustersLay1=0;}
+ if(fhphiClustersLay1) {delete fhphiClustersLay1; fhphiClustersLay1=0;}
+//
if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;}
if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;}
if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;}
if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;}
if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;}
if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;}
+ if(fhClustersInChip) {delete fhClustersInChip; fhClustersInChip=0;}
+ if(fhClustersInModuleLay1) {
+ for (Int_t i=0; i<80; i++ ) delete fhClustersInModuleLay1[i];
+ delete [] fhClustersInModuleLay1; fhClustersInModuleLay1=0;
+ }
+ if(fhClustersInModuleLay2) {
+ for (Int_t i=0; i<160; i++ ) delete fhClustersInModuleLay2[i];
+ delete [] fhClustersInModuleLay2; fhClustersInModuleLay2=0;
+ }
}
//_______________________________________________________________
+Bool_t AliITSTrackleterSPDEff::IsReconstructableAt(Int_t layer,Int_t iC,Int_t ipart,
+ const Float_t* vtx, const AliStack *stack, TTree *ref) {
+// This (private) method can be used only for MC events, where both AliStack and the TrackReference
+// are available.
+// It is used to asses whether a tracklet prediction is reconstructable or not at the other layer
+// Input:
+// - Int_t layer (either 0 or 1): layer which you want to check if the tracklete can be
+// reconstructed at
+// - Int_t iC : cluster index used to build the tracklet prediction
+// if layer=0 ==> iC=iC2 ; elseif layer=1 ==> iC=iC1
+// - Float_t* vtx: actual event vertex
+// - stack: pointer to Stack
+// - ref: pointer to TTRee of TrackReference
+Bool_t ret=kFALSE; // returned value
+Float_t trefLayExtr[3]; // equivalent to fClustersLay1/fClustersLay2 but for the track reference
+if(!fMC) {AliError("This method works only if SetMC() has been called"); return ret;}
+if(!stack) {AliError("null pointer to MC stack"); return ret;}
+if(!ref) {AliError("null pointer to TrackReference Tree"); return ret;}
+if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return ret;}
+if(layer<0 || layer>1) {AliError("You can extrapolate either at lay 0 or at lay 1"); return ret;}
+
+AliTrackReference *tref=0x0;
+Int_t imatch=-100; // index of the track in TrackReference which matches with ipart
+Int_t nentries = (Int_t)ref->GetEntries();
+TClonesArray *tcaRef = new TClonesArray("AliTrackReference");
+TBranch *br = ref->GetBranch("TrackReferences");
+br->SetAddress(&tcaRef);
+for(Int_t itrack=0;itrack<nentries;itrack++) { // loop over all Tracks in TrackReference to match the ipart one
+ br->GetEntry(itrack);
+ Int_t nref=tcaRef->GetEntriesFast();
+ if(nref>0) { //it is enough to look at the first one
+ tref=(AliTrackReference*)tcaRef->At(0); // it is enough to look at the first one
+ if(tref->GetTrack()==ipart) {imatch=itrack; break;}
+ }
+}
+if(imatch<0) {AliWarning(Form("Could not find AliTrackReference for particle %d",ipart)); return kFALSE;}
+br->GetEntry(imatch); // redundant, nevertheless ...
+Int_t nref=tcaRef->GetEntriesFast();
+for(Int_t iref=0;iref<nref;iref++) { // loop over all the refs of the matching track
+ tref=(AliTrackReference*)tcaRef->At(iref);
+ if(tref->R()>10) continue; // not SPD ref
+ if(layer==0 && tref->R()>5) continue; // ref on SPD outer layer
+ if(layer==1 && tref->R()<5) continue; // ref on SPD inner layer
+
+// compute the proper quantities for this tref, as was done for fClustersLay1/2
+ Float_t x = tref->X() - vtx[0];
+ Float_t y = tref->Y() - vtx[1];
+ Float_t z = tref->Z() - vtx[2];
+
+ Float_t r = TMath::Sqrt(x*x + y*y +z*z);
+
+ trefLayExtr[0] = TMath::ACos(z/r); // Store Theta
+ trefLayExtr[1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
+ trefLayExtr[2] = z; // Store z
+
+ if(layer==1) { // try to see if it is reconstructable at the outer layer
+// find the difference in angles
+ Float_t dPhi = TMath::Abs(trefLayExtr[1] - fClustersLay1[iC][1]);
+ // take into account boundary condition
+ if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
+
+ // find the difference in z (between linear projection from layer 1
+ // and the actual point: Dzeta= z1/r1*r2 -z2)
+ Float_t r2 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
+ Float_t dZeta = TMath::Cos(fClustersLay1[iC][0])*r2 - trefLayExtr[2];
+
+ // make "elliptical" cut in Phi and Zeta!
+ Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 +
+ dZeta*dZeta/fZetaWindowL2/fZetaWindowL2);
+ if (d<1) {ret=kTRUE; break;}
+ }
+ if(layer==0) { // try to see if it is reconstructable at the inner layer
+
+ // find the difference in angles
+ Float_t dPhi = TMath::Abs(fClustersLay2[iC][1] - trefLayExtr[1]);
+ // take into account boundary condition
+ if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
+
+ // find the difference in z (between linear projection from layer 2
+ // and the actual point: Dzeta= z2/r2*r1 -z1)
+ Float_t r1 = trefLayExtr[2]/TMath::Cos(trefLayExtr[0]);
+ Float_t dZeta = TMath::Cos(fClustersLay2[iC][0])*r1 - trefLayExtr[2];
+
+ // make "elliptical" cut in Phi and Zeta!
+ Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 +
+ dZeta*dZeta/fZetaWindowL1/fZetaWindowL1);
+ if (d<1) {ret=kTRUE; break;};
+ }
+}
+delete tcaRef;
+return ret;
+}
+//_________________________________________________________________________
+void AliITSTrackleterSPDEff::ReflectClusterAroundZAxisForLayer(Int_t ilayer){
+//
+// this method apply a rotation by 180 degree around the Z (beam) axis to all
+// the RecPoints in a given layer to be used to build tracklets.
+// **************** VERY IMPORTANT:: ***************
+// It must be called just after LoadClusterArrays, since afterwards the datamember
+// fClustersLay1[iC1][0] and fClustersLay1[iC1][1] are redefined using polar coordinate
+// instead of Cartesian
+//
+if(ilayer<0 || ilayer>1) {AliInfo("Input argument (ilayer) should be either 0 or 1: nothing done"); return ;}
+AliDebug(3,Form("Applying a rotation by 180 degree around z axiz to all clusters on layer %d",ilayer));
+if(ilayer==0) {
+ for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
+ fClustersLay1[iC1][0]*=-1;
+ fClustersLay1[iC1][1]*=-1;
+ }
+}
+if(ilayer==1) {
+ for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
+ fClustersLay2[iC2][0]*=-1;
+ fClustersLay2[iC2][1]*=-1;
+ }
+}
+return;
+}
+//____________________________________________________________________________
+Int_t AliITSTrackleterSPDEff::Clusters2Tracks(AliESDEvent *esd){
+// This method is used to find the tracklets.
+// It is called from AliReconstruction
+// The vertex is supposed to be associated to the Tracker (i.e. to this) already
+// The cluster is supposed to be associated to the Tracker already
+// In case Monte Carlo is required, the appropriate linking to Stack and TrackRef is attempted
+//
+ Int_t rc=1;
+ // apply cuts on the vertex quality
+ const AliESDVertex *vertex = esd->GetVertex();
+ if(vertex->GetNContributors()<fMinContVtx) return 0;
+ //
+ AliRunLoader* runLoader = AliRunLoader::Instance();
+ if (!runLoader) {
+ Error("Clusters2Tracks", "no run loader found");
+ return rc;
+ }
+ AliStack *pStack=0x0; TTree *tRefTree=0x0;
+ if(GetMC()) {
+ runLoader->LoadKinematics("read");
+ runLoader->LoadTrackRefs("read");
+ pStack= runLoader->Stack();
+ tRefTree= runLoader->TreeTR();
+ }
+ Reconstruct(pStack,tRefTree);
+
+ if (GetLightBkgStudyInParallel()) {
+ AliStack *dummy1=0x0; TTree *dummy2=0x0;
+ ReflectClusterAroundZAxisForLayer(1);
+ Reconstruct(dummy1,dummy2,kTRUE);
+ }
+ return 0;
+}
+//____________________________________________________________________________
+Int_t AliITSTrackleterSPDEff::PostProcess(AliESDEvent *){
+//
+// It is called from AliReconstruction
+//
+//
+//
+//
+ Int_t rc=0;
+ if(GetMC()) SavePredictionMC("TrackletsMCpred.root");
+ if(GetHistOn()) rc=(Int_t)WriteHistosToFile();
+ if(GetLightBkgStudyInParallel()) {
+ TString name="AliITSPlaneEffSPDtrackletBkg.root";
+ TFile* pefile = TFile::Open(name, "RECREATE");
+ rc*=fPlaneEffBkg->Write();
+ pefile->Close();
+ }
+ return rc;
+}
+//____________________________________________________________________
+void
+AliITSTrackleterSPDEff::LoadClusterArrays(TTree* itsClusterTree) {
+ // This method
+ // - gets the clusters from the cluster tree
+ // - convert them into global coordinates
+ // - store them in the internal arrays
+ // - count the number of cluster-fired chips
+
+ //AliDebug(1,"Loading clusters and cluster-fired chips ...");
+
+ fNClustersLay1 = 0;
+ fNClustersLay2 = 0;
+
+ TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
+ TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
+
+ itsClusterBranch->SetAddress(&itsClusters);
+
+ Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
+ Float_t cluGlo[3]={0.,0.,0.};
+
+ // loop over the its subdetectors
+ for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
+
+ if (!itsClusterTree->GetEvent(iIts))
+ continue;
+
+ Int_t nClusters = itsClusters->GetEntriesFast();
+
+ // number of clusters in each chip of the current module
+ Int_t layer = 0;
+
+ // loop over clusters
+ while(nClusters--) {
+ AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
+
+ layer = cluster->GetLayer();
+ if (layer>1) continue;
+
+ cluster->GetGlobalXYZ(cluGlo);
+ Float_t x = cluGlo[0];
+ Float_t y = cluGlo[1];
+ Float_t z = cluGlo[2];
+
+ if (layer==0) {
+ fClustersLay1[fNClustersLay1][0] = x;
+ fClustersLay1[fNClustersLay1][1] = y;
+ fClustersLay1[fNClustersLay1][2] = z;
+
+ for (Int_t i=0; i<3; i++)
+ fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
+ fNClustersLay1++;
+ if(fHistOn) {
+ Int_t det=cluster->GetDetectorIndex();
+ if(det<0 || det>79) {AliError("Cluster with det. index out of boundaries"); return;}
+ fhClustersInModuleLay1[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
+ }
+ }
+ if (layer==1) {
+ fClustersLay2[fNClustersLay2][0] = x;
+ fClustersLay2[fNClustersLay2][1] = y;
+ fClustersLay2[fNClustersLay2][2] = z;
+
+ for (Int_t i=0; i<3; i++)
+ fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
+ fNClustersLay2++;
+ if(fHistOn) {
+ Int_t det=cluster->GetDetectorIndex();
+ if(det<0 || det>159) {AliError("Cluster with det. index out of boundaries"); return;}
+ fhClustersInModuleLay2[det]->Fill((Double_t)cluster->GetDetLocalZ(),(Double_t)cluster->GetDetLocalX());
+ }
+ }
+
+ }// end of cluster loop
+
+ } // end of its "subdetector" loop
+ if (itsClusters) {
+ itsClusters->Delete();
+ delete itsClusters;
+ itsClusters = 0;
+ }
+ AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
+}
+//_________________________________________________________________________
+void
+AliITSTrackleterSPDEff::SetLightBkgStudyInParallel(Bool_t b) {
+// This method:
+// - set Bool_t fLightBackgroundStudyInParallel = b
+// a) if you set this kTRUE, then the estimation of the
+// SPD efficiency is done as usual for data, but in
+// parallel a light (i.e. without control histograms, etc.)
+// evaluation of combinatorial background is performed
+// with the usual ReflectClusterAroundZAxisForLayer method.
+// b) if you set this kFALSE, then you would not have a second
+// container for PlaneEfficiency statistics to be used for background
+// (fPlaneEffBkg=0). If you want to have a full evaluation of the
+// background (with all control histograms and additional data
+// members referring to the background) then you have to call the
+// method SetReflectClusterAroundZAxisForLayer(kTRUE) esplicitily
+ fLightBkgStudyInParallel=b;
+ if(fLightBkgStudyInParallel) {
+ if(!fPlaneEffBkg) fPlaneEffBkg = new AliITSPlaneEffSPD();
+ }
+ else {
+ delete fPlaneEffBkg;
+ fPlaneEffBkg=0;
+ }
+}
+//______________________________________________________________
+void AliITSTrackleterSPDEff::SetReflectClusterAroundZAxisForLayer(Int_t ilayer,Bool_t b){
+//
+// method to study residual background:
+// Input b= KTRUE --> reflect the clusters
+// ilayer (either 0 or 1) --> which SPD layers should be reflected
+//
+ if(b) {AliInfo(Form("All clusters on layer %d will be rotated by 180 deg around z",ilayer));
+ SetLightBkgStudyInParallel(kFALSE);}
+ if(ilayer==0) fReflectClusterAroundZAxisForLayer0=b; // a rotation by 180degree around the Z axis
+ else if(ilayer==1) fReflectClusterAroundZAxisForLayer1=b; // (x->-x; y->-y) to all RecPoints on a
+ else AliInfo("Nothing done: input argument (ilayer) either 0 or 1"); // given layer is applied. In such a way
+ }
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff