ITS/AliITSFindClustersV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.2.2.1  2002/10/14 13:14:08  hristov
  19 Updating VirtualMC to v3-09-02
  20
  21 Revision 1.2  2002/09/23 13:25:52  hristov
  22 Typo corrected (HP)
  23
  24 Revision 1.1  2002/09/09 17:36:05  nilsen
  25 new TTask to replace non-working AliITSFindClusterV2.C macro.
  26
  27 */
  28 #include <TROOT.h>
  29 #include <TFile.h>
  30 #include <TTree.h>
  31 #include <TBranch.h>
  32 #include <TClonesArray.h>
  33 #include <TParticle.h>
  34
  35 #include "AliRun.h"
  36 #include "AliHeader.h"
  37
  38 #include "AliITS.h"
  39 #include "AliITSRecPoint.h"
  40 #include "AliITSFindClustersV2.h"
  41 #include "AliITSclusterV2.h"
  42 #include "AliITSgeom.h"
  43
  44 ClassImp(AliITSFindClustersV2)
  45
  46 //______________________________________________________________________
  47 AliITSFindClustersV2::AliITSFindClustersV2(){
  48     // Default constructor.
  49     // Inputs:
  50     //  none.
  51     // Outputs:
  52     //   none.
  53     // Return:
  54     //    A zero-ed constructed AliITSFindClustersV2 class.
  55
  56     fAr          = 0;
  57     fDeletfAr     = kFALSE; // fAr=0 dont't delete it.
  58     fGeom        = 0;
  59     fInFileName  = 0;
  60     fOutFileName = 0;
  61     fIn          = 0;
  62     fOut         = 0;
  63     fSlowFast    = kFALSE;  // slow simulation
  64     fInit        = kFALSE;  // Init failed
  65 }
  66 //______________________________________________________________________
  67 AliITSFindClustersV2::AliITSFindClustersV2(const TString infile,
  68                                            const TString outfile){
  69     // Standard constructor.
  70     // Inputs:
  71     //  const TString infile   Input file name where the RecPoints are
  72     //                         to be read from.
  73     //  const TString outfile  Output file where V2 tracking clulsters
  74     //                         are to be written. if =="" writen to the
  75     //                         same file as input.
  76     // Outputs:
  77     //   none.
  78     // Return:
  79     //    A standardly constructed AliITSFindClustersV2 class.
  80
  81     fAr          = 0;
  82     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
  83     fGeom        = 0;
  84     fInFileName  = 0;
  85     fOutFileName = 0;
  86     fIn          = 0;
  87     fOut         = 0;
  88     fSlowFast    = kFALSE;  // slow simulation
  89     fInit        = kFALSE;  // Init failed
  90
  91     fInFileName = new TString(infile);
  92     if(outfile.CompareTo("")!=0){
  93         fOutFileName = new TString(outfile);
  94     } // end if outfile.CompareeTo("")!=0
  95
  96     if(fOutFileName!=0){
  97         fIn = new TFile(fInFileName->Data(),"READ");
  98         fOut = new TFile(fOutFileName->Data(),"UPDATE");
  99     }else{ // open fIn file for update
 100         fIn = new TFile(fInFileName->Data(),"UPDATE");
 101     } // end if
 102
 103     fAr  = (AliRun*) fIn->Get("gAlice");
 104     if(!fAr){
 105         Warning("AliITSFindClusterV2",
 106                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
 107         return;
 108     } // end if !fAr
 109     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
 110
 111     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 112     if(!its){
 113         Warning("AliITSFindClusterV2",
 114                 "Can't fine the ITS in gAlice. Aborting.");
 115         return;
 116     } // end if !its
 117     fGeom = its->GetITSgeom();
 118     if(!fGeom){
 119         Warning("AliITSFindClusterV2",
 120                 "Can't fine the ITS geometry in gAlice. Aborting.");
 121         return;
 122     } // end if !fGeom
 123
 124     if(fOut) fOut->cd();
 125     fInit = kTRUE;
 126 }
 127 //______________________________________________________________________
 128 AliITSFindClustersV2::AliITSFindClustersV2(TFile *in,
 129                                            TFile *out){
 130     // Standard constructor.
 131     // Inputs:
 132     //  const TFile *in   Input file where the RecPoints are
 133     //                         to be read from.
 134     //  const Tfile *out  Output file where V2 tracking clulsters
 135     //                         are to be written. if ==0 writen to the
 136     //                         same file as input.
 137     // Outputs:
 138     //   none.
 139     // Return:
 140     //    A standardly constructed AliITSFindClustersV2 class.
 141
 142     fAr          = 0;
 143     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
 144     fGeom        = 0;
 145     fInFileName  = 0;
 146     fOutFileName = 0;
 147     fIn          = 0;
 148     fOut         = 0;
 149     fSlowFast    = kFALSE;  // slow simulation
 150     fInit        = kFALSE;  // Init failed
 151
 152     fIn  = in;
 153     fOut = out;
 154     fAr  = (AliRun*) fIn->Get("gAlice");
 155     if(!fAr){
 156         Warning("AliITSFindClusterV2",
 157                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
 158         return;
 159     } // end if !fAr
 160     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
 161     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 162     if(!its){
 163         Warning("AliITSFindClusterV2",
 164                 "Can't fine the ITS in gAlice. Aborting.");
 165         return;
 166     } // end if !its
 167     fGeom = its->GetITSgeom();
 168     if(!fGeom){
 169         Warning("AliITSFindClusterV2",
 170                 "Can't fine the ITS geometry in gAlice. Aborting.");
 171         return;
 172     } // end if !fGeom
 173
 174     if(fOut) fOut->cd();
 175     fInit = kTRUE;
 176 }
 177 //______________________________________________________________________
 178 AliITSFindClustersV2::AliITSFindClustersV2(AliRun *ar,
 179                                            const TString outfile){
 180     // Standard constructor.
 181     // Inputs:
 182     //  const TString infile   Input file name where the RecPoints are
 183     //                         to be read from.
 184     //  const TString outfile  Output file where V2 tracking clulsters
 185     //                         are to be written. if =="" writen to the
 186     //                         same file as input.
 187     // Outputs:
 188     //   none.
 189     // Return:
 190     //    A standardly constructed AliITSFindClustersV2 class.
 191
 192     fAr          = 0;
 193     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
 194     fGeom        = 0;
 195     fInFileName  = 0;
 196     fOutFileName = 0;
 197     fIn          = 0;
 198     fOut         = 0;
 199     fSlowFast    = kFALSE;  // slow simulation
 200     fInit        = kFALSE;  // Init failed
 201
 202     if(outfile.CompareTo("")!=0){
 203         fOutFileName = new TString(outfile);
 204     } // end if outfile.CompareeTo("")!=0
 205
 206     if(fOutFileName!=0){
 207         fOut = new TFile(fOutFileName->Data(),"UPDATE");
 208     } // end if
 209
 210     fAr  = ar;
 211     if(!fAr){
 212         Warning("AliITSFindClusterV2",
 213                 "ar==0. Aborting.");
 214         return;
 215     } // end if !fAr
 216     AliITS *its = (AliITS*) fAr->GetModule("ITS");
 217     if(!its){
 218         Warning("AliITSFindClusterV2",
 219                 "Can't fine the ITS in gAlice. Aborting.");
 220         return;
 221     } // end if !its
 222     fGeom = its->GetITSgeom();
 223     if(!fGeom){
 224         Warning("AliITSFindClusterV2",
 225                 "Can't fine the ITS geometry in gAlice. Aborting.");
 226         return;
 227     } // end if !fGeom
 228
 229     if(fOut) fOut->cd();
 230     fInit = kTRUE;
 231 }
 232 //______________________________________________________________________
 233 AliITSFindClustersV2::~AliITSFindClustersV2(){
 234     // Default constructor.
 235     // Inputs:
 236     //  none.
 237     // Outputs:
 238     //   none.
 239     // Return:
 240     //    A destroyed AliITSFindclustersV2 class.
 241
 242     fGeom = 0; // Deleted when AliRun/ITS is deleted.
 243     if(fInFileName!=0){ // input file opened by AliITSFindClustersV2
 244         if(fIn!=0){
 245             if(fIn->IsOpen()) fIn->Close();
 246             delete fIn;
 247             fIn = 0;
 248         } // end if
 249         delete fInFileName;
 250         fInFileName = 0;
 251     } // end if
 252
 253     if(fOutFileName!=0){ // input file opened by AliITSFindClustersV2
 254         if(fOut!=0){
 255             if(fOut->IsOpen()) fOut->Close();
 256             delete fOut;
 257             fOut = 0;
 258         } // end if
 259         delete fOutFileName;
 260         fOutFileName = 0;
 261     } // end if
 262     if(fDeletfAr && !fAr){
 263         cout << "deleting fAr."<< endl;
 264         delete fAr;
 265         fAr = 0;
 266         cout << "fAr deleted OK."<< endl;
 267     } // end if fDeletfAr
 268 }
 269 //______________________________________________________________________
 270 void AliITSFindClustersV2::Exec(const Option_t *opt){
 271     // Main FindclustersV2 function.
 272     // Inputs:
 273     //      Option_t * opt   list of subdetector to digitize. =0 all.
 274     // Outputs:
 275     //      none.
 276     // Return:
 277     //      none.
 278     Char_t name[50];
 279
 280     if(!fInit){
 281         Warning("Exec","Initilization not succesfull. Aborting.");
 282         return;
 283     } // end if !fInit
 284
 285     fGeom->Write();
 286
 287     fAr->GetEvent(0);
 288     TTree *pTree = fAr->TreeR();
 289     if(!pTree){
 290         Warning("Exec","Error getting TreeR. TreeR=%p",pTree);
 291         return;
 292     } // end if !pTree
 293     TBranch *branch = 0;
 294     if(fSlowFast) sprintf(name,"ITSRecPointsF");
 295     else sprintf(name,"ITSRecPoints");
 296     branch = pTree->GetBranch(name);
 297     if(!branch){
 298         Warning("Exec","Can't find branch %s in TreeR fSlowFast=%d",
 299                 name,fSlowFast);
 300         return;
 301     } // end if !branch
 302     TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
 303     branch->SetAddress(&points);
 304     Int_t nentr = (Int_t) branch->GetEntries();
 305
 306     if(fOut) fOut->cd();
 307     TClonesArray *cluster = new TClonesArray("AliITSclusterV2",10000);
 308     sprintf(name,"TreeC_ITS_%d",fAr->GetHeader()->GetEvent());
 309     TTree *cTree = new TTree(name,"ITS clusters");
 310     cTree->Branch("Clusters",&cluster);
 311     TClonesArray &cl = *cluster;
 312
 313     Float_t lp[5];
 314     Int_t lab[6];
 315     Int_t i,j,lay,lad,det,nclusters=0,ncl;
 316     Float_t kmip,x,y,zshift,yshift;
 317     Double_t rot[9];
 318     AliITSRecPoint *p;
 319     TParticle *part;
 320
 321     for(i=0;i<nentr;i++){
 322         points->Clear();
 323         branch->GetEvent(i);
 324         ncl = points->GetEntriesFast();
 325         if(ncl<=0) {
 326             cTree->Fill();
 327             continue;
 328         } // end if ncl<=0
 329         nclusters += ncl;
 330         fGeom->GetModuleId(i,lay,lad,det);
 331         fGeom->GetTrans(i,x,y,zshift);
 332         fGeom->GetRotMatrix(i,rot);
 333         yshift = x*rot[0] + y*rot[1];
 334         kmip = 1.0; // fGeom->GetModuletype(i)==kSPD
 335         if(fGeom->GetModuleType(i)==kSDD) kmip = 280.0;
 336         if(fGeom->GetModuleType(i)==kSSD) kmip = 38.0;
 337         for(j=0;j<ncl;j++){
 338             p = (AliITSRecPoint*)(points->UncheckedAt(j));
 339             lp[0] = - (p->GetX() + yshift);
 340             if(lay==1) lp[0] = -lp[0];
 341             lp[1] = p->GetZ() + zshift;
 342             lp[2] = p->GetSigmaX2();
 343             lp[3] = p->GetSigmaZ2();
 344             lp[4] = p->GetQ()/kmip;
 345             lab[0] = p->GetLabel(0);
 346             lab[1] = p->GetLabel(1);
 347             lab[2] = p->GetLabel(2);
 348             lab[3] = i;
 349             lad = lab[0];
 350             if(lad>=0){
 351                 part = (TParticle*) fAr->Particle(lad);
 352                 lad = -3;
 353                 while(part->P() < 0.005){
 354                     if(part->GetFirstMother()<0){
 355                         Warning("Exec","Primary momentum: %e",part->P());
 356                         break;
 357                     } // end if part->GetFirstMother()<0
 358                     lad = part->GetFirstMother();
 359                     part = (TParticle*) fAr->Particle(lad);
 360                 } // end while part->P() < 0.005
 361                 if(lab[1]<0) lab[1] = lad;
 362                 else if(lab[2]<0) lab[2] = lad;
 363                 else Warning("Exec","No empty lables!");
 364             } // end if lab>=0
 365             new(cl[j]) AliITSclusterV2(lab,lp);
 366         } // end for j
 367         cTree->Fill();
 368         cluster->Delete();
 369         points->Delete();
 370     } // end for i
 371     cTree->Write();
 372
 373     delete cTree;
 374     delete cluster;
 375     delete points;
 376 }