MUON/MUONlist2digits.C

   1 void MUONlist2digits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and do some analysis.
   7 //
   8 /////////////////////////////////////////////////////////////////////////
   9
  10 // Dynamically link some shared libs
  11
  12
  13     if (gClassTable->GetID("AliRun") < 0) {
  14         gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  15         gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  16         gSystem->Load("$ROOTSYS/lib/libEG.so");
  17         gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  18         gSystem->Load("libGeant3Dummy.so");    //a dummy version of Geant3
  19         gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
  20         gSystem->Load("libgalice.so");         // the standard Alice classes
  21     }
  22
  23 // Connect the Root Galice file containing Geometry, Kine and Hits
  24
  25    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  26    if (file) file->Close();
  27    file = new TFile("galice.root","UPDATE");
  28    file->ls();
  29    file->Map();
  30
  31    printf ("I'm after Map \n");
  32
  33 // Get AliRun object from file or create it if not on file
  34
  35    if (!gAlice) {
  36       gAlice = (AliRun*)file->Get("gAlice");
  37       if (gAlice) printf("AliRun object found on file\n");
  38       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  39    }
  40    printf ("I'm after gAlice \n");
  41
  42        AliMUON *MUON  = gAlice->GetDetector("MUON");
  43
  44    //   AliMUONchamber&  iChamber;
  45    AliMUONchamber*  iChamber;
  46    AliMUONsegmentation*  segmentation;
  47
  48    const Float_t rmin[10] = {17.5,17.5,23.5,23.5,33.5,33.5,43,43,50,50};
  49    const Float_t rmax[10] = {91.5,91.5,122.5,122.5,173,173,221,221,256.5,256.5};
  50
  51
  52    Int_t Npx[10];
  53    Int_t Npy[10];
  54
  55
  56    for ( int i=0;i<10;i++) {
  57         Npx[i]=(Int_t)rmax[i]/0.75;
  58         Npy[i]=(Int_t)rmax[i]/0.5;
  59    }
  60
  61        Int_t nxmax=1026;
  62        Int_t nymax=1026;
  63        AliMUONlist *elem[10][1026][1026];
  64        TObjArray *obj=new TObjArray;
  65
  66        Int_t digits[3];
  67 //
  68 //   Loop over events
  69 //
  70
  71    Int_t Nh=0;
  72    Int_t Nh1=0;
  73    for (int nev=0; nev<= evNumber2; nev++) {
  74        Int_t nparticles = gAlice->GetEvent(nev);
  75        cout << "nev         " << nev <<endl;
  76        cout << "nparticles  " << nparticles <<endl;
  77        if (nev < evNumber1) continue;
  78        if (nparticles <= 0) return;
  79
  80        TTree *TH = gAlice->TreeH();
  81        Int_t ntracks = TH->GetEntries();
  82        Int_t Nc=0;
  83        //
  84        Int_t counter=0;
  85
  86     for (int icat=0;icat<nCathode;icat++) {
  87
  88        //initialize the arrray of pointers !!!!!
  89        for (int i=0; i<10; i++) {
  90          for (int j=0; j<nymax; j++) {
  91            for (int k=0; k<nxmax; k++) {
  92                elem[i][j][k]=0;
  93            }
  94          }
  95        }
  96
  97
  98 //
  99 //   Loop over events
 100 //
 101        for (Int_t track=0; track<ntracks;track++) {
 102            gAlice->ResetHits();
 103            Int_t nbytes += TH->GetEvent(track);
 104            if (MUON)  {
 105                for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
 106                    mHit;
 107                    mHit=(AliMUONhit*)MUON->NextHit())
 108                {
 109                    Int_t   nch   = mHit->fChamber;  // chamber number
 110                    Float_t x     = mHit->fX;        // x-pos of hit
 111                    Float_t y     = mHit->fY;        // y-pos
 112
 113                    iChamber = &(MUON->Chamber(nch-1));
 114                    response=iChamber->GetResponseModel();
 115
 116                    Int_t nsec=iChamber->Nsec();
 117                    //                   Int_t rmin = (Int_t)iChamber->frMin;
 118                    //                   Int_t rmax = (Int_t)iChamber->frMax;
 119 //
 120 //
 121                    for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit);
 122                         mPad;
 123                         mPad=(AliMUONcluster*)MUON->NextPad())
 124                    {
 125                         Int_t cathode     = mPad->fCathode;   // chamber number
 126                         Int_t nhit     = mPad->fHitNumber; // hit number
 127                         Int_t qtot     = mPad->fQ;         // charge
 128                         Int_t ipx      = mPad->fPadX;      // pad number on X
 129                         Int_t ipy      = mPad->fPadY;      // pad number on Y
 130                         Int_t iqpad    = mPad->fQpad;      // charge per pad
 131                         Int_t rpad     = mPad->fRpad;      // r-pos of pad
 132 //
 133 //
 134         if (cathode != (icat+1)) continue;
 135         segmentation=iChamber->GetSegmentationModel(cathode);
 136         //        Int_t Npx[nch-1]  = segmentation->Npx();
 137         //        Int_t Npy[nch-1]  = segmentation->Npy();
 138
 139         if (cathode==1) {
 140              Int_t npx=Npx[nch-1];
 141              Int_t npy=Npy[nch-1];
 142         } else {
 143              Int_t npx=Npy[nch-1];
 144              Int_t npy=Npx[nch-1];
 145         }
 146
 147         //         printf("nch, cathode, rmin, rmax, npx, npy %d %d %d %d %d %d\n",
 148         //                        nch,cathode,rmin,rmax,npx,npy);
 149
 150         //             printf("icat, iqpad, ipx, ipy %d %d %d %d \n",
 151         //                               icat, iqpad, ipx, ipy);
 152
 153
 154         // check boundaries
 155         if (rpad < rmin[nch-1] || iqpad ==0 || rpad > rmax[nch-1]) continue;
 156         //         if (rpad < rmin[nch-1] || iqpad ==0 || rpad > 81.65) continue;
 157          //         if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
 158
 159           // fill the info array
 160           TVector *trinfo;
 161            trinfo=new TVector(2);
 162
 163           trinfo(0)=(Float_t)track;
 164           trinfo(1)=(Float_t)iqpad;
 165
 166           //          Int_t digits[3];
 167           digits[0]=ipx;
 168           digits[1]=ipy;
 169           digits[2]=iqpad;
 170
 171
 172           Int_t trk[50];
 173           Int_t chtrk[50];
 174
 175     // build the list of fired pads and update the info
 176     AliMUONlist *pdigit=elem[nch-1][ipy+npy][ipx+npx];
 177     if (pdigit==0) {
 178        obj->AddAtAndExpand(new AliMUONlist(rpad,digits),counter);
 179        counter++;
 180        //       Int_t nentrobj=obj->GetEntriesFast();
 181        Int_t last=obj->GetLast();
 182        pdigit=(AliMUONlist*)obj->At(last);
 183        elem[nch-1][ipy+npy][ipx+npx]=pdigit;
 184        //       Int_t q=pdigit->fSignal;
 185        // list of tracks
 186        TObjArray *trlist=(TObjArray*)pdigit->TrackList();
 187        trlist->Add(trinfo);
 188      } else {
 189        // update charge
 190        (*pdigit).fSignal+=iqpad;
 191        // update list of tracks
 192        TObjArray* trlist=(TObjArray*)pdigit->TrackList();
 193        Int_t nptracks=trlist->GetEntriesFast();
 194         for (Int_t tr=0;tr<nptracks;tr++) {
 195            TVector *ptrk=(TVector*)trlist->At(tr);
 196            Int_t entries=ptrk->GetNrows();
 197            TVector &vtrk = *ptrk;
 198            trk[tr]=(Int_t)vtrk(0);
 199            chtrk[tr]=(Int_t)vtrk(1);
 200            if (trk[tr]==track) {
 201                 chtrk[tr]+=iqpad;
 202                 trlist->RemoveAt(tr);
 203                 trinfo(0)=trk[tr];
 204                 trinfo(1)=chtrk[tr];
 205                 trlist->AddAt(trinfo,tr);
 206            } else {
 207                 trlist->Add(trinfo);
 208            }
 209         } //end loop over list of tracks for one pad
 210      }  //  end if pdigit
 211
 212
 213                    } //end loop over clust
 214
 215                } // hit loop
 216            } // if MUON
 217        } // track loop
 218
 219     Int_t tracks[10];
 220     Int_t charges[10];
 221     cout<<"start filling digits "<<endl;
 222     //    Int_t digits[3];
 223     const Int_t zero_supm = 6;
 224     const Int_t adc_satm = 1024;
 225     Int_t nd=0;
 226
 227
 228     // start filling the digits
 229  for (Int_t id=0;id<10;id++) {
 230              Int_t nd=0;
 231              //             Int_t npx=Npx[id];
 232              //             Int_t npy=Npy[id];
 233
 234              //      Int_t nx=2*npx;
 235              //      Int_t ny=2*npx;
 236    for (Int_t iy=0;iy<nymax;iy++) {
 237      for (Int_t ix=0;ix<nxmax;ix++) {
 238         AliMUONlist *address=elem[id][iy][ix];
 239         if (address==0) continue;
 240         // do zero-suppression and signal truncation
 241         Int_t rpad=address->fRpad;
 242         Int_t q=address->fSignal;
 243         if ( q <= zero_supm  || rpad > 55) continue;
 244         //        if ( q <= zero_supm ) continue;
 245         if ( q > adc_satm)  q=adc_satm;
 246         digits[0]=(*address).fPadX;
 247         digits[1]=(*address).fPadY;
 248         digits[2]=address->fSignal;
 249
 250         TObjArray* trlist=(TObjArray*)address->TrackList();
 251         Int_t nptracks=trlist->GetEntriesFast();
 252         // this should be changed to accomodate the real number of tracks
 253         if (nptracks > 10) {
 254             cout<<"Attention - nptracks > 3 "<<nptracks<<endl;
 255             nptracks=10;
 256         }
 257         for (Int_t tr=0;tr<nptracks;tr++) {
 258            TVector *pp=(TVector*)trlist->At(tr);
 259            Int_t entries=pp->GetNrows();
 260            TVector &v2 = *pp;
 261            tracks[tr]=(Int_t)v2(0);
 262            charges[tr]=(Int_t)v2(1);
 263         } //end loop over list of tracks for one pad
 264
 265         // fill digits
 266          MUON->AddDigits(id+1,tracks,charges,digits);
 267          nd++;
 268
 269      }  // end loop over ix
 270    }  // end loop over iy
 271     cout<<"I'm out of the loops over pads"<<endl;
 272     cout<<" id  nd  "<<id<<" "<<nd<<endl;
 273
 274     //    CathodeIndex[id+icat*10] = nd;     // update cathodes index
 275
 276
 277  } //loop over chambers
 278     cout<<"I'm out of the loops for digitisation"<<endl;
 279
 280
 281         gAlice->TreeD()->Fill();
 282         TTree *TD=gAlice->TreeD();
 283         int ndig=TD->GetEntries();
 284         cout<<"number of digits  "<<ndig<<endl;
 285
 286         int ndig1=(MUON->Dch1())->GetEntriesFast();
 287         cout<<"number of digits  1   "<<ndig1<<endl;
 288         int ndig2=(MUON->Dch2())->GetEntriesFast();
 289         cout<<"number of digits  2   "<<ndig2<<endl;
 290         int ndig3=(MUON->Dch3())->GetEntriesFast();
 291         cout<<"number of digits  3   "<<ndig3<<endl;
 292         int ndig4=(MUON->Dch4())->GetEntriesFast();
 293         cout<<"number of digits  4   "<<ndig4<<endl;
 294         int ndig5=(MUON->Dch5())->GetEntriesFast();
 295         cout<<"number of digits  5   "<<ndig5<<endl;
 296         int ndig6=(MUON->Dch6())->GetEntriesFast();
 297         cout<<"number of digits  6   "<<ndig6<<endl;
 298         int ndig7=(MUON->Dch7())->GetEntriesFast();
 299         cout<<"number of digits  7   "<<ndig7<<endl;
 300         int ndig8=(MUON->Dch8())->GetEntriesFast();
 301         cout<<"number of digits  8   "<<ndig8<<endl;
 302         int ndig9=(MUON->Dch9())->GetEntriesFast();
 303         cout<<"number of digits  9   "<<ndig9<<endl;
 304         int ndig10=(MUON->Dch10())->GetEntriesFast();
 305         cout<<"number of digits  10   "<<ndig10<<endl;
 306
 307         MUON->ResetDigits();
 308
 309
 310         //        char hname[30];
 311         //        sprintf(hname,"TreeD%d",nev);
 312         //        gAlice->TreeD()->Write(hname);
 313
 314     } //end loop over cathodes
 315
 316         char hname[30];
 317         sprintf(hname,"TreeD%d",nev);
 318         gAlice->TreeD()->Write(hname);
 319     //    MUON->ResetDigits();
 320     //    if (CathodeIndex) CathodeIndex->Clear();
 321
 322     //    file->Write();
 323
 324      file->ls();
 325      file->Map();
 326    } // event loop
 327
 328    //     delete [] CathodeIndex;
 329
 330      file->Close();
 331      cout<<"END  digitisation "<<endl;
 332
 333 }