/[PAMELA software]/PamelaLevel2/src/PamLevel2.cpp
ViewVC logotype

Diff of /PamelaLevel2/src/PamLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Fri Jun 16 16:43:55 2006 UTC revision 1.6 by pam-fi, Tue Oct 24 10:24:27 2006 UTC
# Line 24  PamTrack::PamTrack(TrkTrack* t, CaloTrkV Line 24  PamTrack::PamTrack(TrkTrack* t, CaloTrkV
24      tof_track  = this->ToFTrkVar::GetToFTrkVar();      tof_track  = this->ToFTrkVar::GetToFTrkVar();
25      if(t) *trk_track  = *t;      if(t) *trk_track  = *t;
26      if(c) *calo_track = *c;      if(c) *calo_track = *c;
27      if(o) *tof_track  = *o;      if(o) *tof_track  = *o;    
28  };  };
29    
30  //--------------------------------------  //--------------------------------------
# Line 35  PamTrack::PamTrack(TrkTrack* t, CaloTrkV Line 35  PamTrack::PamTrack(TrkTrack* t, CaloTrkV
35   * Constructor   * Constructor
36   */   */
37  PamLevel2::PamLevel2(){  PamLevel2::PamLevel2(){
38      trk_obj  = this->TrkLevel2::GetTrkLevel2();          
39      calo_obj = this->CaloLevel2::GetCaloLevel2();      trk_l1_obj  = this->TrkLevel1::GetTrkLevel1();      
40      tof_obj  = this->ToFLevel2::GetToFLevel2();      trk_obj     = this->TrkLevel2::GetTrkLevel2();
41      trig_obj = this->TrigLevel2::GetTrigLevel2();      calo_obj    = this->CaloLevel2::GetCaloLevel2();
42      s4_obj   = this->S4Level2::GetS4Level2();      tof_obj     = this->ToFLevel2::GetToFLevel2();
43      nd_obj   = this->NDLevel2::GetNDLevel2();      trig_obj    = this->TrigLevel2::GetTrigLevel2();
44      ac_obj   = this->AcLevel2::GetAcLevel2();      s4_obj      = this->S4Level2::GetS4Level2();
45        nd_obj      = this->NDLevel2::GetNDLevel2();
46        ac_obj      = this->AcLevel2::GetAcLevel2();
47        orb_obj     = this->OrbitalInfo::GetOrbitalInfo();
48    
49        run_obj     = new GL_RUN();
50        
51        sorted_tracks = new TRefArray();
52        
53        CAL = true;
54        TRK = true;
55        TRG = true;
56        TOF = true;
57        S4  = true;
58        ND  = true;
59        AC  = true;
60        ORB = true;
61        
62        TRK_L1 = false;
63        
64        Tout  = NULL;
65    };
66    /**
67     * Destructor
68     */
69    PamLevel2::~PamLevel2(){
70                    
71        TrkLevel1::Clear();
72        
73        TrkLevel2::Clear();
74        CaloLevel2::Clear();
75        ToFLevel2::Clear();
76        TrigLevel2::Clear();
77        S4Level2::Clear();
78        NDLevel2::Clear();
79        AcLevel2::Clear();
80        OrbitalInfo::Clear();
81        
82        delete run_obj;
83        
84    //      delete sorted_tracks;
85        sorted_tracks->Delete(); // clean the reference array
86            
87        if(Tout)Tout->Delete(); // delete loaded tree from memory
88            
89    };
90    /**
91     * Clear the event
92     */
93    void PamLevel2::Clear(){
94            
95        TrkLevel1::Clear();
96            
97        TrkLevel2::Clear();
98        CaloLevel2::Clear();
99        ToFLevel2::Clear();
100        TrigLevel2::Clear();
101        S4Level2::Clear();
102        NDLevel2::Clear();
103        AcLevel2::Clear();
104        OrbitalInfo::Clear();
105        
106        sorted_tracks->Delete(); // clean the reference array
107        
108  };  };
109    
110    
111    //--------------------------------------
112    //
113    //
114    //--------------------------------------
115    /**
116     * Retrieves the calorimeter track matching the seqno-th tracker stored track.
117     * (If seqno = -1 retrieves the self-trigger calorimeter track)
118     */
119     CaloTrkVar *PamLevel2::GetCaloStoredTrack(int seqno){
120            
121         if( CaloLevel2::ntrk()==0 ){
122             cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no Calorimeter tracks are stored"<<endl;
123             return NULL;
124         };
125        
126         CaloTrkVar *c = 0;
127         Int_t it_calo=0;
128        
129         do{
130             c = CaloLevel2::GetCaloTrkVar(it_calo);
131             it_calo++;
132         } while( c && seqno != c->trkseqno && it_calo < CaloLevel2::ntrk());      
133        
134         if(!c || seqno != c->trkseqno){
135             c = 0;
136             if(seqno!=-1)cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match Calorimeter stored tracks"<<endl;
137         };
138         return c;
139        
140     };
141    //--------------------------------------
142     //
143     //
144    //--------------------------------------
145    /**
146      * Retrieves the ToF track matching the seqno-th tracker stored track.
147      * (If seqno = -1 retrieves the tracker-independent tof track)
148     */
149     ToFTrkVar *PamLevel2::GetToFStoredTrack(int seqno){
150            
151         if( ToFLevel2::ntrk()==0 ){
152             cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no ToF tracks are stored"<<endl;
153             return NULL;
154         };
155        
156         ToFTrkVar *c = 0;
157         Int_t it_tof=0;
158        
159         do{
160             c = ToFLevel2::GetToFTrkVar(it_tof);
161             it_tof++;
162         } while( c && seqno != c->trkseqno && it_tof < ToFLevel2::ntrk());
163        
164         if(!c || seqno != c->trkseqno){
165             c = 0;
166             if(seqno!=-1)cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match ToF stored tracks"<<endl;
167         };
168         return c;
169        
170     };
171    
172    //--------------------------------------
173     //
174     //
175    //--------------------------------------
176    /**
177     * Give the pamela track associated to a tracker track, retrieving related calorimeter and tof track information.
178     */
179     PamTrack* PamLevel2::GetPamTrackAlong(TrkTrack* t){
180            
181         CaloTrkVar *c = 0;
182         ToFTrkVar  *o = 0;
183        
184         if(CAL) c = GetCaloStoredTrack(t->GetSeqNo());
185         if(TOF) o = GetToFStoredTrack(t->GetSeqNo());
186        
187    //    if(t && c && o)track = new PamTrack(t,c,o);
188         PamTrack *track = new PamTrack(t,c,o);
189        
190         return track;
191    
192     };
193  //--------------------------------------  //--------------------------------------
194  //  //
195  //  //
# Line 55  PamLevel2::PamLevel2(){ Line 202  PamLevel2::PamLevel2(){
202    
203  PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){  PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){
204            
205      // retrieve itrk-th tracker stored track      PamTrack *track = 0;
206      TrkTrack   *t = 0;      
207      CaloTrkVar *c = 0;      if( itrk >=0 && itrk < TrkLevel2::ntrk() ){
208      ToFTrkVar  *o = 0;          
209            TrkTrack *t = TrkLevel2::GetStoredTrack(itrk);
210      if( itrk >=0 && itrk < TrkLevel2::ntrk() ){          track = GetPamTrackAlong(t);
211          t = TrkLevel2::GetStoredTrack(itrk);              
212          c = CaloLevel2::GetCaloTrkVar(t->GetSeqNo());      }else{
213          o = ToFLevel2::GetToFTrkVar(t->GetSeqNo());          cout << "PamLevel2::GetStoredTrack(int) : tracker track SeqNo "<< itrk <<" does not exist (GetNTracks() = "<<TrkLevel2::GetNTracks()<<")"<<endl;
214      };      };
     //  retrieve related ToF-track  
215            
     // hence create a "PamTrack" object  
     PamTrack *track = 0;  
     if(t && c && o)track = new PamTrack(t,c,o);  
216      return track;      return track;
217            
218  }  }
# Line 78  PamTrack* PamLevel2::GetStoredTrack(Int_ Line 221  PamTrack* PamLevel2::GetStoredTrack(Int_
221  //  //
222  //--------------------------------------  //--------------------------------------
223  /**  /**
224   * Sort physical (tracker) tracks and stores them in a TObjectArray (of TrkTrack objects).   * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is  PamLevel2::sorted_tracks.
225   * The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it).   * @param how String to set the sorting cryterium (es: "CAL" or "TRK+CAL+TOF" ecc...).
226   * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.   * Sorting cryteria:
227   * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.   * TRK: lower chi**2
228     * CAL: lower Y spatial residual on the first calorimeter plane
229     * TOF: bigger numebr of hit PMTs along the track, on S12 S21 S32 (where paddles are along the Y axis).
230     * The default sorting cryterium is "TOF+CAL".
231     *
232     * The total number of physical tracks is always given by GetNTracks() and the it-th physical track can be retrieved by means of the methods GetTrack(int it) and GetTrack(int it, TString how).
233   */   */
234  TClonesArray *PamLevel2::GetTracks(){  void PamLevel2::SortTracks(TString how){
   
     // define new array for sorted tracker tracks  
     TClonesArray *aa = new TClonesArray("TrkTrack");  
     TClonesArray &sorted = *aa;  
235    
236      //Save current Object count
237        Int_t ObjectNumber = TProcessID::GetObjectCount();
238        
239        sorted_tracks->Delete(); //temporaneo???
240            
241      // loop over the tracks sorted by the tracker      // loop over the tracks sorted by the tracker
242        Bool_t use_TRK = how.Contains("TRK", TString::kIgnoreCase);
243        Bool_t use_CAL = how.Contains("CAL", TString::kIgnoreCase);
244        Bool_t use_TOF = how.Contains("TOF", TString::kIgnoreCase);
245        
246        if( !CAL &&  use_CAL) use_CAL = false;
247        if( !TOF &&  use_TOF) use_TOF = false;
248        
249        if( !TRK ){
250            cout << "SortTracks() : without tracker does not work!!! (not yet)" << endl;
251            return;
252        };
253        
254      for(Int_t i=0; i < TrkLevel2::GetNTracks(); i++){      for(Int_t i=0; i < TrkLevel2::GetNTracks(); i++){
255            
256          TrkTrack *ts=0;          TrkTrack *ts = 0;
257            
258          // get tracker tracks          // get tracker tracks
259          TrkTrack   *tp = TrkLevel2::GetTrack(i);                    //tracker          TrkTrack   *tp = TrkLevel2::GetTrack(i);                    //tracker
260          CaloTrkVar *cp = CaloLevel2::GetCaloTrkVar(tp->GetSeqNo()); //calorimeter  //              CaloTrkVar *cp = GetCaloStoredTrack(tp->GetSeqNo());
261          ToFTrkVar  *op = ToFLevel2::GetToFTrkVar(tp->GetSeqNo());   //tof  //              ToFTrkVar  *op = GetToFStoredTrack(tp->GetSeqNo());
262            CaloTrkVar *cp = 0;
263            ToFTrkVar  *op = 0;
264            
265          // if track has an image, check image selection          // if track has an image, check image selection
266          if(tp->HasImage()){          if(tp->HasImage()){
267                
268              TrkTrack   *ti = TrkLevel2::GetTrackImage(i);              //tracker (image)              TrkTrack   *ti = TrkLevel2::GetTrackImage(i);              //tracker (image)
269              CaloTrkVar *ci = CaloLevel2::GetCaloTrkVar(ti->GetSeqNo());//calorimeter (image)  //                      CaloTrkVar *ci = GetCaloStoredTrack(ti->GetSeqNo());
270              ToFTrkVar  *oi = ToFLevel2::GetToFTrkVar(ti->GetSeqNo());  //tof (image)  //                      ToFTrkVar  *oi = GetToFStoredTrack(ti->GetSeqNo());
271                CaloTrkVar *ci = 0;
272                ToFTrkVar  *oi = 0;
273                
274              //assign starting scores              //assign starting scores
275              Int_t tp_score = 1;              Int_t tp_score = 0;  //main track sorted by the tracker
276              Int_t ti_score = 0;              Int_t ti_score = 0;  //image track
277                
278              // ------------------------              // ------------------------
279              // calorimeter check              // calorimeter check
280              // ------------------------              // ------------------------
281                // check the Y spatial residual on the first calorimeter plane
282                // (cut on calorimeter variables from Emiliano)
283              if(              if(
284                        npcfit[1] > 3   &&  //no. of fit planes on Y view                  use_CAL            &&
285  //              varcfit[1] < 50.&&  //fit variance on Y view                  npcfit[1] > 15     &&   //no. of fit planes on Y view
286                        true){                  varcfit[1] < 1000. &&  //fit variance on Y view
287                    true){
288                  Float_t resy_p = cp->tbar[0][1] - cbar[0][1];if(resy_p < 0)resy_p= - resy_p;                  
289                  Float_t resy_i = ci->tbar[0][1] - cbar[0][1];if(resy_i < 0)resy_i= - resy_i;                  cp = GetCaloStoredTrack(tp->GetSeqNo());
290                    ci = GetCaloStoredTrack(ti->GetSeqNo());
291                    
292                    Float_t resy_p = cp->tbar[0][1] - cbar[0][1]; if(resy_p < 0)resy_p= - resy_p;
293                    Float_t resy_i = ci->tbar[0][1] - cbar[0][1]; if(resy_i < 0)resy_i= - resy_i;
294                    
295                  if(resy_p <= resy_i) tp_score++;                  if(resy_p <= resy_i) tp_score++;
296                  else                 ti_score++;                  else                 ti_score++;
297                        };              };
298              // ------------------------              // ------------------------
299              // TOF check              // TOF check
300              // ------------------------                  // ------------------------    
301                // check the number of hit pmts along the track
302                // on S12 S21 and S32, where paddles are parallel to Y axis
303                if( use_TOF ){
304                    
305                    Int_t nphit_p =0;
306                    Int_t nphit_i =0;
307                    
308                    op = GetToFStoredTrack(tp->GetSeqNo());
309                    oi = GetToFStoredTrack(ti->GetSeqNo());
310                    
311    /*                              cout << "track: npmtadc "<< op->npmtadc << endl;
312                                    cout << "track: npmttdc "<< op->npmttdc << endl;
313                                    cout << "image: npmtadc "<< oi->npmtadc << endl;
314                                    cout << "image: npmttdc "<< oi->npmttdc << endl;*/
315                    
316                    for (Int_t ih=0; ih < op->npmtadc; ih++){
317                        Int_t pl = GetPlaneIndex( (op->pmtadc).At(ih) );
318                        if(pl == 1 || pl == 2 || pl == 5)nphit_p++;
319                    };
320                    
321                    for (Int_t ih=0; ih < oi->npmtadc; ih++){
322                        Int_t pl = GetPlaneIndex( (oi->pmtadc).At(ih) );
323                        if(pl == 1 || pl == 2 || pl == 5)nphit_i++;
324                    };
325                    
326                    if(
327                        use_TOF            &&
328                        (nphit_p+nphit_i) !=0 &&    
329                        true){
330                        
331                        if( nphit_p >= nphit_i) tp_score++;
332                        else ti_score++;
333                    };
334                };
335                if(tp_score == ti_score) use_TRK = true;
336                // ------------------------
337                // tracker check
338                // ------------------------
339                // chi**2 difference is not always large enough to distinguish among
340                // the real track and its image.
341                // Tracker check will be applied always when calorimeter and tof information is ambiguous.
342                if(use_TRK){
343                    if(      tp->chi2 > 0 && tp->chi2 < ti->chi2 ) tp_score++ ;
344                    else if( ti->chi2 > 0 && ti->chi2 < tp->chi2 ) ti_score++ ;
345                };
346                
347              // ------------------------              // ------------------------
348              // the winner is....              // the winner is....
349              // ------------------------                  // ------------------------    
350                        if(tp_score > ti_score) ts = tp;//the track sorted by the tracker!!              if      (tp_score > ti_score) ts = tp;//the track sorted by the tracker!!
351                        else                    ts = ti;//its image!!              else if (tp_score < ti_score) ts = ti;//its image!!
352                else {
353                    ts = tp;
354    //                              cout << "Warning - track image ambiguity not solved" << endl;
355    //                              cout << ts->GetNtot() << " "<< ts->chi2 << " " << npcfit[1] << " "<< nphit_p << endl;
356                };
357                
358          }else{          }else{
359              ts = tp;              ts = tp;
360          };          };
361                    
362          new(sorted[i]) TrkTrack(*ts); //save the track in the sorted array  //              cout <<" SortTracks() "<<i<<" -- "<<ts<<endl;
363            sorted_tracks->Add(ts);//save the track in the sorted array
364    //              cout << "SortTracks:: sorted_tracks->Add(it) "<<i<<" "<<ts<<endl;
365            
366      };      };
367        //Restore Object count
368      return aa;      //To save space in the table keeping track of all referenced objects
369        //we assume that our events do not address each other. We reset the
370        //object count to what it was at the beginning of the event.
371        TProcessID::SetObjectCount(ObjectNumber);
372        
373  };  };
374  //--------------------------------------  //--------------------------------------
375  //  //
376  //  //
377  //--------------------------------------  //--------------------------------------
378  /**  /**
379     * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.
380     * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.
381     */
382     TRefArray *PamLevel2::GetTracks(){
383             //
384    //      SortTracks("+CAL+TOF");
385             SortTracks("+CAL+TRK");//TEMP!!!!!!!!!!!!!
386            //
387            return sorted_tracks;
388     };
389    //--------------------------------------
390     //
391     //
392    //--------------------------------------
393    /**
394   * Retrieves the it-th Pamela "physical" track.   * Retrieves the it-th Pamela "physical" track.
395   * It override TrkLevel2::GetTrack(int it).   * It override TrkLevel2::GetTrack(int it).
396   * @param it Track number, ranging from 0 to GetNTracks().   * @param it Track number, ranging from 0 to GetNTracks().
397   */   */
398  PamTrack *PamLevel2::GetTrack(int it){  PamTrack *PamLevel2::GetTrack(int it){
399            
400    //  *-*-*-*-*-*-*
401    //      SortTracks("+CAL+TOF");
402        SortTracks("+CAL+TRK");//TEMP!!!!!!!!!!!!!
403    //  *-*-*-*-*-*-*
404    
405      TClonesArray *aa = this->GetTracks();      PamTrack *track = 0;
406      TClonesArray &sorted = *aa;      
407        if( it >=0 && it < TrkLevel2::GetNTracks() && it<sorted_tracks->GetEntriesFast() ){
408      TrkTrack   *t = 0;          TrkTrack   *t = (TrkTrack*)sorted_tracks->At(it);
409      CaloTrkVar *c = 0;          track = GetPamTrackAlong(t);
410      ToFTrkVar  *o = 0;      }else{
411            cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<TrkLevel2::GetNTracks()<<")"<<endl;
     if( it >=0 && it < TrkLevel2::GetNTracks() ){  
         t = (TrkTrack*)sorted[it];  
         c = CaloLevel2::GetCaloTrkVar(t->GetSeqNo());  
         o = ToFLevel2::GetToFTrkVar(t->GetSeqNo());  
412      };      };
413            
     // hence create a "PamTrack" object  
     PamTrack *track = 0;  
     if(t && c && o)track = new PamTrack(t,c,o);  
414      return track;      return track;
415        
416  };  };
417    
418  //--------------------------------------  //--------------------------------------
419  //  //
420  //  //
421  //--------------------------------------  //--------------------------------------
422  /**  /**
423   * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::GetTracks().   * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::SortTracks().
424   * @param it Track number, ranging from 0 to GetNTracks().   * @param it Track number, ranging from 0 to GetNTracks().
425   */   */
426  PamTrack *PamLevel2::GetTrackImage(int it){  PamTrack *PamLevel2::GetTrackImage(int it){
427    
428      TClonesArray *aa = this->GetTracks();  //      SortTracks("+CAL+TOF");
429      TClonesArray &sorted = *aa;          SortTracks("+CAL+TRK");//TEMP!!!!!!!!!!!!!
430            
431            PamTrack *image = 0;
432    
433      TrkTrack   *t = 0;          if( it >=0 && it < TrkLevel2::GetNTracks() ){
434      CaloTrkVar *c = 0;                  TrkTrack *temp = (TrkTrack*)sorted_tracks->At(it);
435      ToFTrkVar  *o = 0;                  if( temp->HasImage() ){
436                            TrkTrack   *t = TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());
437      if( it >=0 && it < TrkLevel2::GetNTracks() ){                          image = GetPamTrackAlong(t);
438          TrkTrack *temp = (TrkTrack*)sorted[it];                  }else{
439          if( temp->HasImage() ){                          cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
440              t = TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());                  };
441              c = CaloLevel2::GetCaloTrkVar(temp->GetImageSeqNo());          }else{
442              o = ToFLevel2::GetToFTrkVar(temp->GetImageSeqNo());                  cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<TrkLevel2::GetNTracks()<<")"<<endl;
443          };          };
     };  
       
     // hence create a "PamTrack" object  
     PamTrack *image = 0;  
     if(t && c && o)image = new PamTrack(t,c,o);  
     return image;  
444            
445            return image;
446  }  }
447    
448  //--------------------------------------  //--------------------------------------
# Line 209  PamTrack *PamLevel2::GetTrackImage(int i Line 450  PamTrack *PamLevel2::GetTrackImage(int i
450  //  //
451  //--------------------------------------  //--------------------------------------
452  /**  /**
453   * Get the Pamela detector trees and make them friends.   * Get the Pamela detector trees in a single file and make them friends.
454     * @param f TFile pointer
455     * @param detlist String to select trees to be included
456     * @return Pointer to a TTree
457   */   */
458  TTree *PamLevel2::LoadPamTrees(TFile *f){  TTree *PamLevel2::GetPamTree(TFile *f, TString detlist="+ALL"){
459        
460      TTree *Tout =0;  //      cout << "WARNING!!! -- obsolete method -- \n";
461    //      cout << "better use TChain *PamLevel2::GetPamTree(TList*, TString) \n";
462            
463        TTree *Trout =0;
464            
465      // Tracker      SetWhichTrees(detlist);
466    // Tracker
467      TTree *T = (TTree*)f->Get("Tracker");      TTree *T = (TTree*)f->Get("Tracker");
468      if(T) {      if(T && TRK) {
469          T->SetBranchAddress("TrkLevel2", GetPointerToTrk());          if(TRK_L1){
470          cout << "Tracker      : set branch address TrkLevel2"<<endl;              T->SetBranchAddress("TrkLevel1", GetPointerToTrk(1));
471          if(!Tout)Tout=T;              cout << "Tracker      : set branch address TrkLevel1"<<endl;
472            };
473            T->SetBranchAddress("TrkLevel2", GetPointerToTrk(2));
474            cout << "Tracker      : set branch address TrkLevel2"<<endl;
475            Trout=T;
476      }else{      }else{
477          cout << "Tracker      : missing tree"<<endl;          cout << "Tracker      : missing tree"<<endl;
478      };      };
479      // Calorimeter      // Calorimeter
480      TTree *C = (TTree*)f->Get("Calorimeter");      TTree *C = (TTree*)f->Get("Calorimeter");
481      if(C) {      if(C && CAL) {
482          C->SetBranchAddress("CaloLevel2", GetPointerToCalo());          C->SetBranchAddress("CaloLevel2", GetPointerToCalo());
483          cout << "Calorimeter  : set branch address CaloLevel2"<<endl;          cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
484          if(!Tout)Tout=C;          if(!Trout)Trout=C;
485          else Tout->AddFriend(C);          else Trout->AddFriend(C);
486      }else{      }else{
487          cout << "Calorimeter  : missing tree"<<endl;          cout << "Calorimeter  : missing tree"<<endl;
488      };      };
489      // ToF          // ToF    
490      TTree *O = (TTree*)f->Get("ToF");      TTree *O = (TTree*)f->Get("ToF");
491      if(O) {      if(O && TOF) {
492          O->SetBranchAddress("ToFLevel2", GetPointerToToF());          O->SetBranchAddress("ToFLevel2", GetPointerToToF());
493          cout << "ToF          : set branch address ToFLevel2"<<endl;          cout << "ToF          : set branch address ToFLevel2"<<endl;
494          if(!Tout)Tout=O;          if(!Trout)Trout=O;
495          else Tout->AddFriend(O);          else Trout->AddFriend(O);
496      }else{      }else{
497          cout << "ToF         : missing tree"<<endl;          cout << "ToF         : missing tree"<<endl;
498      };      };
499      // Trigger      // Trigger
500      TTree *R = (TTree*)f->Get("Trigger");      TTree *R = (TTree*)f->Get("Trigger");
501      if(R) {      if(R && TRG) {
502          R->SetBranchAddress("TrigLevel2", GetPointerToTrig());          R->SetBranchAddress("TrigLevel2", GetPointerToTrig());
503          cout << "Trigger      : set branch address TrigLevel2"<<endl;          cout << "Trigger      : set branch address TrigLevel2"<<endl;
504          if(!Tout)Tout=O;          if(!Trout)Trout=O;
505          else Tout->AddFriend(R);          else Trout->AddFriend(R);
506      }else{      }else{
507          cout << "Trigger      : missing tree"<<endl;          cout << "Trigger      : missing tree"<<endl;
508      };      };
509      // S4      // S4
510      TTree *S = (TTree*)f->Get("S4");      TTree *S = (TTree*)f->Get("S4");
511      if(S) {      if(S && S4) {
512          S->SetBranchAddress("S4Level2", GetPointerToS4());          S->SetBranchAddress("S4Level2", GetPointerToS4());
513          cout << "S4           : set branch address S4Level2"<<endl;          cout << "S4           : set branch address S4Level2"<<endl;
514          if(!Tout)Tout=O;          if(!Trout)Trout=O;
515          else Tout->AddFriend(S);          else Trout->AddFriend(S);
516      }else{      }else{
517          cout << "S4           : missing tree"<<endl;          cout << "S4           : missing tree"<<endl;
518      };      };
519      // Neutron Detector      // Neutron Detector
520      TTree *N = (TTree*)f->Get("NeutronD");      TTree *N = (TTree*)f->Get("NeutronD");
521      if(N) {      if(N && ND) {
522          N->SetBranchAddress("NDLevel2", GetPointerToND());          N->SetBranchAddress("NDLevel2", GetPointerToND());
523          cout << "NeutronD     : set branch address NDLevel2"<<endl;          cout << "NeutronD     : set branch address NDLevel2"<<endl;
524          if(!Tout)Tout=O;          if(!Trout)Trout=O;
525          else Tout->AddFriend(N);          else Trout->AddFriend(N);
526      }else{      }else{
527          cout << "NeutronD     : missing tree"<<endl;          cout << "NeutronD     : missing tree"<<endl;
528      };      };
529      // Anticounters      // Anticounters
530      TTree *A = (TTree*)f->Get("Anticounter");      TTree *A = (TTree*)f->Get("Anticounter");
531      if(A) {      if(A && AC) {
532          A->SetBranchAddress("AcLevel2", GetPointerToAc());          A->SetBranchAddress("AcLevel2", GetPointerToAc());
533          cout << "Anticounter  : set branch address AcLevel2"<<endl;          cout << "Anticounter  : set branch address AcLevel2"<<endl;
534          if(!Tout)Tout=O;          if(!Trout)Trout=O;
535          else Tout->AddFriend(A);          else Trout->AddFriend(A);
536      }else{      }else{
537          cout << "Anticounter  : missing tree"<<endl;          cout << "Anticounter  : missing tree"<<endl;
538      };      };
539        // Orbital Info
540            TTree *B = (TTree*)f->Get("OrbitalInfo");
541            if(B && ORB) {
542                    B->SetBranchAddress("OrbitalInfo", GetPointerToOrb());
543                    cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
544                    if(!Trout)Trout=O;
545                    else Trout->AddFriend(B);
546            }else{
547                    cout << "OrbitalInfo  : missing tree"<<endl;
548            };
549        
550            cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
551            
552            return Trout;
553        
554    }
555    //--------------------------------------
556    //
557    //
558    //--------------------------------------
559    /**
560     * Get list of Level2 files.
561     * @param ddir Level2 data directory.
562     * @param flisttxt Name of txt file containing file list.
563     * @return Pointer to a TList of TSystemFiles
564     * If no input file list is given , all the Level2 files inside the directory are processed.
565     */
566    TList*  PamLevel2::GetListOfLevel2Files(TString ddir, TString flisttxt = ""){
567            
568        TString wdir = gSystem->WorkingDirectory();
569        
570        if(ddir=="")ddir = wdir;
571    //      TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
572        cout << "Level2 data directory : "<<  ddir << endl;
573        
574        TList *contents  = new TList; // create output list
575        contents->SetOwner();
576        
577        char *fullpath;
578        
579        // if no input file list is given:  
580        if ( flisttxt != "" ){
581            
582            if( !gSystem->IsFileInIncludePath(flisttxt,&fullpath) )return 0;
583            
584    //              flisttxt = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
585            flisttxt = fullpath;
586            
587            if( !gSystem->ChangeDirectory(ddir) )return 0;
588            
589            cout <<"Input file list : " << flisttxt <<endl;
590            ifstream in;
591            in.open(flisttxt, ios::in); //open input file list
592            while (1) {
593                TString file;
594                in >> file;
595                if (!in.good()) break;
596                if( gSystem->IsFileInIncludePath(file,&fullpath) ){
597    //                              contents->Add(new TSystemDirectory(fullpath,ddir)); // add file to the list
598    //                              contents->Add(new TSystemFile(fullpath,ddir)); // add file to the list
599                    contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
600                }else{
601                    cout << "warning! --- File "<<file<<" does not exists"<< endl;
602                };
603            };
604            in.close();
605            
606        }else{
607            
608            cout << "No input file list given."<<endl;
609            cout << "Check for existing root files."<<endl;
610    //              cout << "Warking directory: "<< gSystem->WorkingDirectory()<< endl;
611            
612            TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
613            TList *temp = datadir->GetListOfFiles();
614    //              temp->Print();
615    //              cout << "*************************************" << endl;
616            
617            TIter next(temp);
618            TSystemFile *questo = 0;
619            
620            
621            while ( (questo = (TSystemFile*) next()) ) {
622                TString name =  questo-> GetName();
623                if( name.EndsWith(".root") ){
624                    char *fullpath;
625                    gSystem->IsFileInIncludePath(name,&fullpath);
626                    contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));
627                };
628            }
629            delete temp;
630            delete datadir;
631            
632        };
633        gSystem->ChangeDirectory(wdir); // back to the working directory
634    //      cout << endl << "Selected files:" << endl;
635    //      contents->Print();
636        cout << contents->GetEntries()<<" files selected\n";
637    //      cout << endl;
638    //      cout << "Working directory: "<< gSystem->WorkingDirectory()<< endl;
639        return contents;
640    };
641    //--------------------------------------
642    //
643    //
644    //--------------------------------------
645    /**
646     * Get the Pamela detector chains from a list of files and make them friends.
647     * @param fl Pointer to a TList of TSystemFiles
648     * @param detlist String to select trees to be included
649     * @return Pointer to a TChain
650     */
651    TChain *PamLevel2::GetPamTree(TList *fl, TString detlist="+ALL"){
652    
653    //      TChain *Tout=0;
654        if(Tout)Tout->Delete();
655        Tout = NULL;
656        
657        SetWhichTrees(detlist);
658        
659        TChain *T = 0;      
660        TChain *C = 0;
661        TChain *O = 0;
662        TChain *R = 0;
663        TChain *S = 0;
664        TChain *N = 0;
665        TChain *A = 0;
666        TChain *B = 0;
667        
668        if(TRK) T = new TChain("Tracker");  
669        if(CAL) C = new TChain("Calorimeter");
670        if(TOF) O = new TChain("ToF");
671        if(TRG) R = new TChain("Trigger");
672        if(S4)  S = new TChain("S4");
673        if(ND)  N = new TChain("NeutronD");
674        if(AC)  A = new TChain("Anticounter");
675        if(ORB) B = new TChain("OrbitalInfo");
676        
677        // loop over files and create chains        
678        TIter next(fl);
679        TSystemFile *questo = 0;
680        while ( (questo = (TSystemFile*) next()) ) {
681            TString name =  questo->GetName();
682    //              cout << "File: "<< name << endl;
683            if( CheckLevel2File(name) ){
684                if(TRK) T->Add(name);
685                if(CAL) C->Add(name);
686                if(TOF) O->Add(name);
687                if(TRG) R->Add(name);
688                if(S4)  S->Add(name);
689                if(ND)  N->Add(name);
690                if(AC)  A->Add(name);
691                if(ORB) B->Add(name);
692            };
693        }
694        
695        cout << "done chain \n";
696    
697        // Tracker
698        if(TRK) if(!Tout)Tout=T;
699        // Calorimeter
700        if(CAL) {          
701            if(!Tout)Tout=C;
702            else Tout->AddFriend("Calorimeter");
703        };
704        // ToF    
705        if(TOF) {
706            if(!Tout)Tout=O;
707            else Tout->AddFriend("ToF");
708        };
709        // Trigger
710        if(TRG) {
711            if(!Tout)Tout=O;
712            else Tout->AddFriend("Trigger");
713        };
714        // S4
715        if(S4) {
716            if(!Tout)Tout=O;
717            else Tout->AddFriend("S4");
718        };
719        // Neutron Detector
720        if(ND) {
721            if(!Tout)Tout=O;
722            else Tout->AddFriend("NeutronD");
723        };
724        // Anticounters
725        if(AC) {
726            if(!Tout)Tout=O;
727            else Tout->AddFriend("Anticounter");
728        };
729        // OrbitalInfo
730        if(ORB) {
731            if(!Tout)Tout=O;
732            else Tout->AddFriend("OrbitalInfo");
733        };
734        
735    //    cout<<endl<<" Number of entries: "<<Tout->GetEntries()<<endl<<endl;
736        
737        if( Tout->GetEntries() )PamLevel2::SetBranchAddress();
738            
739      return Tout;      return Tout;
740    }
741    //--------------------------------------
742    //
743    //
744    //--------------------------------------
745    /**
746     * Set branch addresses for Pamela friend trees
747     */
748    void PamLevel2::SetBranchAddress(){
749        // Tracker
750        if(TRK) {
751            if(TRK_L1){
752                Tout->SetBranchAddress("TrkLevel1", GetPointerToTrk(1));
753                cout << "Tracker      : set branch address TrkLevel1"<<endl;
754            };
755            Tout->SetBranchAddress("TrkLevel2", GetPointerToTrk(2));
756            cout << "Tracker      : set branch address TrkLevel2"<<endl;
757        };
758        
759        // Calorimeter
760        if(CAL) {
761            Tout->SetBranchAddress("CaloLevel2", GetPointerToCalo());
762            cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
763        };
764        
765        // ToF    
766        if(TOF) {
767            Tout->SetBranchAddress("ToFLevel2", GetPointerToToF());
768            cout << "ToF          : set branch address ToFLevel2"<<endl;
769        };
770        // Trigger
771        if(TRG) {
772            Tout->SetBranchAddress("TrigLevel2", GetPointerToTrig());
773            cout << "Trigger      : set branch address TrigLevel2"<<endl;
774        };
775        // S4
776        if(S4) {
777            Tout->SetBranchAddress("S4Level2", GetPointerToS4());
778            cout << "S4           : set branch address S4Level2"<<endl;
779        };
780        // Neutron Detector
781        if(ND) {
782            Tout->SetBranchAddress("NDLevel2", GetPointerToND());
783            cout << "NeutronD     : set branch address NDLevel2"<<endl;
784        };
785        // Anticounters
786        if(AC) {
787            Tout->SetBranchAddress("AcLevel2", GetPointerToAc());
788            cout << "Anticounter  : set branch address AcLevel2"<<endl;
789        };
790        // OrbitalInfo
791        if(ORB) {
792            Tout->SetBranchAddress("OrbitalInfo", GetPointerToOrb());
793            cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
794        };
795        
796    }
797    
798    //--------------------------------------
799    //
800    //
801    //--------------------------------------
802    void *PamLevel2::GetPointerTo(const char* c ){
803    //      cout << "objname "<< objname << endl;
804        TString objname = c;
805        if(!objname.CompareTo("TrkLevel1"))return &trk_l1_obj;
806        if(!objname.CompareTo("TrkLevel2"))return &trk_obj;
807        
808        return NULL;
809    };
810    //--------------------------------------
811    //
812    //
813    //--------------------------------------
814    /**
815     * Get the Run tree chain from a list of files.
816     * @param fl Pointer to a TList of TSystemFiles
817     * @return Pointer to a TChain
818     */
819    TChain *PamLevel2::GetRunTree(TList *fl){
820            
821        TChain *R = new TChain("Run");
822        
823        // loop over files and create chains        
824        TIter next(fl);
825        TSystemFile *questo = 0;
826        while ( (questo = (TSystemFile*) next()) ) {
827            TString name =  questo->GetName();
828    //              cout << "File: "<< name << endl;
829            if( CheckLevel2File(name) ){
830                R->Add(name);
831            };
832        }
833        
834        R->SetBranchAddress("RunInfo", GetPointerToRun());
835        cout << "Run         : set branch address RunInfo"<<endl;
836        
837        return R;
838        
839    }
840    //--------------------------------------
841    //
842    //
843    //--------------------------------------
844    /**
845     * Get the Run tree  from a file.
846     * @param f Pointer to a TFile
847     * @return Pointer to a TTree
848     */
849    TTree *PamLevel2::GetRunTree(TFile *f){
850    
851    
852        TTree *R = (TTree*)f->Get("Run");
853        
854        R->SetBranchAddress("RunInfo", GetPointerToRun());
855        cout << "Run         : set branch address RunInfo"<<endl;
856        
857        return R;
858            
859  }  }
860    //--------------------------------------
861    //
862    //
863    //--------------------------------------
864    /**
865     * Set which trees should be analysed
866     * @param detlist TString containing the sequence of trees required
867    */
868    void PamLevel2::SetWhichTrees(TString detlist){
869            
870        if(detlist.Contains("+ALL", TString::kIgnoreCase)){
871            CAL = true;
872            TRK = true;
873            TRK_L1 = false;
874            TRG = true;
875            TOF = true;
876            S4  = true;
877            ND  = true;
878            AC  = true;
879            ORB = true;
880        }else if( detlist.Contains("-ALL", TString::kIgnoreCase) ){
881            CAL = false;
882            TRK = false;
883            TRK_L1 = false;
884            TRG = false;
885            TOF = false;
886            S4  = false;
887            ND  = false;
888            AC  = false;
889            ORB = false;
890        };
891        
892        if( detlist.Contains("-CAL", TString::kIgnoreCase) )CAL = false;
893        else if( detlist.Contains("+CAL", TString::kIgnoreCase) )CAL = true;
894        
895        if( detlist.Contains("-TRK", TString::kIgnoreCase) )TRK = false;
896        else if( detlist.Contains("+TRK", TString::kIgnoreCase) )TRK = true;
897        
898        if( detlist.Contains("-TRK1", TString::kIgnoreCase) )TRK_L1 = false;
899        else if( detlist.Contains("+TRK1", TString::kIgnoreCase) )TRK_L1 = true;
900        
901        if( detlist.Contains("-TRG", TString::kIgnoreCase) )TRG = false;
902        else if( detlist.Contains("+TRG", TString::kIgnoreCase) )TRG = true;
903        
904        if( detlist.Contains("-TOF", TString::kIgnoreCase) )TOF = false;
905        else if( detlist.Contains("+TOF", TString::kIgnoreCase) )TOF = true;
906        
907        if( detlist.Contains("-S4",  TString::kIgnoreCase) )S4  = false;
908        else if( detlist.Contains("+S4",  TString::kIgnoreCase) )S4  = true;
909        
910        if( detlist.Contains("-ND",  TString::kIgnoreCase) )ND  = false;
911        else if( detlist.Contains("+ND",  TString::kIgnoreCase) )ND  = true;
912        
913        if( detlist.Contains("-AC",  TString::kIgnoreCase) )AC  = false;
914        else if( detlist.Contains("+AC",  TString::kIgnoreCase) )AC  = true;
915        
916        if( detlist.Contains("-ORB", TString::kIgnoreCase) )ORB = false;
917        else if( detlist.Contains("+ORB", TString::kIgnoreCase) )ORB = true;
918        
919    };
920    //--------------------------------------
921    //
922    //
923    //--------------------------------------
924    /**
925     * Check if a file contains selected Pamela Level2 trees.
926     * @param name File name
927     * @return true if the file is ok.
928     */
929    Bool_t  PamLevel2::CheckLevel2File(TString name){
930            
931        Bool_t CAL__ok = false;    
932        Bool_t TRK__ok = false;    
933        Bool_t TRK_L1__ok = false;    
934        Bool_t TRG__ok = false;    
935        Bool_t TOF__ok = false;    
936        Bool_t S4__ok  = false;    
937        Bool_t ND__ok  = false;    
938        Bool_t AC__ok  = false;    
939        Bool_t ORB__ok = false;    
940        
941        Bool_t RUN__ok = false;
942        
943    //    cout << "Checking file: "<<name<<endl;
944        TFile *f = new TFile(name.Data());
945        if( !f || f->IsZombie() ){
946            cout << "File: "<< f->GetName() <<" discarded ---- Non valid root file"<< endl; return false;
947        }
948        //   cout << "Get list of keys: "<<f<<endl;
949        TList *lk = f->GetListOfKeys();
950    //      lk->Print();
951        TIter next(lk);
952        TKey *key =0;
953        Int_t nev =0;
954        Int_t nev_previous =0;
955        while( (key = (TKey*)next()) ){
956            
957            nev_previous = nev;
958            nev = ((TTree*)f->Get(key->GetName()))->GetEntries();
959            if( key->GetClassName()=="TTree" && nev_previous && nev != nev_previous ){
960                cout << "File: "<< f->GetName() <<" discarded ---- "<< key->GetName() << " tree: n.entries does not match "<<nev<<" "<<nev_previous<< endl;
961                return false;
962            };
963    
964            if( !strcmp(key->GetName(),"Calorimeter") )CAL__ok = true;
965            if( !strcmp(key->GetName(),"Trigger"    ) )TRG__ok = true;
966            if( !strcmp(key->GetName(),"ToF"        ) )TOF__ok = true;
967            if( !strcmp(key->GetName(),"S4"         ) )S4__ok = true;
968            if( !strcmp(key->GetName(),"NeutronD"   ) )ND__ok = true;
969            if( !strcmp(key->GetName(),"Anticounter") )AC__ok = true;
970            if( !strcmp(key->GetName(),"OrbitalInfo") )ORB__ok = true;
971            if( !strcmp(key->GetName(),"Run"        ) )RUN__ok = true;      
972            if( !strcmp(key->GetName(),"Tracker"    ) ){
973                TRK__ok = true;
974                TTree *T = (TTree*)f->Get("Tracker");
975                for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
976                    TString name = T->GetListOfBranches()->At(i)->GetName();
977                    if( !name.CompareTo("TrkLevel1") )TRK_L1__ok=true;
978                };
979            };
980            
981        };
982        
983        lk->Delete();
984        f->Close();
985        
986        Bool_t FLAG = true;
987        if(!RUN__ok) {
988            cout << "File: "<< f->GetName() <<" discarded ---- Missing RunInfo tree"<< endl; return false;
989    //      FLAG = false;
990        };
991        if(CAL && !CAL__ok){
992            cout << "File: "<< f->GetName() <<" discarded ---- Missing Calorimeter tree"<< endl; return false;
993    //      FLAG = false;
994        };
995        if(TRK && !TRK__ok){
996            cout << "File: "<< f->GetName() <<" discarded ---- Missing Tracker tree"<< endl; return false;
997    //      FLAG = false;
998        };
999        if(TRK_L1 && !TRK_L1__ok){
1000            cout << "File: "<< f->GetName() <<" discarded ---- Missing Tracker Level1 Branch"<< endl; return false;
1001    //      FLAG = false;
1002        };
1003        if(TRG && !TRG__ok){
1004            cout << "File: "<< f->GetName() <<" discarded ---- Missing Trigger tree"<< endl; return false;
1005    //      FLAG = false;
1006        };
1007        if(TOF && !TOF__ok){
1008            cout << "File: "<< f->GetName() <<" discarded ---- Missing ToF tree"<< endl; return false;
1009    //      FLAG = false;
1010        };
1011        if(S4 && !S4__ok){
1012            cout << "File: "<< f->GetName() <<" discarded ---- Missing S4 tree"<< endl; return false;
1013    //      FLAG = false;
1014        };
1015        if(ND && !ND__ok){
1016            cout << "File: "<< f->GetName() <<" discarded ---- Missing ND tree"<< endl; return false;
1017    //      FLAG = false;
1018        };
1019        if(AC && !AC__ok){
1020            cout << "File: "<< f->GetName() <<" discarded ---- Missing AC tree"<< endl; return false;
1021    //      FLAG = false;
1022        };
1023        if(ORB && !ORB__ok){
1024            cout << "File: "<< f->GetName() <<" discarded ---- Missing ORB tree"<< endl; return false;
1025    //      FLAG = false;
1026        };
1027        
1028        return FLAG;
1029            
1030    };
1031    
1032    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23