MFT/MergeClustersMFT.C

   1 #include "AliMFTCluster.h"
   2 #include "TFile.h"
   3 #include "AliMFTCluster.h"
   4 #include "TObjArray.h"
   5 #include "TClonesArray.h"
   6 #include "TTree.h"
   7 #include "AliMFTConstants.h"
   8 #include "TBranch.h"
   9 #include "TRandom.h"
  10 #include "TMath.h"
  11
  12 TFile *fileInPileUp=0, *fileInSignal, *fileInUnderlying, *fileOutPileUp=0, *fileOutTotal;
  13
  14 enum {kMergePileUpClusters_pp, kMergePileUpClusters_PbPb, kMergeAllClusters_pp, kMergeAllClusters_PbPb};
  15 enum {k_pp, k_PbPb};
  16
  17 const Char_t *fPileUpDir[2] = {"$UnderlyingDirPP/spreadVtx",      // local directory with pile-up event p-p
  18                                "$UnderlyingDirPbPb/spreadVtx"};   // local directory with pile-up event Pb-Pb
  19
  20 const Int_t fNRunsPileUp_pp   = 100;
  21 const Int_t fNRunsPileUp_PbPb =   3;
  22
  23 const Int_t fNRunsPileUpAvailable[2] = {500, 500};    // pp, PbPb
  24
  25 TClonesArray *fRecPointsPerPlaneIn[AliMFTConstants::fNMaxPlanes] = {0};
  26 TClonesArray *fRecPointsPerPlaneOut[AliMFTConstants::fNMaxPlanes] = {0};
  27
  28 void MergePileUpClusters(Int_t *runs, Int_t nRunsPileUp, const Char_t *pileUpDir, Int_t nEvents);
  29 void MergeAllClusters();
  30
  31 //====================================================================================================================================================
  32
  33 void MergeClustersMFT(Int_t seed = 12345,
  34                       Int_t option = 9999,
  35                       Int_t nEvents = 9999,
  36                       Int_t maxPileUp = 9999) {
  37
  38   gRandom->SetSeed(seed);
  39
  40   if (option==kMergePileUpClusters_pp) {
  41
  42     Int_t runs_pp[fNRunsPileUp_pp] = {0};
  43     Int_t nFiles_pp = 0;
  44     while (nFiles_pp<TMath::Min(fNRunsPileUp_pp,maxPileUp)) {
  45       Int_t run = gRandom->Integer(fNRunsPileUpAvailable[k_pp]);
  46       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",fPileUpDir[k_pp],run));
  47       if (!fileInPileUp->IsOpen() || !fileInPileUp->cd("Event0")) {
  48         delete fileInPileUp;
  49         continue;
  50       }
  51       Bool_t runAlreadyLoaded=kFALSE;
  52       for (Int_t iAddedRun=0; iAddedRun<nFiles_pp; iAddedRun++) {
  53         if (run==runs_pp[iAddedRun]) {
  54           runAlreadyLoaded=kTRUE;
  55           break;
  56         }
  57       }
  58       if (!runAlreadyLoaded) runs_pp[nFiles_pp++] = run;
  59       delete fileInPileUp;
  60     }
  61     MergePileUpClusters(runs_pp, TMath::Min(fNRunsPileUp_pp,maxPileUp), fPileUpDir[k_pp], nEvents);
  62   }
  63
  64   else if (option==kMergePileUpClusters_PbPb) {
  65     Int_t runs_PbPb[fNRunsPileUp_PbPb] = {0};
  66     Int_t nFiles_PbPb = 0;
  67     while (nFiles_PbPb<TMath::Min(fNRunsPileUp_PbPb,maxPileUp)) {
  68       Int_t run = gRandom->Integer(fNRunsPileUpAvailable[k_PbPb]);
  69       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",fPileUpDir[k_PbPb],run));
  70       if (!fileInPileUp->IsOpen() || !fileInPileUp->cd("Event0")) {
  71         delete fileInPileUp;
  72         continue;
  73       }
  74       Bool_t runAlreadyLoaded=kFALSE;
  75       for (Int_t iAddedRun=0; iAddedRun<nFiles_PbPb; iAddedRun++) {
  76         if (run==runs_PbPb[iAddedRun]) {
  77           runAlreadyLoaded=kTRUE;
  78           break;
  79         }
  80       }
  81       if (!runAlreadyLoaded) runs_PbPb[nFiles_PbPb++] = run;
  82       delete fileInPileUp;
  83     }
  84     MergePileUpClusters(runs_PbPb, TMath::Min(fNRunsPileUp_PbPb,maxPileUp), fPileUpDir[k_PbPb], nEvents);
  85   }
  86
  87   if (option>=kMergeAllClusters_pp) MergeAllClusters();
  88
  89 }
  90
  91 //====================================================================================================================================================
  92
  93 void MergePileUpClusters(Int_t *runs, Int_t nRunsPileUp, const Char_t *pileUpDir, Int_t nEvents) {
  94
  95   fileOutPileUp = new TFile("./MFT.RecPoints.PileUp.root", "recreate");
  96
  97   Bool_t eventExist = kTRUE;
  98   Int_t iEv=0;
  99
 100   while (eventExist && iEv<nEvents) {
 101
 102     for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
 103       fRecPointsPerPlaneOut[iPlane] = new TClonesArray("AliMFTCluster");
 104     }
 105
 106     if (iEv) fileOutPileUp = new TFile("./MFT.RecPoints.PileUp.root", "update");
 107
 108     printf("Merging Event %d\n",iEv);
 109
 110     fileOutPileUp-> cd();
 111     TTree *treeOutPileUp = new TTree("TreeR_OutPileUp", "Reconstructed Points Container");
 112
 113     Int_t nFilesAdded = 0;
 114     while (nFilesAdded<nRunsPileUp) {
 115       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",pileUpDir,runs[nFilesAdded++]));
 116       printf("Merging file %s\n", fileInPileUp->GetName());
 117       if (!fileInPileUp->cd(Form("Event%d",iEv))) {
 118         eventExist = kFALSE;
 119         break;
 120       }
 121       for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
 122         fRecPointsPerPlaneIn[iPlane] = new TClonesArray("AliMFTCluster");
 123       }
 124       TTree *treeInPileUp = (TTree*) gDirectory->Get("TreeR");
 125       treeInPileUp -> SetName("TreeR_InPileUp");
 126       Int_t iPlane=0;
 127       while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
 128         //      printf("Plane %02d\n",iPlane);
 129         treeInPileUp ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
 130         TBranch *branch = treeOutPileUp->GetBranch(Form("Plane_%02d",iPlane));
 131         if (!branch) treeOutPileUp->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane]));
 132         iPlane++;
 133       }
 134       iPlane=0;
 135       treeInPileUp -> GetEntry(0);
 136       while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
 137         Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
 138         for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 139           //      printf("Cluster %4d\n",iCluster);
 140           AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
 141           new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 142         }
 143         iPlane++;
 144       }
 145       for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 146         delete fRecPointsPerPlaneIn[jPlane];
 147       }
 148       delete treeInPileUp;
 149       fileInPileUp -> Close();
 150       delete fileInPileUp;
 151     }
 152     if (eventExist) {
 153       treeOutPileUp -> Fill();
 154       fileOutPileUp -> mkdir(Form("Event%d",iEv));
 155       fileOutPileUp -> cd(Form("Event%d",iEv));
 156       treeOutPileUp -> SetName("TreeR");
 157       treeOutPileUp -> Write();
 158       delete treeOutPileUp;
 159       fileOutPileUp -> Close();
 160       delete fileOutPileUp;
 161       for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
 162         delete fRecPointsPerPlaneOut[iPlane];
 163       }
 164       iEv++;
 165     }
 166   }
 167
 168 }
 169
 170 //====================================================================================================================================================
 171
 172 void MergeAllClusters() {
 173
 174   fileOutTotal = new TFile("./MFT.RecPoints.root", "recreate");
 175
 176   Bool_t eventExist = kTRUE;
 177   Int_t iEv=0;
 178
 179   while (eventExist) {
 180
 181     for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
 182       fRecPointsPerPlaneOut[iPlane] = new TClonesArray("AliMFTCluster");
 183     }
 184
 185     if (iEv) fileOutTotal = new TFile("./MFT.RecPoints.root", "update");
 186
 187     printf("Merging Event %d\n",iEv);
 188
 189     fileOutTotal-> cd();
 190     TTree *treeOutTotal = new TTree("TreeR_Out", "Reconstructed Points Container");
 191
 192     // Adding Pile-Up Clusters
 193
 194     fileInPileUp = new TFile("./MFT.RecPoints.PileUp.root");
 195     printf("Merging file %s\n", fileInPileUp->GetName());
 196     if (!fileInPileUp->cd(Form("Event%d",iEv))) {
 197       eventExist = kFALSE;
 198       break;
 199     }
 200     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 201       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
 202     }
 203     TTree *treeInPileUp = (TTree*) gDirectory->Get("TreeR");
 204     treeInPileUp -> SetName("TreeR_InPileUp");
 205     Int_t iPlane=0;
 206     while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
 207       //        printf("Plane %02d\n",iPlane);
 208       treeInPileUp ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
 209       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
 210       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane]));
 211       iPlane++;
 212     }
 213     iPlane=0;
 214     treeInPileUp -> GetEntry(0);
 215     while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
 216       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
 217       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 218         //        printf("Cluster %4d\n",iCluster);
 219         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
 220         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 221       }
 222       iPlane++;
 223     }
 224     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 225       delete fRecPointsPerPlaneIn[jPlane];
 226     }
 227     delete treeInPileUp;
 228     fileInPileUp -> Close();
 229     delete fileInPileUp;
 230
 231     // Adding Underlying Event Clusters
 232
 233     fileInUnderlying = new TFile("./MFT.RecPoints.Underlying.root");
 234     printf("Merging file %s\n", fileInUnderlying->GetName());
 235     if (!fileInUnderlying->cd(Form("Event%d",iEv))) {
 236       eventExist = kFALSE;
 237       break;
 238     }
 239     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 240       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
 241     }
 242     TTree *treeInUnderlying = (TTree*) gDirectory->Get("TreeR");
 243     treeInUnderlying -> SetName("TreeR_InUnderlying");
 244     iPlane=0;
 245     while (treeInUnderlying->GetBranch(Form("Plane_%02d",iPlane))) {
 246       //        printf("Plane %02d\n",iPlane);
 247       treeInUnderlying ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
 248       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
 249       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane]));
 250       iPlane++;
 251     }
 252     iPlane=0;
 253     treeInUnderlying -> GetEntry(0);
 254     while (treeInUnderlying->GetBranch(Form("Plane_%02d",iPlane))) {
 255       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
 256       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 257         //        printf("Cluster %4d\n",iCluster);
 258         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
 259         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 260       }
 261       iPlane++;
 262     }
 263     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 264       delete fRecPointsPerPlaneIn[jPlane];
 265     }
 266     delete treeInUnderlying;
 267     fileInUnderlying -> Close();
 268     delete fileInUnderlying;
 269
 270     // Adding Signal Clusters
 271
 272     fileInSignal = new TFile("./MFT.RecPoints.Signal.root");
 273     printf("Merging file %s\n", fileInSignal->GetName());
 274     if (!fileInSignal->cd(Form("Event%d",iEv))) {
 275       eventExist = kFALSE;
 276       break;
 277     }
 278     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 279       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
 280     }
 281     TTree *treeInSignal = (TTree*) gDirectory->Get("TreeR");
 282     treeInSignal -> SetName("TreeR_InSignal");
 283     iPlane=0;
 284     while (treeInSignal->GetBranch(Form("Plane_%02d",iPlane))) {
 285       //        printf("Plane %02d\n",iPlane);
 286       treeInSignal ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
 287       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
 288       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane]));
 289       iPlane++;
 290     }
 291     iPlane=0;
 292     treeInSignal -> GetEntry(0);
 293     while (treeInSignal->GetBranch(Form("Plane_%02d",iPlane))) {
 294       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
 295       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
 296         //        printf("Cluster %4d\n",iCluster);
 297         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
 298         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
 299       }
 300       iPlane++;
 301     }
 302     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 303       delete fRecPointsPerPlaneIn[jPlane];
 304     }
 305     delete treeInSignal;
 306     fileInSignal -> Close();
 307     delete fileInSignal;
 308
 309     if (eventExist) {
 310       treeOutTotal -> Fill();
 311       fileOutTotal -> mkdir(Form("Event%d",iEv));
 312       fileOutTotal -> cd(Form("Event%d",iEv));
 313       treeOutTotal -> SetName("TreeR");
 314       treeOutTotal -> Write();
 315       delete treeOutTotal;
 316       fileOutTotal -> Close();
 317       delete fileOutTotal;
 318       for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
 319         delete fRecPointsPerPlaneOut[jPlane];
 320       }
 321       iEv++;
 322     }
 323
 324   }
 325
 326 }
 327
 328 //====================================================================================================================================================
 329
 330 //  LocalWords:  nFiles