]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALTracker.cxx
added macros to create and read TOF TriggerMask
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.cxx
index 07f7d8dc316da677c672d763e82bd8ed6da63349..c381ff377d464098ad9e12756e0f583a1f331210 100644 (file)
@@ -53,6 +53,7 @@
 #include "AliEMCALRecParam.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
+#include "AliEMCALReconstructor.h"
 
 #include "AliEMCALTracker.h"
 
@@ -156,20 +157,22 @@ void AliEMCALTracker::InitParameters()
        
   // Check if the instance of AliEMCALRecParam exists, 
   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
+
   if(!recParam){
     AliFatal("Reconstruction parameters for EMCAL not set!");
   }
-  
-  fCutX =  recParam->GetTrkCutX();
-  fCutY =  recParam->GetTrkCutY();
-  fCutZ =  recParam->GetTrkCutZ();
-  fMaxDist =  recParam->GetTrkCutR();
-  fCutAngle =  recParam->GetTrkCutAngle();
-  fCutAlphaMin =  recParam->GetTrkCutAlphaMin();
-  fCutAlphaMax =  recParam->GetTrkCutAlphaMax();
-  fCutNITS = recParam->GetTrkCutNITS();
-  fCutNTPC = recParam->GetTrkCutNTPC();
-
+  else{
+    fCutX =  recParam->GetTrkCutX();
+    fCutY =  recParam->GetTrkCutY();
+    fCutZ =  recParam->GetTrkCutZ();
+    fMaxDist =  recParam->GetTrkCutR();
+    fCutAngle =  recParam->GetTrkCutAngle();
+    fCutAlphaMin =  recParam->GetTrkCutAlphaMin();
+    fCutAlphaMax =  recParam->GetTrkCutAlphaMax();
+    fCutNITS = recParam->GetTrkCutNITS();
+    fCutNTPC = recParam->GetTrkCutNTPC();
+  }
+       
 }
 //
 //------------------------------------------------------------------------------
@@ -192,7 +195,7 @@ TTree* AliEMCALTracker::SearchTrueMatches()
        outTree->Branch("indexT", &indexT, "indexT/I");
        outTree->Branch("label",  &label , "label/I");
        
-       Double_t dist;
+       Double_t dist=0.;
        Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
        Int_t it, nTracks = fTracks->GetEntries();
        
@@ -262,6 +265,9 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
        
        Clear("CLUSTERS");
 
+       cTree->SetBranchStatus("*",0); //disable all branches
+       cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
+
        TBranch *branch = cTree->GetBranch("EMCALECARP");
        if (!branch) {
                AliError("Can't get the branch with the EMCAL clusters");
@@ -271,13 +277,15 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
        TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
        branch->SetAddress(&clusters);
        
-       cTree->GetEvent(0);
+       //cTree->GetEvent(0);
+       branch->GetEntry(0);
        Int_t nClusters = (Int_t)clusters->GetEntries();
-       fClusters = new TObjArray(0);
+       if(fClusters) fClusters->Delete();
+       else fClusters = new TObjArray(0);
        for (Int_t i = 0; i < nClusters; i++) {
                AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
                if (!cluster) continue;
-               if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
+               if (cluster->GetClusterType() != AliVCluster::kEMCALClusterv1) continue;
                AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
                fClusters->AddLast(matchCluster);
        }
@@ -286,9 +294,9 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
         clusters->Delete();
         delete clusters;
         if (fClusters->IsEmpty())
-           AliWarning("No clusters collected");
+           AliDebug(1,"No clusters collected");
 
-       AliInfo(Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
+       AliDebug(1,Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
 
        return 0;
 }
@@ -305,8 +313,8 @@ Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
        // make sure that tracks/clusters collections are empty
        Clear("CLUSTERS");
        
-       Int_t start = esd->GetFirstEMCALCluster();
-       Int_t nClustersEMC = esd->GetNumberOfEMCALClusters();
+       Int_t start = 0;
+       Int_t nClustersEMC = esd->GetNumberOfCaloClusters();
        Int_t end = start + nClustersEMC;
        
        fClusters = new TObjArray(0);
@@ -315,14 +323,14 @@ Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
        for (i = start; i < end; i++) {
                AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
                if (!cluster) continue;
-               if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
+        if (!cluster->IsEMCAL()) continue ; 
                AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
                fClusters->AddLast(matchCluster);
        }
         if (fClusters->IsEmpty())
-           AliWarning("No clusters collected");
+           AliDebug(1,"No clusters collected");
        
-       AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
+       AliDebug(1,Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 
        return 0;
 }
@@ -340,9 +348,9 @@ Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
        Int_t nTracks = esd->GetNumberOfTracks();
        fTracks = new TObjArray(0);
        
-       Int_t i, j;
-       Bool_t isKink;
-       Double_t alpha; 
+       Int_t i=0, j=0;
+       Bool_t isKink=kFALSE;
+       Double_t alpha=0.
        for (i = 0; i < nTracks; i++) {
                AliESDtrack *esdTrack = esd->GetTrack(i);
                // set by default the value corresponding to "no match"
@@ -371,10 +379,10 @@ Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
                fTracks->AddLast(track);
        }
        if (fTracks->IsEmpty()) {
-               AliWarning("No tracks collected");
+               AliDebug(1,"No tracks collected");
        }
        
-       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
+       AliDebug(1,Form("Collected %d tracks", fTracks->GetEntries()));
 
        return 0;
 }
@@ -419,17 +427,17 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
        // IF no clusters lie within the maximum allowed distance, no matches are assigned.
        Int_t nMatches = CreateMatches();
        if (!nMatches) {
-               AliInfo(Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
+               AliDebug(1,Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
                return 0;
        }
        else {
-               AliInfo(Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
+               AliDebug(1,Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
        }
        
        // step 4:
        // when more than 1 track share the same matched cluster, only the closest one is kept.
        Int_t nRemoved = SolveCompetitions();
-       AliInfo(Form("Removed %d duplicate matches", nRemoved));
+       AliDebug(1,Form("Removed %d duplicate matches", nRemoved));
        if (nRemoved >= nMatches) {
                AliError("Removed ALL matches! Check the algorithm or data. Nothing to save");
                return 5;
@@ -437,7 +445,7 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
        
        // step 5:
        // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
-       Int_t nSaved = 0, trackID, nGood = 0, nFake = 0;
+       Int_t nSaved = 0, trackID;
        TListIter iter(fMatches);
        AliEMCALMatch *match = 0;
        while ( (match = (AliEMCALMatch*)iter.Next()) ) {
@@ -451,15 +459,8 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
                // cut on its and tpc track hits
                if(esdTrack->GetNcls(0)<=fCutNITS)continue;
                if(esdTrack->GetNcls(1)<=fCutNTPC)continue;
-
-               if (TMath::Abs(esdTrack->GetLabel()) == cluster->Label()) {
-                       esdTrack->SetEMCALcluster(cluster->Index());
-                       nGood++;
-               }
-               else {
-                       esdTrack->SetEMCALcluster(-cluster->Index());
-                       nFake++;
-               }
+             
+               esdTrack->SetEMCALcluster(cluster->Index());
                nSaved++;
        }
        /*
@@ -476,18 +477,13 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
                else {
                        AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
                        if (!cluster) continue;
-                       if (esdTrack->GetLabel() == cluster->Label()) {
-                               nGood++;
-                               esdTrack->SetEMCALcluster(cluster->Index());
-                       }
-                       else {
-                               esdTrack->SetEMCALcluster(-cluster->Index());
-                       }
+               
+                       esdTrack->SetEMCALcluster(cluster->Index());
                        nSaved++;
                }
        }
        */
-       AliInfo(Form("Saved %d matches [%d good + %d fake]", nSaved, nGood, nFake));
+       AliDebug(1,Form("Saved %d matches", nSaved));
 
        return 0;
 }
@@ -559,7 +555,7 @@ Double_t AliEMCALTracker::CheckPair
        //
        
        // TEMP
-       Bool_t isTrue = kFALSE;
+       //Bool_t isTrue = kFALSE;
 //     if (tr->GetSeedLabel() == cl->Label()) {
 //             isTrue = kTRUE;
 //     }
@@ -606,11 +602,11 @@ Double_t AliEMCALTracker::CheckPair
                x0 = 0.0;
        }
        if (fNPropSteps) {
-               Int_t i;
-               Double_t r;
+               Int_t i=0;
+               Double_t r=0.;
                cout.setf(ios::fixed);
                cout.precision(5);
-               if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
+               //if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
                for (i = 0; i < fNPropSteps; i++) {
                        r = rt + (rc - rt) * ((Double_t)(i+1)/(Double_t)fNPropSteps);
                        if (!tr->PropagateTo(r, x0, rho)){
@@ -618,9 +614,9 @@ Double_t AliEMCALTracker::CheckPair
                          return distance;
                        }
                        tr->GetXYZ(pos);
-                       if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
+               //      if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
                }
-               if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
+               //if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
        }
        else {
                // when no steps are used, no correction makes sense
@@ -656,10 +652,9 @@ Double_t AliEMCALTracker::CheckPair
        TVector3 vt(x, y, z);
        Double_t angle = TMath::Abs(vc.Angle(vt)) * TMath::RadToDeg();
        // check: where is the track?
-       Double_t r, phiT, phiC;
-       r = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
-       phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
-       phiC = vc.Phi() * TMath::RadToDeg();
+//     Double_t r    = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+//     Double_t phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
+//     Double_t phiC = vc.Phi() * TMath::RadToDeg();
        //cout << "Propagated R, phiT, phiC = " << r << ' ' << phiT << ' ' << phiC << endl;
        
        if (angle > fCutAngle) {
@@ -722,7 +717,7 @@ Double_t AliEMCALTracker::CheckPairV2
        
        Double_t distance = 2.0 * fMaxDist;
        
-       Double_t x0, rho;
+       Double_t x0=0., rho=0.;
        if (fTrackCorrMode == kTrackCorrMMB) {
                Double_t pos1[3], pos2[3], param[6];
                tr->GetXYZ(pos1);
@@ -754,7 +749,7 @@ Double_t AliEMCALTracker::CheckPairV2
        TVector3 vc(cl->X(), cl->Y(), cl->Z());
        // rotate the vector in order to put all clusters on a plane intersecting 
        // vertically the X axis; the angle depends on the sector
-       Double_t clusterRot, clusterPhi = vc.Phi() * TMath::RadToDeg();
+       Double_t clusterRot=0., clusterPhi = vc.Phi() * TMath::RadToDeg();
        if (clusterPhi < 0.0) clusterPhi += 360.0;
        if (clusterPhi < 100.0) {
                clusterRot = -90.0;
@@ -780,7 +775,7 @@ Double_t AliEMCALTracker::CheckPairV2
        // compute the 'phi' coordinate of the intersection point to 
        // the EMCAL surface
        Double_t x = vc.X();
-       Double_t y;
+       Double_t y = 0.;
        track->GetYAt(vc.X(), track->GetBz(), y);
        Double_t tmp = x*TMath::Cos(track->GetAlpha()) - y*TMath::Sin(track->GetAlpha());
        y = x*TMath::Sin(track->GetAlpha()) + y*TMath::Cos(track->GetAlpha());
@@ -826,12 +821,12 @@ Double_t AliEMCALTracker::CheckPairV3
        
        AliEMCALTrack tr(*track);
        
-       Int_t    sector;
+       Int_t    sector=-1;
        Double_t distance = 2.0 * fMaxDist;
-       Double_t dx, dy, dz;
-       Double_t phi, alpha, slope, tgtXnum, tgtXden, sectorWidth = 20.0 * TMath::DegToRad();
-       Double_t xcurr, xprop, param[6] = {0., 0., 0., 0., 0., 0.}, x0, rho, bz;
-       Double_t x[3], x1[3], x2[3];
+       Double_t dx=0., dy=0., dz=0.;
+       Double_t phi=0., alpha=0., slope=0., tgtXnum=0., tgtXden=0., sectorWidth = 20.0 * TMath::DegToRad();
+       Double_t xcurr=0., xprop=0., param[6] = {0., 0., 0., 0., 0., 0.}, x0=0., rho=0., bz=0.;
+       Double_t x[3]= {0., 0., 0.}, x1[3]= {0., 0., 0.}, x2[3]= {0., 0., 0.};
        
        // get initial track position
        xcurr = tr.GetX();
@@ -853,7 +848,7 @@ Double_t AliEMCALTracker::CheckPairV3
        if (!tr.GetXYZAt(xprop, bz, x2)) return distance;
        //AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
        rho = param[0]*param[4];
-       x0 = param[1];
+       x0  = param[1];
        if (!tr.PropagateTo(xprop, x0, rho)) return distance;
        //if (!tr.PropagateTo(xprop, 0.0, 0.0)) return distance;
        
@@ -877,8 +872,8 @@ Bool_t AliEMCALTracker::PropagateToEMCAL(AliEMCALTrack *tr)
        // Propagates the track to the proximity of the EMCAL surface
        //
        
-       Double_t xcurr, xtemp, xprop = 438.0, step = 10.0, param[6], x0, rho, bz;
-       Double_t x1[3], x2[3];
+       Double_t xcurr=0., xtemp=0., xprop = 438.0, step = 10.0, param[6]= {0., 0., 0., 0., 0., 0.}, x0=0., rho=0., bz=0.;
+       Double_t x1[3]= {0., 0., 0.}, x2[3]= {0., 0., 0.};
        
        // get initial track position
        xcurr = tr->GetX();
@@ -922,11 +917,11 @@ Int_t AliEMCALTracker::CreateMatches()
        
        // initialize counters and indexes
        Int_t count = 0;
-       Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
-       Int_t it, nTracks = fTracks->GetEntries();
+       Int_t ic=0, nClusters = (Int_t)fClusters->GetEntries();
+       Int_t it=0, nTracks = fTracks->GetEntries();
        
        // external loop on clusters, internal loop on tracks
-       Double_t dist;
+       Double_t dist=0.;
        for (ic = 0; ic < nClusters; ic++) {
                AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
                for (it = 0; it < nTracks; it++) {
@@ -944,35 +939,6 @@ Int_t AliEMCALTracker::CreateMatches()
                        }
                }
        }
-               
-       /*
-       // loop on clusters and tracks
-       Int_t icBest;
-       Double_t dist, distBest;
-       for (it = 0; it < nTracks; it++) {
-               AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
-               if (!track) continue;
-               icBest = -1;
-               distBest = fMaxDist;
-               for (ic = 0; ic < nClusters; ic++) {
-                       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
-                       if (!cluster) continue;
-                       dist = CheckPair(track, cluster);
-                       if (dist < distBest) {
-                               distBest = dist;
-                               icBest = ic;
-                       }
-               }
-               if (icBest >= 0) {
-                       track->SetMatchedClusterIndex(icBest);
-                       track->SetMatchedClusterDist(distBest);
-                       count++;
-               }
-               else {
-                       track->SetMatchedClusterIndex(-1);
-               }
-       }
-       */
        
        return count;
 }
@@ -998,8 +964,8 @@ Int_t AliEMCALTracker::SolveCompetitions()
        Int_t count = 0;
        
        // initialize flags to check repetitions
-       Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
-       Int_t it, nTracks = fTracks->GetEntries();
+       Int_t ic=0, nClusters = (Int_t)fClusters->GetEntries();
+       Int_t it=0, nTracks = fTracks->GetEntries();
        Bool_t *usedC = new Bool_t[nClusters];
        Bool_t *usedT = new Bool_t[nTracks];
        for (ic = 0; ic < nClusters; ic++) usedC[ic] = kFALSE;
@@ -1021,29 +987,6 @@ Int_t AliEMCALTracker::SolveCompetitions()
                }
        }
        
-       /*
-       Int_t it1, it2, nTracks = (Int_t)fTracks->GetEntries();
-       AliEMCALTrack *track1 = 0, *track2 = 0;
-       for (it1 = 0; it1 < nTracks; it1++) {
-               track1 = (AliEMCALTrack*)fTracks->At(it1);
-               if (!track1) continue;
-               if (track1->GetMatchedClusterIndex() < 0) continue;
-               for (it2 = it1+1; it2 < nTracks; it2++) {
-                       track2 = (AliEMCALTrack*)fTracks->At(it2);
-                       if (!track2) continue;
-                       if (track2->GetMatchedClusterIndex() < 0) continue;
-                       if (track1->GetMatchedClusterIndex() != track2->GetMatchedClusterIndex()) continue;
-                       count++;
-                       if (track1->GetMatchedClusterDist() < track2->GetMatchedClusterDist()) {
-                               track2->SetMatchedClusterIndex(-1);
-                       }
-                       else if (track2->GetMatchedClusterDist() < track1->GetMatchedClusterDist()) {
-                               track1->SetMatchedClusterIndex(-1);
-                       }
-               }
-       }
-       */
-       
        delete [] usedC;
        delete [] usedT;
 
@@ -1062,12 +1005,29 @@ void AliEMCALTracker::UnloadClusters()
        
        Clear();
 }
+//
+//------------------------------------------------------------------------------
+//
+TVector3 AliEMCALTracker::FindExtrapolationPoint(Double_t x,Double_t y,Double_t z, AliESDtrack *track)
+{
+  //Method to determine extrapolation point of track at location x,y,z
+  AliEMCALTrack *tr = new AliEMCALTrack(*track);
+  TVector3 error(-100.,-100.,-100.);
+  if (!tr->PropagateToGlobal(x,y,z, 0.0, 0.0)) {
+    return error;
+  }
+  Double_t pos[3]= {0., 0., 0.};
+  tr->GetXYZ(pos);
+  TVector3 ExTrPos(pos[0],pos[1],pos[2]);
+  return ExTrPos;
+}
+
 //
 //------------------------------------------------------------------------------
 //
 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
   : fIndex(index),
-    fLabel(recPoint->GetPrimaryIndex()),
+    fLabel(recPoint->GetPrimaryIndex()), //wrong!  fixed below
     fX(0.),
     fY(0.),
     fZ(0.)
@@ -1082,6 +1042,17 @@ AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCA
        fX = clpos.X();
        fY = clpos.Y();
        fZ = clpos.Z();
+
+       //AliEMCALRecPoint stores the track labels in the parents
+       //list, sorted according to the fractional contribution to the
+       //RecPoint.  The zeroth parent gave the highest contribution
+       Int_t multparent = 0;
+        Int_t *parents = recPoint->GetParents(multparent);
+        if(multparent > 0)
+          fLabel = parents[0];
+        else
+          fLabel = -1;
+
 }
 //
 //------------------------------------------------------------------------------
@@ -1097,7 +1068,7 @@ AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDC
        // Translates an AliESDCaloCluster object into the internal format.
        // Index of passed cluster in its native array must be specified.
        //
-       Float_t clpos[3];
+       Float_t clpos[3]= {0., 0., 0.};
        caloCluster->GetPosition(clpos);
        
        fX = (Double_t)clpos[0];